def pythag_triples(n):
i = 0
start = time.time()
for x in range(1, int(sqrt(n) + sqrt(n)) + 1, 2):
for m in range(x+2,int(sqrt(n) + sqrt(n)) + 1, 2):
if gcd(x, m) == 1:
# q = x*m
# l = (m**2 - x**2)/2
c = (m**2 + x**2)/2
# trips.append((q,l,c))
if c < n:
i += 1
end = time.time()
return i, end-start
print(pythag_triples(3141592653589793))
I'm trying to calculate primitive pythagorean triples using the idea that all triples are generated from using m, n that are both odd and coprime. I already know that the function works up to 1000000 but when doing it to the larger number its taken longer than 24 hours. Any ideas on how to speed this up/ not brute force it. I am trying to count the triples.
This new answer brings the total time for
big_n
down to 4min 6s.An profiling of my initial answer revealed these facts:
In contrast, generating all primes from
3
tosqrt(2*N - 1)
takes only 38.5s (using Atkin's sieve).I therefore decided to try a version where we generate all numbers
m
as known products of prime numbers. That is, the generator yields the number itself as well as the distinct prime factors involved. No factorization needed.The result is still. Edit: after correction of the500_000_000_002_841
, off by 4 as @Koder noticed. I do not know yet where that problem comes fromxmax
bound (isqrt(2*N - m**2)
instead ofisqrt(2*N - m**2 - 1)
, since we do want to include triangles with hypothenuse equal toN
), we now get the correct result.The code for the primes generator is included at the end. Basically, I used Atkin's sieve, adapted (without spending much time on it) to Python. I am quite sure it could be sped up (e.g. using
numpy
and perhaps evennumba
).To generate integers from primes (which we know we can do thanks to the Fundamental theorem of arithmetic), we just need to iterate through all the possible products
prod(p_i**k_i)
wherep_i
is thei^th
prime number andk_i
is any non-negative integer.The easiest formulation is a recursive one:
Unfortunately, we quickly run into memory constraints (and recursion limit). So here is a non-recursive version which uses no extra memory aside from the list of primes themselves. Essentially, the current value of
q
(the integer in process of being generated) and an index in the list are all the information we need to generate the next integer. Of course, the values come unsorted, but that doesn't matter, as long as they are all covered.Example-
As you can see, we don't need to factorize any number, as they conveniently come along with the distinct primes involved.
To get only odd numbers, simply remove
2
from the list of primes:In order to exploit this number-and-factors presentation, we need to amend a bit the function given in the original answer:
With all this, we can now rewrite our counting functions:
And now:
Addendum: Atkin's sieve
As promised, here is my version of Atkin's sieve. It can definitely be sped up.