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

prime factors geeks for geeks

Question

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)

Answer 1

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

Answer 2

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.

Answer 3

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

Answer 4

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.

Answer 5

Pollard-Brent en implementación en Python:

https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/

Answer 6

Probablemente deberías hacer una detección primaria que podrías mirar aquí. Algoritmo rápido para encontrar números primos.

Aunque deberías leer todo el blog, hay algunos algoritmos que enumera para probar la primalidad.

Answer 7

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

Answer 8

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}