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.
Thanks to Pierre I found a much faster solution.
Here is my new code meshed with Pierre's for anyone wanting it.
def sieve_factors(n):
s = [0] * (n+1)
s[1] = 1
for i in range(2, n+1, 2):
s[i] = 2
for i in range(3, n+1, 2):
if s[i] == 0:
s[i] = i
for j in range(i, n + 1, i):
if s[j] == 0:
s[j] = i
return s
Q = sieve_factors(2*(isqrt(2 * 3141592653589793) + 1))
def findfactors(n):
global Q
yield Q[n]
last = Q[n]
while n > 1:
if Q[n] != last and Q[n] != 1:
last = Q[n]
yield Q[n]
n //= Q[n]
def products_of(p_list, upto):
for i, p in enumerate(p_list):
if p > upto:
break
yield -p
for q in products_of(p_list[i+1:], upto=upto // p):
yield -p * q
def phi(n, upto=None):
if upto is not None and upto < n:
cnt = upto
p_list = list(findfactors(n))
for q in products_of(p_list, upto):
cnt += upto // q if q > 0 else -(upto // -q)
return cnt
cnt = n
for p in findfactors(n):
cnt *= (1 - 1/p)
return int(cnt)
def countprimtrips(n):
cnt = 0
for m in range(3, int(sqrt(2*n)) + 1, 2):
xmax = int(sqrt(2*n - m**2))
cnt += phi(2*m, upto=xmax) if xmax < m else phi(2*m) // 2
return cnt
print(countprimtrips(3141592653589793))
As mentioned in the answer above most of the time was spent factoring so I took his code and added a sieve of all the numbers up to the x-max with each number being the index that yields their lowest prime factor. It finds the answer in 4 minutes and 44 seconds. (284.6727148 seconds). Thanks for the help Pierre.