4
\$\begingroup\$

Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n). If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.

My solution to the 21 problem on project Euler is very slow (with pure brute force it took 30 mins to find the solution). Any ways that I can improve my code or should I scrap it and think about something else ?

Here's my code :

import time, math

start = time.time()

amicable_nums = set()

def sum_div(n):
    divisors = []
    for x in range(1, int(math.sqrt(n) + 1)):
        if n % x == 0:
            divisors.append(x)

            if x * x != n and x != 1:
                divisors.append(int(n / x))
    return sum(divisors)

for i in range(1, 10000):
    for j in range(1, 10000):
        if sum_div(i) == j and sum_div(j) == i and i != j:
            amicable_nums.update([i, j])

print(sum(amicable_nums))
print("It took " + str(time.time() - start) + " seconds")
\$\endgroup\$
4
\$\begingroup\$

In your sum_div method, 1 is always a divisor, and n/1 is never a “proper” divisor of n. So you can remove the “1” case from your loop, and get rid of the x != 1 test, for a slight speed increase.

Also, if n % x == 0, then the other divisor is n // x; no call to int() is required.

Finally, no list creation is required. Just maintain a running total.

def sum_div(n):
    total = 1
    for x in range(2, int(math.sqrt(n) + 1)):
        if n % x == 0:
            total += x
            y = n // x
            if y > x:
                total += y
    return total

In your main loop, your search for i,j both in range(1,10000), and add both i, j (if different) to your set to ensure you don’t have duplicates.

You can cut down your search space by half by recognizing that if you find one pair i, j, you have also found the other pair j, i, and adding the other pair to the set is a no-op. You only need to search the cases where (say) i < j is true, which eliminates just over half of the pairs.

for i in range(1, 10000):
    for j in range(i+1, 10000):
        if sum_div(i) == j and sum_div(j) == i:
            amicable_nums.update([i, j])

How many times are you computing sum_div(i)??? Does it change in the inner loop? Maybe you could only compute it once per loop?

for i in range(1, 10000):
    sdi = sum_div(i)
    for j in range(i+1, 10000):
        if sdi == j and sum_div(j) == i:
            amicable_nums.update([i, j])

Do you really need the inner loop? sdi == j will only be true for exactly one value in the loop.

\$\endgroup\$
4
\$\begingroup\$

Another way to speed this up is to keep a record of all the sum of divisors seen so far (by using a caching decorator on the function, like functools.lru_cache) and realize that you only need a single loop if you make this a generator and use the fact that a number is amicable if \$d(d(a)) = a\$ and \$d(a) \neq a\$:

import math
from functools import lru_cache

@lru_cache(None)
def sum_div(n):
    # Taken from AJNeufeld's answer
    total = 1
    for x in range(2, int(math.sqrt(n) + 1)):
        if n % x == 0:
            total += x
            y = n // x
            if y > x:
                total += y
    return total

def amicable_numbers(limit):
    for a in range(limit):
        b = sum_div(a)
        if a != b and sum_div(b) == a:
            yield a

print(sum(amicable_numbers(10000)))

This runs in a bit more than 4 milli-seconds on my computer.


As for timing the runtime, I usually either prefer using ipythons magic command %timeit, or writing a small context manager:

from time import perf_counter

class Timer:
    def __init__(self, name=""):
        self.name = name

    def __enter__(self):
        self.start = perf_counter()

    def __exit__(self, *args):
        runtime = perf_counter() - self.start

        # get it in nice units
        units = ["s", "ms", "μs"]
        for unit in units:
            if runtime > 1:
                break
            runtime *= 1000

        if self.name:
            print(f"{self.name}: {runtime:.1f}{unit}")
        else:
            print(f"{runtime:.1f}{unit}")

Which you can use like this:

with Timer("amicable numbers"):
    print(sum(set(amicable_numbers(10000))))
# XXXXX  # I don't want to give away the correct answer
# amicable numbers: 4.1ms

Note that this will not be more precise than micro-seconds due to the time it takes to run the context manager.

\$\endgroup\$
  • \$\begingroup\$ I think you could be getting the wrong answer. If there exists an amicable pair a <= limit, b > limit, you will add a to the total, where as the original solution adds neither. You could fixed this with:if b < a and sum_div(b) == a: yield a; yield b. The summing the output of a generator is unnecessary, and will slow things down. For faster code, use total = 0 total += a + b and return total. \$\endgroup\$ – AJNeufeld Aug 30 '18 at 15:59
  • 1
    \$\begingroup\$ @AJNeufeld I specifically fixed it so that I do include a and don't include b in that case. a is still an amicable number, it's amicable pair is just greater than the limit. You are right on the running total, though. But since the timing is well under a second I stopped optimizing... \$\endgroup\$ – Graipher Aug 30 '18 at 16:04
  • \$\begingroup\$ The question is vague as to whether or not to include a if b is beyond the limit. Curiosity got the better of me. There is no amicable pair (a,b) where a <= limit and b > limit. \$\endgroup\$ – AJNeufeld Aug 30 '18 at 16:11
  • \$\begingroup\$ @AJNeufeld For limit = 10000 yes. I tested my code with limit = 250, which splits the example pair, to get the behavior I consider the best interpretation of the challenge rules. \$\endgroup\$ – Graipher Aug 30 '18 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.