2
\$\begingroup\$

Problem:

Consider equations of the form: \$a^2 + b^2 = N; 0 \leq a \leq b; a, b, N \in \mathbb{N}\$.

For \$N=65\$ there are two solutions:

\$a=1, b=8\$ and \$a=4, b=7\$.

We call \$S(N)\$ the sum of the values of \$a\$ of all solutions of \$a^2 + b^2 = N\$.

Thus \$S(65) = 1 + 4 = 5\$.

Find \$\sum S(N)\$, for all squarefree \$N\$ only divisible by primes of the form \$4k+1\$ with \$4k+1 < 150\$.

My solution is painfully slow:

import math
import itertools
import time

def candidate_range(n):
    cur = 5
    incr = 2
    while cur < n+1:
        yield cur
        cur += incr
        incr ^= 6 # or incr = 6-incr, or however
def sieve(end):
    prime_list = [2, 3]
    sieve_list = [True] * (end+1)
    for each_number in candidate_range(end):
        if sieve_list[each_number]:
            prime_list.append(each_number)
            for multiple in range(each_number*each_number, end+1, each_number):
                sieve_list[multiple] = False
    return prime_list

primes = sieve(150)
goodprimes = []
for prime in primes:
    if prime%(4)==1:
        goodprimes.append(prime)
sum=[]
start_time = time.time()
#get a number that works
print("-------Part 1------")
mi=0
for L in range(1, len(goodprimes)+1):
    sumf=0
    for subset in itertools.combinations(goodprimes, L):
        max=2**L/2
        n=1
        for x in subset:
            n*=x

        for b in range(math.floor(math.sqrt(n/2)), math.floor(math.sqrt(n)+1)):
            a=math.sqrt(n-b*b)
            if a.is_integer() and b>=a:
                sum.append(a)
                mi+=1
                if mi==max:
                    mi=0
                    break
    for num in sum:
        sumf+=num
    print(L,sumf, "--- %s seconds ---" % (time.time() - start_time))

                    #q+=1

print("--- %s seconds ---" % (time.time() - start_time))
sumf=0
for num in sum:
    sumf+=num
print(sumf)
\$\endgroup\$
  • \$\begingroup\$ I find beautiful that N has 65535 values. \$\endgroup\$ – Astrinus Mar 14 at 11:17
  • \$\begingroup\$ .... and that you will need a BigInteger library because the topmost one will require 92 bit to be represented. \$\endgroup\$ – Astrinus Mar 14 at 11:34
  • 1
    \$\begingroup\$ @Astrinus Python has unlimited size integers built in. \$\endgroup\$ – mkrieger1 Mar 14 at 19:05
  • \$\begingroup\$ @mkrieger1 en.wikipedia.org/wiki/… will solve the problem better ;-) \$\endgroup\$ – Astrinus Mar 15 at 10:10
4
\$\begingroup\$

The code

The formatting has a number of PEP8 violations.


It's not obvious from the name what candidate_range does. It seems to be a wheel for the sieve. Normally that would be inlined in the sieve; even if you prefer not to do that, you could place the function inside sieve to make its scope clear.

I don't find sieve_list a very helpful name. In general for sieving I prefer is_composite, inverting the booleans from the way you've done it. Similarly for each_number: it reads well on the first line which uses it, but very oddly on the others.


goodprimes = []
for prime in primes:
    if prime%(4)==1:
        goodprimes.append(prime)

It's more Pythonic to use comprehensions:

goodprimes = [p for p in primes if p % 4 == 1]

#get a number that works

What does this mean? It looks more like noise than a useful comment to me.


for L in range(1, len(goodprimes)+1):
    sumf=0
    for subset in itertools.combinations(goodprimes, L):

I don't know why itertools doesn't have a function to give all subsets, but it seems like the kind of thing which is worth pulling out as a separate function, both for reuse and for readability.


        max=2**L/2

What does this do?


        n=1
        for x in subset:
            n*=x

Consider as an alternative

from functools import reduce
import operator

    n = reduce(operator.mul, subset, 1)

        for b in range(math.floor(math.sqrt(n/2)), math.floor(math.sqrt(n)+1)):
            a=math.sqrt(n-b*b)
            if a.is_integer() and b>=a:

Why floors rather than ceils?

Are you certain that math.sqrt on an integer is never out by 1ULP?

Why is b>=a necessary? (Obviously b==a is impossible, and isn't the point of the range chosen to force b > a?)


                sum.append(a)

Is this for debugging? I can't see why you wouldn't just add a to a total.

NB sum is aliasing the builtin function for adding up the values in a list.


                    #q+=1

??? I can't see any other mention of q.

The algorithm

There are a few Project Euler problems which fall to brute force, but in general you need to find the right mathematical insight. Given the way this question is structured, you probably need to figure out how to find \$S(n)\$ given the prime factorisation of \$n\$.

\$\endgroup\$
0
\$\begingroup\$

sieve_list

Why don't you use a generator here? It can be clearer as generator, and the candidate_range proves you know how those work.

def sieve_OP_gen(end):
    yield 2
    yield 3
    sieve_list = [True] * (end+1)
    for each_number in candidate_range(end):
        if sieve_list[each_number]:
            yield each_number
            for multiple in range(each_number*each_number, end+1, each_number):
                sieve_list[multiple] = False

list slice assignment

instead of:

for multiple in range(each_number*each_number, end+1, each_number):
    sieve_list[multiple] = False

you can do:

sieve_list[each_number::each_number] = [False] * (end // each_number)

this doesn't provide any speedup, but is more clear to me

candidate_range

I don't like incr ^= 6. This can be done a lot clearer with itertools.cycle

def candidate_range_maarten():
    cur = 5
    increments = itertools.cycle((2, 4))
    while True:
        yield cur
        cur += next(increments)

But all in all I think this is a lot of effort to reduce the number of checks in generating the primes sieve by 1/3rd. In fact, it slows down the sieve generation

def sieve2(end):
    yield 2
    sieve_list = [True] * (end + 1)
    for each_number in range(3, end + 1, 2):
        if not sieve_list[each_number]:
            continue
        yield each_number
        sieve_list[each_number::each_number] = [False] * (end // each_number)
sieve_OP(150) == list(sieve2(150))
True

timings:

%timeit sieve_OP(150)
24.5 µs ± 1.63 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit list(sieve2(150))
16.3 µs ± 124 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

filtering primes

to filter the primes, you can use the builtin filter

primes = list(sieve2(150))
goodprimes = list(filter(lambda x: x % 4 == 1, primes))

or a list comprehension: good_primes = [i for i in primes if i % 4 == 1]

functions

The rest of the code would be more clear if you split it in different functions. One to find the different candidates for the products, and another function to generate the pythagorean a

product:

The product of an iterable can be calculated like this:

from functools import reduce
from operator import mul

def prod(iterable):
    return reduce(mul, iterable, 1)

powerset

As @PeterTaylor tnoted, there might be a itertools function to do this. There is,'t but there is an itertools recipe powerset:

def powerset(iterable):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return itertools.chain.from_iterable(
        itertools.combinations(s, r) for r in range(len(s) + 1)
    )

So generating the candidates is as easy as

def candidates(good_primes):
    for subset in powerset(good_primes):
        yield prod(subset)

pythagorean_a

Instead of that nested for-loop which is not very clear what happens, I would split this to another function:

def pythagorean_a(n):
    for a in itertools.count(1):
        try:
            b = sqrt(n - (a ** 2))
        except ValueError:
            return
        if b < a:
            return
        if b.is_integer():
            yield a
list(pythagorean_a(65))
[1, 4]

bringing it together

\$S(N)\$ then becomes: sum(pythagorean_a(i))

and \$\sum S(N)\$: sum(sum(pythagorean_a(i)) for i in candidates(good_primes))

\$\endgroup\$

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.