Modul zur schnellen Primfaktorisierung

78

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)
orlp
quelle
1
@wheaties - dafür wäre das while checker*checker <= numda.
Amber
Sie könnten diesen Thread nützlich finden: stackoverflow.com/questions/1024640/calculating-phik-for-1kn/…
RBarryYoung
Warum sind solche Dinge in der Standardbibliothek nicht verfügbar? Wenn ich suche, finde ich nur eine Million Project Euler-Lösungsvorschläge und andere Leute, die auf Fehler hinweisen. Ist das nicht der Zweck von Bibliotheken und Fehlerberichten?
Endolith
@endolith Außerhalb von Dingen wie Prject Euler gibt es dafür nicht viel Verwendung. Sicherlich nicht genug, um es in die Standardbibliotheken aufzunehmen.
Orlp
@nightcracker: Es gibt keine praktische Verwendung für das Faktorisieren von Zahlen?
Endolith

Antworten:

65

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}
Oberst Panik
quelle
3
Vielen Dank für das Teilen, @Colonel Panic. Dies ist genau das , wonach ich gesucht habe: ein Integer-Faktorisierungsmodul in einer gut gepflegten Bibliothek anstelle von Codeausschnitten.
Programmierer
15

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.

Rozuur
quelle
Genießen Sie Ihre +100 (Sie sind der einzige, der seit dem Kopfgeld geantwortet hat). Ihr bigfactorsist jedoch schrecklich ineffizient, weil factors.extend(bigfactors(factor))Sie zu Big Factors zurückkehren, was einfach falsch ist (was ist, wenn Pollard-Brent den Faktor 234892 findet, möchten Sie das nicht noch einmal mit Pollard-Brent faktorisieren). Wenn Sie zu wechseln factors.extend(bigfactors(factor)), factors.extend(primefactors(factor, sort))ist es in Ordnung.
Orlp
Ein Primfaktor nennt BigFactors, dann ist klar, dass es bei den nächsten von Pollard-Brent erhaltenen Faktoren keine Potenz von Small Prime gibt.
Rozuur
Wenn es ineffizient wäre, hätte ich das nicht beantwortet. Sobald der Anruf von Primfaktoren zu Bigfaktoren geht, gibt es keinen Faktor in n, der weniger als 1000 ist, daher kann Pollard-Brent keine Zahl zurückgeben, deren Faktoren weniger als 1000 sind. Wenn es nicht klar ist, antworte so, dass ich meine Antwort mit weiteren Erklärungen bearbeite
Rozuur
Scheiße NVM natürlich. Wenn N keine Faktoren enthält, die von kleinen Primzahlen gefunden wurden, wird der Faktor F von N auch nicht>. <
orlp
Sie sollten auch smallprimeset anstelle von smallprime verwenden und smallprimeset von Miller-Rabin entfernen.
Rozuur
4

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//
Kabie
quelle
Ich muss checker * checker machen, weil num ständig abnimmt. Ich werde das Überspringen der geraden Zahlen jedoch implementieren. Der divmod verringert die Funktion erheblich (er berechnet das // in jeder Schleife, anstatt nur, wenn der Prüfer n teilt).
Orlp
@ Nacht, müssen Sie nur neu berechnen, swann immer Sie numdann ändern
John La Rooy
Stimmt, dachte mir, dass beim Herumspielen :) Es scheint schneller zu sein, sqrt neu zu berechnen als checker * checker.
Orlp
@nightcracker: Lassen Sie N=n*n+1, ceil(sqrt(N))kostet etwa 2 ~ 4 mal als n*n, numändert sich nicht so häufig.
Kabie
Kennen Sie einen Ceil / Floor-Quadrat-Algorithmus, weil int (num ** .5) + 1 zu viel zu sein scheint (zuerst in Float-Genauigkeit berechnen und dann abhacken).
Orlp
4

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.

Ehtesh Choudhury
quelle
3

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.

Milchmann
quelle
1

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))
Antoni Gual Via
quelle
0

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