Ich suche nach einer Implementierung oder einem klaren Algorithmus , um die Primfaktorisierung von N entweder in Python, Pseudocode oder irgendetwas anderem zu erhalten, das gut lesbar ist. Es gibt einige Forderungen / Fakten:

• N liegt zwischen 1 und ~ 20 Stellen
• Keine vorberechnete Nachschlagetabelle, Memoisierung ist jedoch in Ordnung.
• Muss nicht mathematisch bewiesen sein (könnte sich bei Bedarf auf die Goldbach-Vermutung stützen)
• Muss nicht präzise sein, darf bei Bedarf probabilistisch / deterministisch sein

Ich brauche einen schnellen Primfaktorisierungsalgorithmus, nicht nur für sich selbst, sondern auch für die Verwendung in vielen anderen Algorithmen wie der Berechnung des Euler- Phi (n) .

Ich habe andere Algorithmen von Wikipedia und dergleichen ausprobiert, aber entweder konnte ich sie nicht verstehen (ECM) oder ich konnte keine funktionierende Implementierung aus dem Algorithmus (Pollard-Brent) erstellen.

Ich interessiere mich sehr für den Pollard-Brent-Algorithmus, daher wären weitere Informationen / Implementierungen dazu sehr hilfreich.

Vielen Dank!

BEARBEITEN

Nachdem ich ein wenig herumgespielt habe, habe ich ein ziemlich schnelles Prim / Faktorisierungsmodul erstellt. Es kombiniert einen optimierten Trial-Division-Algorithmus, den Pollard-Brent-Algorithmus, einen Miller-Rabin-Primalitätstest und das schnellste Primesieve, das ich im Internet gefunden habe. gcd ist eine reguläre GCD-Implementierung von Euclid (binäre GCD von Euclid ist viel langsamer als die reguläre).

Kopfgeld

Oh Freude, ein Kopfgeld kann erworben werden! Aber wie kann ich es gewinnen?

• Finden Sie eine Optimierung oder einen Fehler in meinem Modul.
• Bereitstellung alternativer / besserer Algorithmen / Implementierungen.

Die Antwort, die am vollständigsten / konstruktivsten ist, erhält das Kopfgeld.

Und schließlich das Modul selbst:

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)

Wenn Sie wollen das Rad nicht neu erfinden, verwenden Sie die Bibliothek sympy

pip install sympy

Verwenden Sie die Funktion sympy.ntheory.factorint

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

Sie können einige sehr große Zahlen berücksichtigen:

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

Es besteht keine Notwendigkeit zu berechnen smallprimesverwenden primesbelow, Verwendung smallprimesetdafür.

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

Teilen Sie Ihre primefactorsin zwei Funktionen für die Handhabung smallprimesund andere für pollard_brent, dies kann ein paar Iterationen sparen, da alle Potenzen von Smallprimes von n geteilt werden.

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

Indem Sie die verifizierten Ergebnisse von Pomerance, Selfridge, Wagstaff und Jaeschke berücksichtigen, können Sie die Wiederholungen reduzieren, bei isprimedenen der Miller-Rabin-Primalitätstest verwendet wird. Aus dem Wiki .

• wenn n <1.373.653 ist, reicht es aus, a = 2 und 3 zu testen;
• wenn n <9.080.191 ist, reicht es aus, a = 31 und 73 zu testen;
• wenn n <4,759,123,141 ist, reicht es aus, a = 2, 7 und 61 zu testen;
• wenn n <2.152.302.898.747 ist, reicht es aus, a = 2, 3, 5, 7 und 11 zu testen;
• wenn n <3,474,749,660,383 ist, reicht es aus, a = 2, 3, 5, 7, 11 und 13 zu testen;
• Wenn n <341,550,071,728,321 ist, reicht es aus, a = 2, 3, 5, 7, 11, 13 und 17 zu testen.

Bearbeiten 1 : Der Rückruf von wurde korrigiert if-else, um große Faktoren an Faktoren anzuhängen primefactors.

Selbst bei der aktuellen Version sind mehrere Stellen zu beachten.

1. Mach nicht checker*checkerjede Schleife, benutze s=ceil(sqrt(num))undchecher < s
2. checher sollte jedes Mal plus 2 sein und alle geraden Zahlen außer 2 ignorieren
3. Verwenden Sie divmodanstelle von %und//

Es gibt eine Python-Bibliothek mit einer Sammlung von Primalitätstests (einschließlich falscher Tests für das, was nicht zu tun ist). Es heißt Pyprimes . Ich dachte, es ist erwähnenswert für den Zweck der Nachwelt. Ich glaube nicht, dass es die von Ihnen erwähnten Algorithmen enthält.

Sie sollten wahrscheinlich eine Primerkennung durchführen, die Sie hier sehen könnten. Schneller Algorithmus zum Finden von Primzahlen?

Sie sollten jedoch den gesamten Blog lesen. Es gibt einige Algorithmen, die er zum Testen der Primalität auflistet.

Sie können bis zu einem Grenzwert faktorisieren und dann brent verwenden, um höhere Faktoren zu erhalten

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

Ich bin gerade auf einen Fehler in diesem Code gestoßen, als ich die Nummer berücksichtigt habe 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

Dieser Code-Ausschnitt:

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

sollte umgeschrieben werden

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

was bei realistischen Eingaben wahrscheinlich sowieso schneller funktioniert. Die Quadratwurzel ist langsam - im Grunde genommen das Äquivalent vieler Multiplikationen -, hat smallprimesnur ein paar Dutzend Mitglieder, und auf diese Weise entfernen wir die Berechnung von n ** .5aus der engen inneren Schleife, was sicherlich hilfreich ist, wenn Zahlen wie berücksichtigt werden 2**1427. Es gibt einfach keinen Grund zu berechnen sqrt(2**1427), sqrt(2**1426), sqrt(2**1425)etc. etc., wenn alles , was wir über Pflege ist „tut das [Quadrat der] checker überschreiten n“.

Der umgeschriebene Code löst bei großen Zahlen keine Ausnahmen aus und ist laut timeit(bei Beispieleingaben 2und 2**718 * 31) ungefähr doppelt so schnell .

Beachten Sie auch, dass isprime(2)das falsche Ergebnis zurückgegeben wird, aber dies ist in Ordnung, solange wir uns nicht darauf verlassen. IMHO sollten Sie das Intro dieser Funktion umschreiben als

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