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)

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$.

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))