prime number geeksforgeeks geeks for python algorithm prime-factoring

python - number - Módulo de factorización prima rápida



prime factors geeks for geeks (8)

Acabo de encontrar un error en este código al factorizar el número 2**1427 * 31 .

File "buckets.py", line 48, in prettyprime factors = primefactors.primefactors(n, sort=True) File "/private/tmp/primefactors.py", line 83, in primefactors limit = int(n ** .5) + 1 OverflowError: long int too large to convert to float

Este fragmento de código:

limit = int(n ** .5) + 1 for checker in smallprimes: if checker > limit: break while n % checker == 0: factors.append(checker) n //= checker limit = int(n ** .5) + 1 if checker > limit: break

debe ser reescrito en

for checker in smallprimes: while n % checker == 0: factors.append(checker) n //= checker if checker > n: break

lo que probablemente sea más rápido en las entradas realistas de todos modos. La raíz cuadrada es lenta, básicamente el equivalente de muchas multiplicaciones, las smallprimes solo tienen unas pocas docenas de miembros, y de esta manera eliminamos el cálculo de n ** .5 del círculo interno ajustado, lo que es útil cuando se factorizan números como 2**1427 . Simplemente no hay ninguna razón para calcular sqrt(2**1427) , sqrt(2**1426) , sqrt(2**1425) , etc., cuando lo único que nos importa es "hace el [cuadrado del] verificador excede n ".

El código reescrito no arroja excepciones cuando se presentan con números grandes, y es aproximadamente el doble de rápido de acuerdo con el timeit (en las entradas de muestra 2 y 2**718 * 31 ).

También isprime(2) cuenta que isprime(2) devuelve el resultado incorrecto, pero está bien siempre que no dependamos de él. En mi humilde opinión, debe volver a escribir la introducción de esa función como

if n <= 3: return n >= 2 ...

Estoy buscando una implementación o un algoritmo claro para obtener la factorización prima de N en Python, pseudocódigo o cualquier otra cosa que sea bien legible. Hay algunas demandas / hechos:

  • N tiene entre 1 y ~ 20 dígitos
  • Sin embargo, ninguna tabla de búsqueda precalculada está bien.
  • No es necesario demostrar matemáticamente (por ejemplo, podría confiar en la conjetura de Goldbach si fuera necesario)
  • No necesita ser preciso, se le permite ser probabilístico / determinista si es necesario

Necesito un algoritmo de factorización de prima rápida, no solo para sí mismo, sino para el uso en muchos otros algoritmos, como el cálculo de Euler phi (n) .

He intentado con otros algoritmos de Wikipedia y similares, pero o bien no pude entenderlos (ECM) o no pude crear una implementación en funcionamiento a partir del algoritmo (Pollard-Brent).

Estoy realmente interesado en el algoritmo de Pollard-Brent, por lo que cualquier información / implementación en él sería realmente agradable.

¡Gracias!

EDITAR

Después de jugar un poco, he creado un módulo de factorización / factorización bastante rápido. Combina un algoritmo optimizado de división de prueba, el algoritmo Pollard-Brent, una prueba de primalidad Miller-Rabin y el primesieve más rápido que encontré en Internet. gcd es una implementación regular de GCD de Euclid (el GCD binario de Euclid es mucho más lento que el ordinario).

Generosidad

¡Oh alegría, se puede adquirir una recompensa! Pero, ¿cómo puedo ganarlo?

  • Encuentre una optimización o error en mi módulo.
  • Proporcione algoritmos / implementaciones alternativos / mejores.

La respuesta que es la más completa / constructiva obtiene la recompensa.

Y finalmente el módulo en sí:

import random def primesbelow(N): # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188 #""" Input N>=6, Returns a list of primes, 2 <= p < N """ correction = N % 6 > 1 N = {0:N, 1:N-1, 2:N+4, 3:N+3, 4:N+2, 5:N+1}[N%6] sieve = [True] * (N // 3) sieve[0] = False for i in range(int(N ** .5) // 3 + 1): if sieve[i]: k = (3 * i + 1) | 1 sieve[k*k // 3::2*k] = [False] * ((N//6 - (k*k)//6 - 1)//k + 1) sieve[(k*k + 4*k - 2*k*(i%2)) // 3::2*k] = [False] * ((N // 6 - (k*k + 4*k - 2*k*(i%2))//6 - 1) // k + 1) return [2, 3] + [(3 * i + 1) | 1 for i in range(1, N//3 - correction) if sieve[i]] smallprimeset = set(primesbelow(100000)) _smallprimeset = 100000 def isprime(n, precision=7): # http://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time if n < 1: raise ValueError("Out of bounds, first argument must be > 0") elif n <= 3: return n >= 2 elif n % 2 == 0: return False elif n < _smallprimeset: return n in smallprimeset d = n - 1 s = 0 while d % 2 == 0: d //= 2 s += 1 for repeat in range(precision): a = random.randrange(2, n - 2) x = pow(a, d, n) if x == 1 or x == n - 1: continue for r in range(s - 1): x = pow(x, 2, n) if x == 1: return False if x == n - 1: break else: return False return True # https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/ def pollard_brent(n): if n % 2 == 0: return 2 if n % 3 == 0: return 3 y, c, m = random.randint(1, n-1), random.randint(1, n-1), random.randint(1, n-1) g, r, q = 1, 1, 1 while g == 1: x = y for i in range(r): y = (pow(y, 2, n) + c) % n k = 0 while k < r and g==1: ys = y for i in range(min(m, r-k)): y = (pow(y, 2, n) + c) % n q = q * abs(x-y) % n g = gcd(q, n) k += m r *= 2 if g == n: while True: ys = (pow(ys, 2, n) + c) % n g = gcd(abs(x - ys), n) if g > 1: break return g smallprimes = primesbelow(1000) # might seem low, but 1000*1000 = 1000000, so this will fully factor every composite < 1000000 def primefactors(n, sort=False): factors = [] for checker in smallprimes: while n % checker == 0: factors.append(checker) n //= checker if checker > n: break if n < 2: return factors while n > 1: if isprime(n): factors.append(n) break factor = pollard_brent(n) # trial division did not fully factor, switch to pollard-brent factors.extend(primefactors(factor)) # recurse to factor the not necessarily prime factor returned by pollard-brent n //= factor if sort: factors.sort() return factors def factorization(n): factors = {} for p1 in primefactors(n): try: factors[p1] += 1 except KeyError: factors[p1] = 1 return factors totients = {} def totient(n): if n == 0: return 1 try: return totients[n] except KeyError: pass tot = 1 for p, exp in factorization(n).items(): tot *= (p - 1) * p ** (exp - 1) totients[n] = tot return tot def gcd(a, b): if a == b: return a while b > 0: a, b = b, a % b return a def lcm(a, b): return abs((a // gcd(a, b)) * b)


Hay una biblioteca de Python con una colección de pruebas de primalidad (incluidas las incorrectas para lo que no se debe hacer). Se llama pyprimes . Supuse que vale la pena mencionarlo para el propósito de la posteridad. No creo que incluya los algoritmos que mencionaste.


Incluso en el actual, hay varios lugares para notar.

  1. No haga checker*checker cada ciclo, use s=ceil(sqrt(num)) y checher < s
  2. checher debería tener más 2 cada vez, ignorar todos los números pares excepto 2
  3. Use divmod lugar de % y //

No hay necesidad de calcular smallprimes usando primesbelow , use smallprimeset para eso.

smallprimes = (2,) + tuple(n for n in xrange(3,1000,2) if n in smallprimeset)

Divida sus smallprimes en dos funciones para el manejo de las smallprimes y otras para las pollard_brent . Esto puede ahorrar un par de iteraciones, ya que todas las potencias de las pequeñas ganancias se dividirán desde n.

def primefactors(n, sort=False): factors = [] limit = int(n ** .5) + 1 for checker in smallprimes: print smallprimes[-1] if checker > limit: break while n % checker == 0: factors.append(checker) n //= checker if n < 2: return factors else : factors.extend(bigfactors(n,sort)) return factors def bigfactors(n, sort = False): factors = [] while n > 1: if isprime(n): factors.append(n) break factor = pollard_brent(n) factors.extend(bigfactors(factor,sort)) # recurse to factor the not necessarily prime factor returned by pollard-brent n //= factor if sort: factors.sort() return factors

Al considerar los resultados verificados de Pomerance, Selfridge y Wagstaff y Jaeschke, puede reducir las repeticiones en isprime que usa la prueba de primalidad Miller-Rabin. De Wiki .

  • si n <1.373.653, es suficiente para probar a = 2 y 3;
  • si n <9,080,191, es suficiente para probar a = 31 y 73;
  • si n <4,759,123,141, es suficiente probar a = 2, 7 y 61;
  • si n <2.152.302.898.747, es suficiente probar a = 2, 3, 5, 7 y 11;
  • si n <3,474,749,660,383, es suficiente probar a = 2, 3, 5, 7, 11 y 13;
  • si n <341,550,071,728,321, es suficiente probar a = 2, 3, 5, 7, 11, 13 y 17.

Edición 1 : llamada de devolución corregida de if-else para anexar bigfactores a factores en factores primos.




Puedes factorizar hasta un límite y luego usar brent para obtener factores más altos

from fractions import gcd from random import randint def brent(N): if N%2==0: return 2 y,c,m = randint(1, N-1),randint(1, N-1),randint(1, N-1) g,r,q = 1,1,1 while g==1: x = y for i in range(r): y = ((y*y)%N+c)%N k = 0 while (k<r and g==1): ys = y for i in range(min(m,r-k)): y = ((y*y)%N+c)%N q = q*(abs(x-y))%N g = gcd(q,N) k = k + m r = r*2 if g==N: while True: ys = ((ys*ys)%N+c)%N g = gcd(abs(x-ys),N) if g>1: break return g def factorize(n1): if n1==0: return [] if n1==1: return [1] n=n1 b=[] p=0 mx=1000000 while n % 2 ==0 : b.append(2);n//=2 while n % 3 ==0 : b.append(3);n//=3 i=5 inc=2 while i <=mx: while n % i ==0 : b.append(i); n//=i i+=inc inc=6-inc while n>mx: p1=n while p1!=p: p=p1 p1=brent(p) b.append(p1);n//=p1 if n!=1:b.append(n) return sorted(b) from functools import reduce #n= 2**1427 * 31 # n= 67898771 * 492574361 * 10000223 *305175781* 722222227*880949 *908909 li=factorize(n) print (li) print (n - reduce(lambda x,y :x*y ,li))


Si no quiere reinventar la rueda, use la biblioteca sympy

pip install sympy

Usa la función sympy.ntheory.factorint

>>> from sympy.ntheory import factorint >>> factorint(10**20+1) {73: 1, 5964848081: 1, 1676321: 1, 137: 1}

Puedes factorizar algunos números muy grandes:

>>> factorint(10**100+1) {401: 1, 5964848081: 1, 1676321: 1, 1601: 1, 1201: 1, 137: 1, 73: 1, 129694419029057750551385771184564274499075700947656757821537291527196801: 1}