In [77]:

    

n = int(10000)
import pyprimes
p = list(pyprimes.primes_below(n/2))


    

In [78]:

    

from bisect import bisect_left
from collections import Counter
def divisors(m):
    counter = Counter()
    for d in p[:bisect_left(p, m/2)]:
        while not m%d:
            counter[d] += 1
            m = m/d
    return(counter)


    

In [102]:

    

from numpy import prod
from itertools import combinations
from itertools import chain
def d(m):
    d = divisors(m)
    return sum(set(chain.from_iterable( (prod(i) for i in combinations(d.elements(), k)) for k in range(1, sum(d.values())))))


    

In [104]:

    

da = { i:d(i) for i in range(1,n) }


    

In [105]:

    

s = 0
for k , v in da.items():
    if v in da and k == da[v] and k != v:
        print(k, v)
        s += k
print(s)


    

    




220 284
284 220
1184 1210
1210 1184
2620 2924
2924 2620
5020 5564
5564 5020
6232 6368
6368 6232
31626

In [84]:

    

da[1184]


    

    
Out[84]:





118

In [103]:

    

d(1184)


    

    
Out[103]:





1210

In [ ]: