Module de factorisation rapide des nombres premiers

Question

Je suis à la recherche d'un mise en œuvre ou algorithme clair pour obtenir les facteurs premiers de N en python, pseudocode ou tout autre langage lisible. Il y a quelques exigences/contraintes :

• N est compris entre 1 et ~20 chiffres
• Pas de table de recherche pré-calculée, mais la mémorisation est possible.
• N'ont pas besoin d'être prouvées mathématiquement (par exemple, elles pourraient s'appuyer sur la conjecture de Goldbach si nécessaire).
• N'a pas besoin d'être précis, peut être probabiliste/déterministe si nécessaire.

J'ai besoin d'un algorithme de factorisation rapide des nombres premiers, non seulement pour lui-même, mais aussi pour être utilisé dans de nombreux autres algorithmes, comme le calcul de l'indice d'Euler. phi(n) .

J'ai essayé d'autres algorithmes tirés de Wikipedia et d'autres sources, mais soit je n'ai pas pu les comprendre (ECM), soit je n'ai pas pu créer une implémentation fonctionnelle à partir de l'algorithme (Pollard-Brent).

Je suis vraiment intéressé par l'algorithme de Pollard-Brent, donc toute information ou implémentation supplémentaire à ce sujet serait la bienvenue.

Merci !

EDITAR

Après avoir bricolé un peu, j'ai créé un module d'amorçage/factorisation assez rapide. Il combine un algorithme de division d'essai optimisé, l'algorithme de Pollard-Brent, un test de primalité de Miller-rabin et le primieve le plus rapide que j'ai trouvé sur internet. gcd est une implémentation du PGCD d'Euclide régulier (le PGCD d'Euclide binaire est beaucoup plus lent que le normal).

Bounty

Oh joie, une abondance peut être acquise ! Mais comment la gagner ?

• Trouver une optimisation ou un bug dans mon module.
• Fournir des algorithmes/implémentations alternatifs/meilleurs.

La réponse la plus complète/constructive reçoit la prime.

Et enfin le module lui-même :

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

Il n'est pas nécessaire de calculer smallprimes en utilisant primesbelow utiliser smallprimeset pour ça.

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

Divisez votre primefactors en deux fonctions pour le traitement smallprimes et autres pour pollard_brent cela peut permettre d'économiser quelques itérations car toutes les puissances des petits nombres seront divisées par 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

En tenant compte des résultats vérifiés de Pomerance, Selfridge et Wagstaff et Jaeschke, vous pouvez réduire les répétitions en isprime qui utilise le test de primauté de Miller-Rabin. À partir de Wiki .

• si n < 1,373,653, il suffit de tester a = 2 et 3 ;
• si n < 9 080 191, il suffit de tester a = 31 et 73 ;
• si n < 4,759,123,141, il suffit de tester a = 2, 7, et 61 ;
• si n < 2,152,302,898,747, il suffit de tester a = 2, 3, 5, 7 et 11 ;
• si n < 3,474,749,660,383, il suffit de tester a = 2, 3, 5, 7, 11 et 13 ;
• si n < 341,550,071,728,321, il suffit de tester a = 2, 3, 5, 7, 11, 13 et 17.

Edit 1 : Correction de l'appel de retour de if-else pour ajouter les bigfactors aux facteurs dans primefactors .

Answer 2

Pollard-Brent en implémenté en Python :

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

Answer 3

Même sur l'actuelle, il y a plusieurs endroits à remarquer.

1. Ne le faites pas. checker*checker à chaque boucle, utilisez s=ceil(sqrt(num)) y checher < s
2. checher doit être plus 2 à chaque fois, ignorer tous les nombres pairs sauf 2
3. Utilice divmod au lieu de % y //

Answer 4

Vous devriez probablement faire un peu de détection de prime que vous pourriez regarder ici, Algorithme rapide pour trouver des nombres premiers ?

Vous devriez lire ce blog en entier, il y a quelques algorithmes qu'il énumère pour tester la primalité.

Answer 5

Il existe une bibliothèque python contenant une collection de tests de primalité (y compris des tests incorrects sur ce qu'il ne faut pas faire). Elle s'appelle pyprimes . J'ai pensé que cela valait la peine de le mentionner pour la postérité. Je ne pense pas qu'il comprenne les algorithmes que vous avez mentionnés.