Der einfachste Ansatz ist

>>> def computeWeirdPairs(limit):
...     return [(a, b) for a in range(1,limit) for b in range(1,limit)
...             if not a%b]
...
>>> computeWeirdPairs(10)
[(1, 1), (2, 1), (2, 2), (3, 1), (3, 3), (4, 1), (4, 2),
(4, 4), (5, 1), (5, 5), (6, 1), (6, 2), (6, 3), (6, 6),
(7, 1), (7, 7), (8, 1), (8, 2), (8, 4), (8, 8), (9, 1),
(9, 3), (9, 9)]

Allerdings ist das ein quadratischer Algorithmus, d.h. er braucht bei tausendfacher Datenmenge eine Million Mal länger. Fällt euch eine schnellere Lösung ein?