Schnellste Semiprime-Faktorisierung

Schreiben Sie ein Programm, um eine Semi-Primzahl in kürzester Zeit zu zerlegen.

Verwenden Sie zu Testzwecken Folgendes: 38! +1 (523022617466601111760007224100074291200000001)

Es ist gleich: 14029308060317546154181 × 37280713718589679646221

Python (mit PyPy JIT v1.9) ~ 1.9s

Verwenden eines quadratischen Siebs mit mehreren Polynomen . Ich nahm dies als eine Code-Herausforderung an und entschied mich daher, keine externen Bibliotheken zu verwenden (außer der Standardfunktion log, nehme ich an). Beim Timing sollte das PyPy-JIT verwendet werden, da es zu 4-5-mal schnelleren Timings als das von cPython führt .

Update (29.07.2013):
Seit dem ursprünglichen Posting habe ich einige geringfügige, aber signifikante Änderungen vorgenommen, die die Gesamtgeschwindigkeit um den Faktor 2,5 erhöhen.

Update (27.08.2014):
Da dieser Beitrag immer noch Beachtung findet, wurde die my_math.pyKorrektur von zwei Fehlern für alle Nutzer aktualisiert :

• isqrtwar fehlerhaft und erzeugte manchmal eine falsche Ausgabe für Werte, die einem perfekten Quadrat sehr nahe kommen. Dies wurde korrigiert und die Leistung durch die Verwendung eines viel besseren Samens gesteigert.
• is_primewurde aktualisiert. Mein vorheriger Versuch, perfekte quadratische 2-Sprps zu entfernen, war bestenfalls halbherzig. Ich habe eine 3-Sprp-Prüfung hinzugefügt - eine von Mathmatica verwendete Technik - um sicherzustellen, dass der getestete Wert quadratfrei ist.

Update (24.11.2014):
Wenn am Ende der Berechnung keine nicht-trivialen Kongruenzen gefunden werden, siebt das Programm nun zusätzliche Polynome. Dies wurde zuvor im Code als markiert TODO.

mpqs.py

from my_math import *
from math import log
from time import clock
from argparse import ArgumentParser

# Multiple Polynomial Quadratic Sieve
def mpqs(n, verbose=False):
  if verbose:
    time1 = clock()

  root_n = isqrt(n)
  root_2n = isqrt(n+n)

  # formula chosen by experimentation
  # seems to be close to optimal for n < 10^50
  bound = int(5 * log(n, 10)**2)

  prime = []
  mod_root = []
  log_p = []
  num_prime = 0

  # find a number of small primes for which n is a quadratic residue
  p = 2
  while p < bound or num_prime < 3:

    # legendre (n|p) is only defined for odd p
    if p > 2:
      leg = legendre(n, p)
    else:
      leg = n & 1

    if leg == 1:
      prime += [p]
      mod_root += [int(mod_sqrt(n, p))]
      log_p += [log(p, 10)]
      num_prime += 1
    elif leg == 0:
      if verbose:
        print 'trial division found factors:'
        print p, 'x', n/p
      return p

    p = next_prime(p)

  # size of the sieve
  x_max = len(prime)*60

  # maximum value on the sieved range
  m_val = (x_max * root_2n) >> 1

  # fudging the threshold down a bit makes it easier to find powers of primes as factors
  # as well as partial-partial relationships, but it also makes the smoothness check slower.
  # there's a happy medium somewhere, depending on how efficient the smoothness check is
  thresh = log(m_val, 10) * 0.735

  # skip small primes. they contribute very little to the log sum
  # and add a lot of unnecessary entries to the table
  # instead, fudge the threshold down a bit, assuming ~1/4 of them pass
  min_prime = int(thresh*3)
  fudge = sum(log_p[i] for i,p in enumerate(prime) if p < min_prime)/4
  thresh -= fudge

  if verbose:
    print 'smoothness bound:', bound
    print 'sieve size:', x_max
    print 'log threshold:', thresh
    print 'skipping primes less than:', min_prime

  smooth = []
  used_prime = set()
  partial = {}
  num_smooth = 0
  num_used_prime = 0
  num_partial = 0
  num_poly = 0
  root_A = isqrt(root_2n / x_max)

  if verbose:
    print 'sieving for smooths...'
  while True:
    # find an integer value A such that:
    # A is =~ sqrt(2*n) / x_max
    # A is a perfect square
    # sqrt(A) is prime, and n is a quadratic residue mod sqrt(A)
    while True:
      root_A = next_prime(root_A)
      leg = legendre(n, root_A)
      if leg == 1:
        break
      elif leg == 0:
        if verbose:
          print 'dumb luck found factors:'
          print root_A, 'x', n/root_A
        return root_A

    A = root_A * root_A

    # solve for an adequate B
    # B*B is a quadratic residue mod n, such that B*B-A*C = n
    # this is unsolvable if n is not a quadratic residue mod sqrt(A)
    b = mod_sqrt(n, root_A)
    B = (b + (n - b*b) * mod_inv(b + b, root_A))%A

    # B*B-A*C = n <=> C = (B*B-n)/A
    C = (B*B - n) / A

    num_poly += 1

    # sieve for prime factors
    sums = [0.0]*(2*x_max)
    i = 0
    for p in prime:
      if p < min_prime:
        i += 1
        continue
      logp = log_p[i]

      inv_A = mod_inv(A, p)
      # modular root of the quadratic
      a = int(((mod_root[i] - B) * inv_A)%p)
      b = int(((p - mod_root[i] - B) * inv_A)%p)

      k = 0
      while k < x_max:
        if k+a < x_max:
          sums[k+a] += logp
        if k+b < x_max:
          sums[k+b] += logp
        if k:
          sums[k-a+x_max] += logp
          sums[k-b+x_max] += logp

        k += p
      i += 1

    # check for smooths
    i = 0
    for v in sums:
      if v > thresh:
        x = x_max-i if i > x_max else i
        vec = set()
        sqr = []
        # because B*B-n = A*C
        # (A*x+B)^2 - n = A*A*x*x+2*A*B*x + B*B - n
        #               = A*(A*x*x+2*B*x+C)
        # gives the congruency
        # (A*x+B)^2 = A*(A*x*x+2*B*x+C) (mod n)
        # because A is chosen to be square, it doesn't need to be sieved
        val = sieve_val = A*x*x + 2*B*x + C

        if sieve_val < 0:
          vec = set([-1])
          sieve_val = -sieve_val

        for p in prime:
          while sieve_val%p == 0:
            if p in vec:
              # keep track of perfect square factors
              # to avoid taking the sqrt of a gigantic number at the end
              sqr += [p]
            vec ^= set([p])
            sieve_val = int(sieve_val / p)

        if sieve_val == 1:
          # smooth
          smooth += [(vec, (sqr, (A*x+B), root_A))]
          used_prime |= vec
        elif sieve_val in partial:
          # combine two partials to make a (xor) smooth
          # that is, every prime factor with an odd power is in our factor base
          pair_vec, pair_vals = partial[sieve_val]
          sqr += list(vec & pair_vec) + [sieve_val]
          vec ^= pair_vec
          smooth += [(vec, (sqr + pair_vals[0], (A*x+B)*pair_vals[1], root_A*pair_vals[2]))]
          used_prime |= vec
          num_partial += 1
        else:
          # save partial for later pairing
          partial[sieve_val] = (vec, (sqr, A*x+B, root_A))
      i += 1

    num_smooth = len(smooth)
    num_used_prime = len(used_prime)

    if verbose:
      print 100 * num_smooth / num_prime, 'percent complete\r',

    if num_smooth > num_used_prime:
      if verbose:
        print '%d polynomials sieved (%d values)'%(num_poly, num_poly*x_max*2)
        print 'found %d smooths (%d from partials) in %f seconds'%(num_smooth, num_partial, clock()-time1)
        print 'solving for non-trivial congruencies...'

      used_prime_list = sorted(list(used_prime))

      # set up bit fields for gaussian elimination
      masks = []
      mask = 1
      bit_fields = [0]*num_used_prime
      for vec, vals in smooth:
        masks += [mask]
        i = 0
        for p in used_prime_list:
          if p in vec: bit_fields[i] |= mask
          i += 1
        mask <<= 1

      # row echelon form
      col_offset = 0
      null_cols = []
      for col in xrange(num_smooth):
        pivot = col-col_offset == num_used_prime or bit_fields[col-col_offset] & masks[col] == 0
        for row in xrange(col+1-col_offset, num_used_prime):
          if bit_fields[row] & masks[col]:
            if pivot:
              bit_fields[col-col_offset], bit_fields[row] = bit_fields[row], bit_fields[col-col_offset]
              pivot = False
            else:
              bit_fields[row] ^= bit_fields[col-col_offset]
        if pivot:
          null_cols += [col]
          col_offset += 1

      # reduced row echelon form
      for row in xrange(num_used_prime):
        # lowest set bit
        mask = bit_fields[row] & -bit_fields[row]
        for up_row in xrange(row):
          if bit_fields[up_row] & mask:
            bit_fields[up_row] ^= bit_fields[row]

      # check for non-trivial congruencies
      for col in null_cols:
        all_vec, (lh, rh, rA) = smooth[col]
        lhs = lh   # sieved values (left hand side)
        rhs = [rh] # sieved values - n (right hand side)
        rAs = [rA] # root_As (cofactor of lhs)
        i = 0
        for field in bit_fields:
          if field & masks[col]:
            vec, (lh, rh, rA) = smooth[i]
            lhs += list(all_vec & vec) + lh
            all_vec ^= vec
            rhs += [rh]
            rAs += [rA]
          i += 1

        factor = gcd(list_prod(rAs)*list_prod(lhs) - list_prod(rhs), n)
        if factor != 1 and factor != n:
          break
      else:
        if verbose:
          print 'none found.'
        continue
      break

  if verbose:
    print 'factors found:'
    print factor, 'x', n/factor
    print 'time elapsed: %f seconds'%(clock()-time1)
  return factor

if __name__ == "__main__":
  parser =ArgumentParser(description='Uses a MPQS to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long,
      help='the composite number to factor')
  parser.add_argument('--verbose', dest='verbose', action='store_true',
      help="enable verbose output")
  args = parser.parse_args()

  if args.verbose:
    mpqs(args.composite, args.verbose)
  else:
    time1 = clock()
    print mpqs(args.composite)
    print 'time elapsed: %f seconds'%(clock()-time1)

my_math.py

# divide and conquer list product
def list_prod(a):
  size = len(a)
  if size == 1:
    return a[0]
  return list_prod(a[:size>>1]) * list_prod(a[size>>1:])

# greatest common divisor of a and b
def gcd(a, b):
  while b:
    a, b = b, a%b
  return a

# modular inverse of a mod m
def mod_inv(a, m):
  a = int(a%m)
  x, u = 0, 1
  while a:
    x, u = u, x - (m/a)*u
    m, a = a, m%a
  return x

# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow(a, (m-1) >> 1, m)

# modular sqrt(n) mod p
# p must be prime
def mod_sqrt(n, p):
  a = n%p
  if p%4 == 3:
    return pow(a, (p+1) >> 2, p)
  elif p%8 == 5:
    v = pow(a << 1, (p-5) >> 3, p)
    i = ((a*v*v << 1) % p) - 1
    return (a*v*i)%p
  elif p%8 == 1:
    # Shank's method
    q = p-1
    e = 0
    while q&1 == 0:
      e += 1
      q >>= 1

    n = 2
    while legendre(n, p) != p-1:
      n += 1

    w = pow(a, q, p)
    x = pow(a, (q+1) >> 1, p)
    y = pow(n, q, p)
    r = e
    while True:
      if w == 1:
        return x

      v = w
      k = 0
      while v != 1 and k+1 < r:
        v = (v*v)%p
        k += 1

      if k == 0:
        return x

      d = pow(y, 1 << (r-k-1), p)
      x = (x*d)%p
      y = (d*d)%p
      w = (w*y)%p
      r = k
  else: # p == 2
    return a

#integer sqrt of n
def isqrt(n):
  c = n*4/3
  d = c.bit_length()

  a = d>>1
  if d&1:
    x = 1 << a
    y = (x + (n >> a)) >> 1
  else:
    x = (3 << a) >> 2
    y = (x + (c >> a)) >> 1

  if x != y:
    x = y
    y = (x + n/x) >> 1
    while y < x:
      x = y
      y = (x + n/x) >> 1
  return x

# strong probable prime
def is_sprp(n, b=2):
  if n < 2: return False
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in xrange(1, s):
    x = (x * x)%n
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False

# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  P = 1
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = 0
  while s:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t:
    if t&1:
      # U, V of n+1
      U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
      q = (q * Q)%n
      t -= 1
    else:
      # U, V of n*2
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r:
    U, V = (U * V)%n, (V * V - 2 * q)%n
    q = (q * q)%n
    r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0

# primes less than 212
small_primes = set([
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
   31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
   73, 79, 83, 89, 97,101,103,107,109,113,
  127,131,137,139,149,151,157,163,167,173,
  179,181,191,193,197,199,211])

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41,
   43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
   89, 97,101,103,107,109,113,121,127,131,
  137,139,143,149,151,157,163,167,169,173,
  179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

max_int = 2147483647

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    return n in small_primes

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    i = 211
    while i*i < n:
      for o in offsets:
        i += o
        if n%i == 0:
          return False
    return True

  # Baillie-PSW
  # this is technically a probabalistic test, but there are no known pseudoprimes
  if not is_sprp(n, 2): return False

  # idea shamelessly stolen from Mathmatica
  # if n is a 2-sprp and a 3-sprp, n is necessarily square-free
  if not is_sprp(n, 3): return False

  a = 5
  s = 2
  # if n is a perfect square, this will never terminate
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2
  # first odd larger than n
  n = (n + 1) | 1
  if n < 212:
    while True:
      if n in small_primes:
        return n
      n += 2

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  s = 0
  e = 47
  m = 24
  while m != e:
    if indices[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1

  i = int(n + (indices[m] - x))
  # adjust offsets
  offs = offsets[m:] + offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o

Beispiel-E / A:

$ pypy mpqs.py --verbose 94968915845307373740134800567566911
smoothness bound: 6117
sieve size: 24360
log threshold: 14.3081031579
skipping primes less than: 47
sieving for smooths...
144 polynomials sieved (7015680 values)
found 405 smooths (168 from partials) in 0.513794 seconds
solving for non-trivial congruencies...
factors found:
216366620575959221 x 438925910071081891
time elapsed: 0.685765 seconds

$ pypy mpqs.py --verbose 523022617466601111760007224100074291200000001
smoothness bound: 9998
sieve size: 37440
log threshold: 15.2376302725
skipping primes less than: 59
sieving for smooths...
428 polynomials sieved (32048640 values)
found 617 smooths (272 from partials) in 1.912131 seconds
solving for non-trivial congruencies...
factors found:
14029308060317546154181 x 37280713718589679646221
time elapsed: 2.064387 seconds

Hinweis: --verboseWenn Sie diese Option nicht verwenden , erhalten Sie ein etwas besseres Timing:

$ pypy mpqs.py 94968915845307373740134800567566911
216366620575959221
time elapsed: 0.630235 seconds

$ pypy mpqs.py 523022617466601111760007224100074291200000001
14029308060317546154181
time elapsed: 1.886068 seconds

Grundlegendes Konzept

Im Allgemeinen basiert ein quadratisches Sieb auf der folgenden Beobachtung: Jedes ungerade Komposit n kann dargestellt werden als:

Dies ist nicht sehr schwer zu bestätigen. Da n ungerade ist, muss der Abstand zwischen zwei Cofaktoren von n gerade 2d sein , wobei x der Mittelpunkt zwischen ihnen ist. Darüber hinaus gilt die gleiche Beziehung für ein beliebiges Vielfaches von n

Es ist zu beachten, dass wenn ein solches x und d gefunden werden kann, dies sofort zu einem (nicht notwendigerweise primalen) Faktor von n führt , da x + d und x - d beide n durch Definition teilen . Diese Beziehung kann weiter geschwächt werden - als Folge möglicher unbedeutender Kongruenzen - auf die folgende Form:

Wenn wir also zwei perfekte Quadrate finden, die äquivalent zu mod n sind , ist es ziemlich wahrscheinlich, dass wir direkt einen Faktor von n a la gcd (x ± d, n) erzeugen können . Scheint ziemlich einfach, oder?

Außer es ist nicht. Wenn wir beabsichtigen, eine umfassende Suche über alle möglichen x durchzuführen , müssten wir den gesamten Bereich von [ √ n , √ ( 2n ) ] durchsuchen , der geringfügig kleiner ist als die vollständige Testdivision, der jedoch bei is_squarejeder Iteration eine teure Operation erfordert Bestätigen Sie den Wert von d . Sofern nicht im Voraus bekannt ist, dass n Faktoren in der Nähe von √ n aufweist , ist die Teilung des Versuchs wahrscheinlich schneller.

Vielleicht können wir diese Beziehung noch mehr schwächen. Angenommen, wir haben ein x gewählt , so dass für

Eine vollständige Primfaktorisierung von y ist ohne weiteres bekannt. Wenn wir genug solche Beziehungen hätten, sollten wir in der Lage sein, ein adäquates d zu konstruieren , wenn wir eine Anzahl von y so wählen, dass ihr Produkt ein perfektes Quadrat ist; Das heißt, alle Primfaktoren werden gerade oft verwendet. In der Tat, wenn wir mehr haben , so y als die Gesamtzahl der einzelnen Primfaktoren sie enthalten, eine Lösung existieren garantiert; Es wird ein lineares Gleichungssystem. Es stellt sich nun die Frage, wie wir uns für ein solches x entschieden haben . Hier kommt das Sieben ins Spiel.

Das Sieb

Betrachten Sie das Polynom:

Dann gilt für jede Primzahl p und ganze Zahl k :

Dies bedeutet, dass Sie nach dem Auflösen nach den Wurzeln des Polynoms mod p - das heißt, Sie haben ein x so gefunden, dass y (x) ≡ 0 (mod p) , ergo y durch p teilbar ist - dann haben Sie eine unendliche Zahl gefunden von solchen x . Auf diese Weise können Sie über einen Bereich von x sieben, wobei Sie kleine Primfaktoren von y identifizieren und hoffentlich einige finden, für die alle Primfaktoren klein sind. Solche Zahlen sind als k-glatt bekannt , wobei k der größte verwendete Primfaktor ist.

Bei diesem Ansatz gibt es jedoch einige Probleme. Nicht alle Werte von x sind angemessen. Tatsächlich gibt es nur sehr wenige von ihnen, die um √ n zentriert sind . Kleinere Werte werden (aufgrund des -n- Terms) größtenteils negativ , und größere Werte werden zu groß, so dass es unwahrscheinlich ist, dass ihre Primfaktorisierung nur aus kleinen Primzahlen besteht. Es gibt eine Reihe solcher x , aber es ist sehr unwahrscheinlich, dass Sie genügend Glättungen finden, um zu einer Faktorisierung zu führen, es sei denn, das zu faktorisierende Composite ist sehr klein. Und so wird es für ein größeres n notwendig, mehrere Polynome einer gegebenen Form zu sieben .

Mehrere Polynome

Wir brauchen also mehr Polynome zum Sieben? Wie wäre es damit:

Das wird klappen Beachten Sie, dass A und B buchstäblich ein ganzzahliger Wert sein können und die Mathematik weiterhin gilt. Wir müssen nur ein paar zufällige Werte auswählen, nach der Wurzel des Polynoms suchen und die Werte nahe Null sieben. An diesem Punkt könnten wir es einfach gut genug nennen: Wenn Sie genug Steine ​​in zufällige Richtungen werfen, müssen Sie früher oder später ein Fenster zerbrechen.

Es gibt aber auch ein Problem damit. Wenn die Steigung des Polynoms am x-Achsenabschnitt groß ist, was dann der Fall ist, wenn es nicht relativ flach ist, gibt es nur wenige geeignete Werte zum Sieben pro Polynom. Es wird funktionieren, aber am Ende werden Sie eine ganze Reihe von Polynomen sieben, bevor Sie das bekommen, was Sie brauchen. Können wir es besser machen?

Wir können es besser machen. Eine Beobachtung als Ergebnis von Montgomery lautet wie folgt: Wenn A und B so gewählt sind, dass es ein gewisses C gibt, das zufriedenstellend ist

dann kann das gesamte Polynom als umgeschrieben werden

Wenn außerdem A als perfektes Quadrat gewählt wird, kann der führende A- Term beim Sieben vernachlässigt werden, was zu viel kleineren Werten und einer viel flacheren Kurve führt. Damit eine solche Lösung existiert, muss n ein quadratischer Rest mod √ A sein , der sofort durch Berechnung des Legendre-Symbols erkannt werden kann :
( n | √ A ) = 1 . Beachten Sie, dass zur Lösung von B eine vollständige Primfaktorisierung von √A bekannt sein muss (um die modulare Quadratwurzel √n (mod √A) zu erhalten ), weshalb √A normalerweise als Primzahl gewählt wird.

Es kann dann gezeigt werden, dass wenn , dann für alle Werte von x ∈ [ -M, M ] :

Und jetzt haben wir endlich alle notwendigen Komponenten, um unser Sieb zu implementieren. Oder wir?

Mächte der Primzahlen als Faktoren

Unser Sieb hat, wie oben beschrieben, einen Hauptfehler. Es kann identifizieren, welche Werte von x zu einem durch p teilbaren y führen , es kann jedoch nicht identifizieren, ob dieses y durch eine Potenz von p teilbar ist oder nicht . Um dies festzustellen, müssten wir den zu siebenden Wert probeweise teilen, bis er nicht mehr durch p teilbar ist . Wir schienen einen Stillstand erreicht zu haben: Der ganze Punkt des Siebs war so, dass wir das nicht tun mussten. Zeit, das Spielbuch zu überprüfen.

Das sieht ziemlich nützlich aus. Wenn die Summe der ln aller kleinen Primfaktoren von y nahe am erwarteten Wert von ln (y) liegt , ist es fast selbstverständlich, dass y keine anderen Faktoren hat. Wenn wir den erwarteten Wert ein wenig nach unten korrigieren, können wir außerdem Werte als glatt identifizieren, die mehrere Potenzen von Primzahlen als Faktoren haben. Auf diese Weise können wir das Sieb als "Vorsieb" -Prozess verwenden und nur die Werte berücksichtigen, die wahrscheinlich glatt sind.

Dies hat auch einige andere Vorteile. Beachten Sie, dass kleine Primzahlen sehr wenig zur beitragen ln Summe, aber doch erfordern sie die meiste Zeit Sieb. Das Sieben des Werts 3 erfordert mehr Zeit als 11, 13, 17, 19 und 23 zusammen . Stattdessen können wir einfach die ersten Primzahlen überspringen und den Schwellenwert entsprechend anpassen, vorausgesetzt, ein bestimmter Prozentsatz von ihnen wäre bestanden.

Ein weiteres Ergebnis ist, dass eine Reihe von Werten "durchrutschen" kann, die größtenteils glatt sind, aber einen einzigen großen Cofaktor enthalten. Wir könnten diese Werte einfach verwerfen, aber nehmen wir an, wir haben einen anderen größtenteils glatten Wert mit genau demselben Cofaktor gefunden. Wir können diese beiden Werte dann verwenden, um ein verwendbares y zu konstruieren ; Da ihr Produkt diesen großen Cofaktor im Quadrat enthält, muss er nicht mehr berücksichtigt werden.

Alles zusammen

Das Letzte, was wir tun müssen, ist, diese Werte von y zu verwenden, um ein angemessenes x und d zu konstruieren . Angenommen, wir betrachten nur die nichtquadratischen Faktoren von y , dh die Primfaktoren einer ungeraden Potenz. Dann kann jedes y folgendermaßen ausgedrückt werden:

was in der Matrixform ausgedrückt werden kann:

Das Problem wird dann, einen Vektor v zu finden, so dass vM = ⦳ (mod 2) , wobei ⦳ der Nullvektor ist. Das heißt, nach dem linken Nullraum von M zu lösen . Dies kann auf verschiedene Arten geschehen, wobei die einfachste darin besteht, die Gaußsche Eliminierung an M ^T durchzuführen und die Zeilenadditionsoperation durch eine Zeilenxor- Operation zu ersetzen . Dies führt zu einer Anzahl von Nullraum-Basisvektoren, von denen jede Kombination eine gültige Lösung ergibt.

Die Konstruktion von x ist ziemlich einfach. Es ist einfach das Produkt von Ax + B für jedes der verwendeten y . Die Konstruktion von d ist etwas komplizierter. Wenn wir das Produkt von allen y nehmen , erhalten wir einen Wert mit Zehntausenden, wenn nicht Hunderttausenden von Ziffern, für die wir die Quadratwurzel finden müssen. Diese Berechnung ist unpraktisch teuer. Stattdessen können wir die geraden Potenzen von Primzahlen während des Siebprozesses verfolgen und dann and und xor- Operationen für die Vektoren von nicht quadratischen Faktoren verwenden, um die Quadratwurzel zu rekonstruieren.

Ich habe anscheinend die maximale Zeichenanzahl von 30000 erreicht. Ahh gut, ich nehme an, das ist gut genug.

Nun, Ihre 38! +1 haben mein PHP-Skript gebrochen, nicht sicher warum. Tatsächlich bricht jedes Semi-Prime mit mehr als 16 Ziffern mein Skript.

Unter Verwendung von 8980935344490257 (86028157 * 104395301) verwaltete mein Skript auf meinem Heimcomputer (2,61 GHz AMD Phenom 9950) eine Zeit von 25,963 Sekunden . Viel schneller als mein Arbeitscomputer, der fast 31 Sekunden bei 2,93 GHz Core 2 Duo war.

PHP - 757 Zeichen inkl. neue Zeilen

<?php
function getTime() {
    $t = explode( ' ', microtime() );
    $t = $t[1] + $t[0];
    return $t;
}
function isDecimal($val){ return is_numeric($val) && floor($val) != $val;}
$start = getTime();
$semi_prime = 8980935344490257;
$slice      = round(strlen($semi_prime)/2);
$max        = (pow(10, ($slice))-1);
$i          = 3;
echo "\nFactoring the semi-prime:\n$semi_prime\n\n";

while ($i < $max) {
    $sec_factor = ($semi_prime/$i);
    if (isDecimal($sec_factor) != 1) {
        $mod_f = bcmod($i, 1);
        $mod_s = bcmod($sec_factor, 1);
        if ($mod_f == 0 && $mod_s == 0) {
            echo "First factor = $i\n";
            echo "Second factor = $sec_factor\n";
            $end=getTime();
            $xtime=round($end-$start,4).' seconds';
            echo "\n$xtime\n";
            exit();
        }
    }
    $i += 2;
}
?>

Es würde mich interessieren, diesen Algorithmus in c oder einer anderen kompilierten Sprache zu sehen.