How to make the following python code on finding amicable pairs more efficient

Question

I have written a python program which finds amicable pairs in a specific range. I recon there are plenty of better ways to do this and am looking for some feedback on how to improve my code.

def d(n):
    x = []
    for i in range(1, n):
        if n % i == 0:
            x.append(i)
    return sum(x)

def amicable(z,y):
    if d(z) == y and d(y) == z:
        print(z, y) 

for z in range(0, 10000, 2):
    for y in range(0, 10000, 2):
        if z != y:
            amicable(z, y)

This code actually does what its supposed to, but its not very efficient. I have to wait awhile for the results.

Do you mind to provide an [mcve](/help/mcve)? In particular I'd like to know which one are the inputs for `amicable`. — rpanai, Jun 27 '19 at 15:18
the inputs are z in range(0, 10000, 2) and y in range(0, 10000, 2). The results i get are the amicable pairs i.e (220, 284), (284, 220), (1184, 1210), (1210, 1184) etc — , Jun 27 '19 at 15:21
It's still not clear to me you are defining `a=[]` and then not using it. — rpanai, Jun 27 '19 at 15:31
Disregard the a=[]. I was fiddling around and forgot to remove it. It doesn't affect the rest of the code. — , Jun 27 '19 at 15:32
A bit of numpy magic, you can replace you `d(n)` with this line `d = (np.argwhere((n % np.arange(1,n))==0)+1).sum()` numpy is very efficient, for my small test (n=200) this produced the exact same vector — Rotem Tal, Jun 27 '19 at 15:50

score 0 · Answer 1 · answered Jun 27 '19 at 15:39

Unfortunately, my box with all of my project Euclid stuff is in storage, so I'll mock something up for you instead. I'll test this when I get a chance to compute some primes.

In general, the optimal way to find the divisors of a number is to only check primes. Also, you only need to check up through the square root of a number. Also, it is convenient to store a list of prime numbers, since they take a long time to compute.

You can also half the number of pairs you have to check with a little nested-loop magic.

Here's a general scheme:

def save_primes(primes):
    #save your primes file, they take a long time to compute
    pass
def load_primes(fname):
    #load your primes file
    ....
    return primes

def compute_primes(n):
    #get the prime numbers up to n
    primes = []
    ...
    save_primes(primes, "my_primes")

def get_prime_divisors(n):
    #return prime divisors and their multiplicity
    global PRIMES
    prime_divisors = []
    prime_divisor_multiplicities = []
    j = 0
    if "PRIMES" not in globals():
        PRIMES = load_primes("my_primes")
    while PRIMES[j] < n**0.5:
        prime = PRIMES[j]
        if n % prime == 0:                
            k = 0
            while n % prime == 0:
                n = n//prime
                k += 1
            prime_divisors.append(prime)
            prime_divisor_multiplicities.append(k)
    return prime_divisors, prime_divisor_multiplicities

def d(n):
    import itertools
    prime_divisors, multiplicities = get_prime_divisors(n)
    sum_divisors = 0
    for powers in itertools.product(map(range, multiplicities))):
        for i, prime in enumerate(prime_divisors):
            sum_divisors += prime**powers[i]
    sum_divisors -= n #only proper divisors
    return sum_divisors

def amicable(z,y):
    if d(z) == y and d(y) == z:
        return True
    return False

for z in range(0, 10000-1, 2):
    for y in range(z+1, 10000, 2):
        if amicable(z, y):
            print(z,y)

If you produce a prime factorization, you can reproduce all of the divisors very quickly. His question is related to speeding up the process, and this is one speed up. — Him, Jun 27 '19 at 16:00
see this [answer](https://stackoverflow.com/questions/38094818/what-is-the-most-efficient-way-to-find-amicable-numbers-in-python) — rpanai, Jun 27 '19 at 16:02
[Here](https://stackoverflow.com/questions/3939660/sieve-of-eratosthenes-finding-primes-python) is an efficient implementation of a prime-computing engine. — Him, Jun 27 '19 at 16:02
@rpanai The number of primes less than or equal to `n` [is actually sublinear](https://en.wikipedia.org/wiki/Prime_number#Definition_and_examples), so that answer is not optimal. — Him, Jun 27 '19 at 16:05
I suspect that amicable numbers need to have a large quantity of very small divisors, so there is probably a pruning algorithm to optimize for that.... I'll think about it and get back to you. :) — Him, Jun 27 '19 at 16:07

score 0 · Answer 2 · answered Jun 27 '19 at 15:53

I could speed up a little your code using numpy and avoiding to compute twice the same pair

import numpy as np

def divisors_sum(n):
    x = np.arange(1, n)
    return np.sum(x[(n % x) == 0])

def is_amicable(z,y):
    if divisors_sum(z) == y and divisors_sum(y) == z:
        return True
    return False

def amicable_pairs(N):
    for z in range(0, N, 2):
        for y in range(z+2, N, 2):
            if is_amicable(z, y) is True:
                print(z,y)

If I compare with your code

def d(n):
    x = []
    for i in range(1, n):
        if n % i == 0:
            x.append(i)
    return sum(x)


def is_amicable_old(z,y):
    if d(z) == y and d(y) == z:
        return True
    return False

def amicable_pairs_old(N):
    for z in range(0, N, 2):
        for y in range(0, N, 2):
            if z!=y:
                if is_amicable_old(z, y) is True:
                    print(z, y)

I got

%%time
amicable_pairs(1000)

220 284
CPU times: user 1.69 s, sys: 28 ms, total: 1.72 s
Wall time: 1.69 s

while

%%time
amicable_pairs_old(1000)

220 284
284 220
CPU times: user 7.35 s, sys: 7.85 ms, total: 7.35 s
Wall time: 7.35 s

How to make the following python code on finding amicable pairs more efficient

2 Answers2