En lugar del bucle doble sobre x y m y verificando repetidamente si son coprimos, iteramos solo sobre m (el mayor de los dos) y aplicamos la función totient de Euler o una versión personalizada de la misma para contar directamente el número de valores de x que son primos relativos a m . Esto nos da un método mucho más rápido (la velocidad aún debe cuantificarse con mayor precisión): por ejemplo, 43 ms para n = 100_000_000 en lugar de 30 s con el código del OP (700x de aceleración).

La necesidad de una versión personalizada surge cuando el valor máximo xmax que puede tomar x es menor que m (para satisfacer la desigualdad (m**2 + x**2)/2 <= n ). En ese caso, no se deben contar todos los coprimos de m , sino solo aquellos hasta ese límite.

 def distinct_factors(n): # a variant of the well-known factorization, but that # yields only distinct factors, rather than all of them # (including possible repeats) last = None i = 2 while i * i <= n: if n % i: i += 1 else: n //= i if i != last: yield i last = i if n > 1 and n != last: yield 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): # Euler's totient or "phi" function if upto is not None and upto < n: # custom version: all co-primes of n up to the `upto` bound cnt = upto p_list = list(distinct_factors(n)) for q in products_of(p_list, upto): cnt += upto // q if q > 0 else -(upto // -q) return cnt # standard formulation: all co-primes of n up to n-1 cnt = n for p in distinct_factors(n): cnt *= (1 - 1/p) return int(cnt)

phi(n) es el tociente de Euler o función ϕ(n) .

phi(n, upto=x) es una variante personalizada que cuenta solo los coprimos hasta un valor x dado. Para entenderlo, trabajemos con un ejemplo:

 >>> n = 3*3*3*5 # 135 >>> list(factors(n)) [3, 3, 3, 5] >>> list(distinct_factors(n)) [3, 5] # there are 72 integers between 1 and 135 that are co-primes of 135 >>> phi(n) 72 # ...but only 53 of them are no greater than 100: # 100 - (100//3 + 100//5 - 100//(3*5)) >>> phi(n, upto=100) 53

Al evaluar el número de coprimos de n bajo un valor x , debemos contar todos los números 1 .. x menos los que son múltiplos de cualquiera de los distintos factores de n . Sin embargo, cuando simplemente eliminamos x // p_i para todos los p_i , contamos dos veces los números que son múltiplos de dos factores, por lo que debemos "volver a agregarlos". Sin embargo, al hacerlo, contamos dos veces (sumamos demasiadas veces) los números que son múltiplos de tres factores, por lo que también debemos tenerlos en cuenta, etc. En el ejemplo n = 135 , eliminamos x // 3 y x // 5 , pero luego cuenta dos veces esos enteros que son factores de 3 y 5 (factores de 15), por lo que debemos volver a sumarlos. Para un conjunto más largo de factores, necesitamos:

• tomar x como cuenta inicial;
• restar el número de múltiplos de cada factor p ;
• "des-restar" (sumar) el número de múltiplos de cualquier producto de 2 factores;
• "un-un-subtract" (restar) el número de múltiplos de cualquier producto de 3 factores;
• etc

La respuesta inicial fue hacer esto iterando sobre todas las combinaciones de factores distintos, pero esto está sustancialmente optimizado en esta respuesta por el generador products_of(p_list, upto) , que da los productos de todos los subconjuntos de los factores distintos p_list dados cuyo producto no es mayor que upto . El signo indica cómo contabilizar cada producto: positiva o negativamente dependiendo de si el tamaño del subconjunto es par o impar, respectivamente.

Con phi(n) y phi(n, upto) en la mano, ahora podemos escribir lo siguiente:

 def pyth_m_counts(n): # yield tuples (m, count(x) where 0 < x < m and odd(x) # and odd(m) and coprime(x, m) and m**2 + x**2 <= 2*n). mmax = isqrt(2*n - 1) for m in range(3, mmax + 1, 2): # requirement: (m**2 + x**2) // 2 <= n # and both m and x are odd # (so (m**2 + x**2) // 2 == (m**2 + x**2) / 2) xmax = isqrt(2*n - m**2) cnt_m = phi(2*m, upto=xmax) if xmax < m else phi(2*m) // 2 if cnt_m > 0: yield m, cnt_m

¿Por qué la expresión phi(2*m) // 2 ? Dado que x (y m ) deben ser impares, según el OP, debemos eliminar todos los valores pares. Podemos hacer eso sin modificar phi() , pasando 2*m (que luego tiene 2 como factor, y por lo tanto "matará" todos los valores pares de x ) y luego dividiendo por 2 para obtener el número real de off co -primos a m . Una consideración similar (pero un poco más sutil) se hace con phi(2*m, upto=xmax) -lo dejaremos como ejercicio para el lector...

Ejemplo de ejecución:

 >>> n = 300 >>> list(pyth_m_counts(n)) [(3, 1), (5, 2), (7, 3), (9, 3), (11, 5), (13, 6), (15, 4), (17, 8), (19, 8), (21, 3), (23, 4)]

Eso significa que, en la función del OP, pythag_triples(300) habría devuelto 1 tupla con m==3 , 2 tuplas con m==5 , etc. De hecho, modifiquemos esa función para verificar esto:

 def mod_pythag_triples(n): 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: c = (m**2 + x**2) // 2 if c < n: yield x, m

Entonces:

 >>> n = 300 >>> list(pyth_m_counts(n)) == list(Counter(m for x, m in mod_pythag_triples(n)).items()) True

Lo mismo para cualquier valor positivo de n .

Ahora en la función de conteo real: solo necesitamos resumir los conteos para cada m:

 def pyth_triples_count(n): cnt = 0 mmax = isqrt(2*n - 1) for m in range(3, mmax + 1, 2): # requirement: (m**2 + x**2) // 2 <= n # and both m and x are odd (so (m**2 + x**2) // 2 == (m**2 + x**2) / 2) xmax = isqrt(2*n - m**2) cnt += phi(2*m, upto=xmax) if xmax < m else phi(2*m) // 2 return cnt

Ejecuciones de muestra:

 >>> pyth_triples_count(1_000_000) 159139 >>> pyth_triples_count(100_000_000) 15915492 >>> pyth_triples_count(1_000_000_000) 159154994 >>> big_n = 3_141_592_653_589_793 >>> pyth_triples_count(big_n) 500000000002845

Velocidad:

 %timeit pyth_triples_count(100_000_000) 42.8 ms ± 56.1 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) %timeit pyth_triples_count(1_000_000_000) 188 ms ± 571 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) %%time pyth_triples_count(big_n) CPU times: user 1h 42min 33s, sys: 480 ms, total: 1h 42min 33s Wall time: 1h 42min 33s

Nota: en la misma máquina, el código en la pregunta del OP tarda 30 segundos para n=100_000_000 ; esta versión es 700 veces más rápida para ese n .

Vea también mi otra respuesta para una solución más rápida.

Esta nueva respuesta reduce el tiempo total de big_n a 4 min 6 s .

Un perfil de mi respuesta inicial reveló estos hechos:

• Tiempo total: 1h 42min 33s
• Tiempo dedicado a factorizar números: casi el 100% del tiempo

Por el contrario, generar todos los números primos de 3 a sqrt(2*N - 1) lleva solo 38,5 s (usando el tamiz de Atkin).

Por lo tanto, decidí probar una versión en la que generamos todos los números m como productos conocidos de números primos. Es decir, el generador da el número mismo así como los distintos factores primos involucrados. No se necesita factorización.

El resultado sigue siendo 500_000_000_002_841 , con una diferencia de 4 como notó @Koder. Todavía no sé de dónde viene ese problema . Editar : después de la corrección del límite xmax ( isqrt(2*N - m**2) en lugar de isqrt(2*N - m**2 - 1) , ya que queremos incluir triángulos con hipotenusa igual a N ), ahora obtenemos el resultado correcto.

El código para el generador de números primos se incluye al final. Básicamente, utilicé el tamiz de Atkin , adaptado (sin dedicarle mucho tiempo) a Python. Estoy bastante seguro de que podría acelerarse (por ejemplo, usando numpy y quizás incluso numba ).

Para generar números enteros a partir de números primos (lo cual sabemos que podemos hacer gracias alteorema fundamental de la aritmética ), solo necesitamos iterar a través de todos los productos posibles prod(p_i**k_i) donde p_i es el i^th número primo y k_i es cualquier entero no negativo.

La formulación más fácil es recursiva:

 def gen_ints_from_primes(p_list, upto): if p_list and upto >= p_list[0]: p, *p_list = p_list pk = 1 p_tup = tuple() while pk <= upto: for q, p_distinct in gen_ints_from_primes(p_list, upto=upto // pk): yield pk * q, p_tup + p_distinct pk *= p p_tup = (p, ) else: yield 1, tuple()

Desafortunadamente, rápidamente nos encontramos con restricciones de memoria (y límite de recurrencia). Entonces, aquí hay una versión no recursiva que no usa memoria adicional aparte de la lista de números primos. Esencialmente, el valor actual de q (el entero en proceso de generarse) y un índice en la lista son toda la información que necesitamos para generar el siguiente entero. Por supuesto, los valores vienen sin ordenar, pero eso no importa, siempre y cuando estén todos cubiertos.

 def rem_p(q, p, p_distinct): q0 = q while q % p == 0: q //= p if q != q0: if p_distinct[-1] != p: raise ValueError(f'rem({q}, {p}, ...{p_distinct[-4:]}): p expected at end of p_distinct if q % p == 0') p_distinct = p_distinct[:-1] return q, p_distinct def add_p(q, p, p_distinct): if len(p_distinct) == 0 or p_distinct[-1] != p: p_distinct += (p, ) q *= p return q, p_distinct def gen_prod_primes(p, upto=None): if upto is None: upto = p[-1] if upto >= p[-1]: p = p + [upto + 1] # sentinel q = 1 i = 0 p_distinct = tuple() while True: while q * p[i] <= upto: i += 1 while q * p[i] > upto: yield q, p_distinct if i <= 0: return q, p_distinct = rem_p(q, p[i], p_distinct) i -= 1 q, p_distinct = add_p(q, p[i], p_distinct)

Ejemplo-

 >>> p_list = list(primes(20)) >>> p_list [2, 3, 5, 7, 11, 13, 17, 19] >>> sorted(gen_prod_primes(p_list, 20)) [(1, ()), (2, (2,)), (3, (3,)), (4, (2,)), (5, (5,)), (6, (2, 3)), (7, (7,)), (8, (2,)), (9, (3,)), (10, (2, 5)), (11, (11,)), (12, (2, 3)), (13, (13,)), (14, (2, 7)), (15, (3, 5)), (16, (2,)), (17, (17,)), (18, (2, 3)), (19, (19,)), (20, (2, 5))]

Como puede ver, no necesitamos factorizar ningún número, ya que convenientemente vienen con los distintos números primos involucrados.

Para obtener solo números impares, simplemente elimine 2 de la lista de números primos:

 >>> sorted(gen_prod_primes(p_list[1:]), 20) [(1, ()), (3, (3,)), (5, (5,)), (7, (7,)), (9, (3,)), (11, (11,)), (13, (13,)), (15, (3, 5)), (17, (17,)), (19, (19,))]

Para explotar esta presentación de números y factores, necesitamos modificar un poco la función dada en la respuesta original:

 def phi(n, upto=None, p_list=None): # Euler's totient or "phi" function if upto is None or upto > n: upto = n if p_list is None: p_list = list(distinct_factors(n)) if upto < n: # custom version: all co-primes of n up to the `upto` bound cnt = upto for q in products_of(p_list, upto): cnt += upto // q if q > 0 else -(upto // -q) return cnt # standard formulation: all co-primes of n up to n-1 cnt = n for p in p_list: cnt = cnt * (p - 1) // p return cnt

Con todo esto, ya podemos reescribir nuestras funciones de conteo:

 def pt_count_m(N): # yield tuples (m, count(x) where 0 < x < m and odd(x) # and odd(m) and coprime(x, m) and m**2 + x**2 <= 2*N)) # in this version, m is generated from primes, and the values # are iterated through unordered. mmax = isqrt(2*N - 1) p_list = list(primes(mmax))[1:] # skip 2 for m, p_distinct in gen_prod_primes(p_list, upto=mmax): if m < 3: continue # requirement: (m**2 + x**2) // 2 <= N # note, both m and x are odd (so (m**2 + x**2) // 2 == (m**2 + x**2) / 2) xmax = isqrt(2*N - m*m) cnt_m = phi(m+1, upto=xmax, p_list=(2,) + tuple(p_distinct)) if cnt_m > 0: yield m, cnt_m def pt_count(N, progress=False): mmax = isqrt(2*N - 1) it = pt_count_m(N) if progress: it = tqdm(it, total=(mmax - 3 + 1) // 2) return sum(cnt_m for m, cnt_m in it)

Y ahora:

 %timeit pt_count(100_000_000) 31.1 ms ± 38.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) %timeit pt_count(1_000_000_000) 104 ms ± 299 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) # the speedup is still very moderate at that stage # however: %%time big_n = 3_141_592_653_589_793 N = big_n res = pt_count(N) CPU times: user 4min 5s, sys: 662 ms, total: 4min 6s Wall time: 4min 6s >>> res 500000000002845 # sin embargo: %timeit pt_count(100_000_000) 31.1 ms ± 38.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) %timeit pt_count(1_000_000_000) 104 ms ± 299 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) # the speedup is still very moderate at that stage # however: %%time big_n = 3_141_592_653_589_793 N = big_n res = pt_count(N) CPU times: user 4min 5s, sys: 662 ms, total: 4min 6s Wall time: 4min 6s >>> res 500000000002845

Apéndice: tamiz de Atkin

Como prometí, aquí está mi versión del tamiz de Atkin. Definitivamente se puede acelerar.

 def primes(limit): # Generates prime numbers between 2 and n # Atkin's sieve -- see http://en.wikipedia.org/wiki/Prime_number sqrtLimit = isqrt(limit) + 1 # initialize the sieve is_prime = [False, False, True, True, False] + [False for _ in range(5, limit + 1)] # put in candidate primes: # integers which have an odd number of # representations by certain quadratic forms for x in range(1, sqrtLimit): x2 = x * x for y in range(1, sqrtLimit): y2 = y*y n = 4 * x2 + y2 if n <= limit and (n % 12 == 1 or n % 12 == 5): is_prime[n] ^= True n = 3 * x2 + y2 if n <= limit and (n % 12 == 7): is_prime[n] ^= True n = 3*x2-y2 if n <= limit and x > y and n % 12 == 11: is_prime[n] ^= True # eliminate composites by sieving for n in range(5, sqrtLimit): if is_prime[n]: sqN = n**2 # n is prime, omit multiples of its square; this is sufficient because # composites which managed to get on the list cannot be square-free for i in range(1, int(limit/sqN) + 1): k = i * sqN # k ∈ {n², 2n², 3n², ..., limit} is_prime[k] = False for i, truth in enumerate(is_prime): if truth: yield i

Gracias a Pierre encontré una solución mucho más rápida.

Aquí está mi nuevo código combinado con el de Pierre para cualquiera que lo desee.

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

Como se mencionó en la respuesta anterior, la mayor parte del tiempo se dedicó a la factorización, así que tomé su código y agregué un tamiz de todos los números hasta x-max, siendo cada número el índice que produce su factor primo más bajo. Encuentra la respuesta en 4 minutos y 44 segundos. (284.6727148 segundos). Gracias por la ayuda Pedro.