Faktorisasi semiprime tercepat

Tulis sebuah program untuk memfaktorkan bilangan semi-prima dalam waktu singkat.

Untuk tujuan pengujian, gunakan ini: 38! +1 (523022617466601111760007224100074291200000001)

Itu sama dengan: 14029308060317546154181 × 37280713718589679646221

Python (w / PyPy JIT v1.9) ~ 1.9s

Menggunakan Multiple Polynomial Quadratic Sieve . Saya menganggap ini sebagai tantangan kode, jadi saya memilih untuk tidak menggunakan perpustakaan eksternal (selain logfungsi standar , saya kira). Saat menghitung waktu, JIT PyPy harus digunakan, karena menghasilkan timing 4-5 kali lebih cepat daripada cPython .

Pembaruan (2013-07-29):
Sejak awalnya memposting, saya telah membuat beberapa perubahan kecil, tetapi signifikan yang meningkatkan kecepatan keseluruhan dengan faktor sekitar 2,5x.

Pembaruan (2014-08-27):
Karena posting ini masih menerima perhatian, saya telah memperbarui my_math.pymemperbaiki dua kesalahan, untuk siapa saja yang mungkin menggunakannya:

• isqrtrusak, terkadang menghasilkan output yang salah untuk nilai yang sangat dekat dengan kuadrat sempurna. Ini telah diperbaiki, dan kinerjanya meningkat dengan menggunakan benih yang jauh lebih baik.
• is_primetelah diperbarui. Upaya saya sebelumnya untuk menghapus 2-sprps kuadrat sempurna adalah setengah hati, paling-paling. Saya telah menambahkan cek 3-sprp - teknik yang digunakan oleh Mathmatica - untuk memastikan bahwa nilai yang diuji bebas persegi.

Pembaruan (2014-11-24):
Jika pada akhir perhitungan tidak ditemukan kongruensi non-sepele, progam sekarang menyaring polinomial tambahan. Ini sebelumnya ditandai dalam kode sebagai 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

Sampel I / O:

$ 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

Catatan: tidak menggunakan --verboseopsi akan memberikan timing yang sedikit lebih baik:

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

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

Konsep dasar

Secara umum, ayakan kuadrat didasarkan pada pengamatan berikut: setiap komposit n aneh dapat direpresentasikan sebagai:

Ini tidak terlalu sulit untuk dikonfirmasi. Karena n ganjil, jarak antara dua kofaktor n harus genap 2d , di mana x adalah titik tengah di antara mereka. Selain itu, hubungan yang sama berlaku untuk kelipatan n

Perhatikan bahwa jika ada x dan d tersebut dapat ditemukan, itu akan segera menghasilkan faktor (tidak harus prima) dari n , karena x + d dan x - d keduanya membagi n berdasarkan definisi. Hubungan ini selanjutnya dapat dilemahkan - dengan konsekuensi memungkinkan kongruensi sepele yang potensial - ke bentuk berikut:

Jadi secara umum, jika kita dapat menemukan dua kuadrat sempurna yang setara dengan mod n , maka sangat mungkin bahwa kita dapat langsung menghasilkan faktor n a la gcd (x ± d, n) . Tampaknya cukup sederhana, bukan?

Kecuali tidak. Jika kami bermaksud melakukan pencarian lengkap atas semua kemungkinan x , kita perlu mencari seluruh rentang dari [ √ n , √ ( 2n ) ], yang sedikit lebih kecil dari divisi percobaan penuh, tetapi juga membutuhkan is_squareoperasi mahal setiap iterasi untuk konfirmasi nilai d . Kecuali diketahui sebelumnya bahwa n memiliki faktor yang sangat dekat √ n , divisi percobaan kemungkinan akan lebih cepat.

Mungkin kita bisa melemahkan hubungan ini lebih jauh lagi. Misalkan kita memilih x , sehingga untuk

faktorisasi utama penuh y sudah diketahui. Jika kita memiliki cukup hubungan seperti itu, kita harus dapat membangun sebuah d yang memadai , jika kita memilih sejumlah y sedemikian rupa sehingga produk mereka adalah kuadrat sempurna; yaitu, semua faktor prima digunakan beberapa kali. Bahkan, jika kita memiliki lebih banyak y daripada jumlah total faktor prima unik yang dikandungnya, sebuah solusi dijamin ada; Ini menjadi sistem persamaan linear. Pertanyaannya sekarang adalah, bagaimana kita memilih x seperti itu ? Di situlah pengayakan berperan.

Saringan

Pertimbangkan polinomialnya:

Kemudian untuk setiap prime p dan integer k , yang berikut ini benar:

Ini berarti bahwa setelah menyelesaikan akar dari polinomial mod p - yaitu, Anda telah menemukan x sehingga y (x) ≡ 0 (mod p) , ergo y dapat dibagi oleh p - maka Anda telah menemukan angka tak hingga dari x tersebut . Dengan cara ini, Anda dapat menyaring rentang x , mengidentifikasi faktor prima kecil y , semoga menemukan beberapa yang semua faktor prima kecil. Angka-angka seperti itu dikenal sebagai k-smooth , di mana k adalah faktor utama terbesar yang digunakan.

Namun, ada beberapa masalah dengan pendekatan ini. Tidak semua nilai x memadai, pada kenyataannya, hanya ada sangat sedikit dari mereka, berpusat di sekitar √ n . Nilai yang lebih kecil akan menjadi sebagian besar negatif (karena istilah -n ), dan nilai yang lebih besar akan menjadi terlalu besar, sehingga tidak mungkin bahwa faktorisasi utama mereka hanya terdiri dari bilangan prima kecil. Akan ada sejumlah x seperti itu , tetapi kecuali jika komposit yang Anda anjak sangat kecil, sangat tidak mungkin bahwa Anda akan menemukan cukup smooth untuk menghasilkan faktorisasi. Jadi, untuk n yang lebih besar , menjadi perlu untuk menyaring beberapa polinomial dari bentuk yang diberikan.

Banyak Polinomial

Jadi kita perlu lebih banyak polinomial untuk diayak? Bagaimana dengan ini:

Itu akan berhasil. Perhatikan bahwa A dan B secara harfiah bisa berupa nilai integer apa pun, dan matematika masih berlaku. Yang perlu kita lakukan adalah memilih beberapa nilai acak, menyelesaikan untuk akar polinomial, dan menyaring nilai mendekati nol. Pada titik ini kita bisa menyebutnya cukup baik: jika Anda melempar cukup banyak batu ke arah acak, Anda pasti akan merusak jendela cepat atau lambat.

Kecuali, ada masalah dengan itu juga. Jika kemiringan polinomial besar pada intersep x, yang akan terjadi jika tidak relatif datar, hanya akan ada beberapa nilai yang sesuai untuk ayakan per polinomial. Ini akan berhasil, tetapi pada akhirnya Anda akan menyaring banyak polinomial sebelum mendapatkan yang Anda butuhkan. Bisakah kita berbuat lebih baik?

Kita bisa melakukan yang lebih baik. Pengamatan, sebagai hasil dari Montgomery adalah sebagai berikut: jika A dan B dipilih sedemikian rupa sehingga ada beberapa C yang memuaskan

maka seluruh polinomial dapat ditulis ulang sebagai

Lebih lanjut, jika A dipilih menjadi kuadrat sempurna, suku A yang memimpin dapat diabaikan sementara pengayakan, menghasilkan nilai yang jauh lebih kecil, dan kurva yang jauh lebih rata. Agar solusi seperti itu ada, n harus merupakan mod residu kuadrat √ A , yang dapat segera diketahui dengan menghitung simbol Legendre : ( n | √A ) = 1 . Perhatikan bahwa untuk menyelesaikan B , faktorisasi lengkap lengkap √A perlu diketahui (untuk mengambil akar kuadrat modular √ n (mod √A) ), itulah sebabnya √A biasanya dipilih untuk menjadi prima.

Kemudian dapat ditunjukkan bahwa jika , maka untuk semua nilai x ∈ [ -M, M ] :

Dan sekarang, akhirnya, kita memiliki semua komponen yang diperlukan untuk menerapkan saringan kita. Atau apakah kita?

Kekuatan Primes sebagai Faktor

Saringan kami, seperti dijelaskan di atas, memiliki satu kelemahan utama. Hal ini dapat mengidentifikasi nilai-nilai x akan menghasilkan y dibagi dengan p , tetapi tidak dapat mengidentifikasi apakah atau tidak ini y habis dibagi oleh kekuatan dari p . Untuk menentukan itu, kita perlu melakukan pembagian percobaan pada nilai yang akan disaring, sampai tidak lagi dapat dibagi oleh p . Kami sepertinya telah mencapai jalan buntu: inti dari ayakan itu adalah agar kami tidak harus melakukan itu. Saatnya memeriksa buku pedoman.

Itu terlihat sangat berguna. Jika jumlah ln dari semua faktor prima kecil y mendekati nilai yang diharapkan dari ln (y) , maka hampir diberikan bahwa y tidak memiliki faktor lain. Selain itu, jika kita menyesuaikan sedikit nilai yang diharapkan, kita juga dapat mengidentifikasi nilai-nilai yang halus yang memiliki beberapa kekuatan bilangan prima sebagai faktor. Dengan cara ini, kita dapat menggunakan saringan sebagai proses 'pra-penyaringan', dan hanya memfaktorkan nilai-nilai yang cenderung lancar.

Ini memiliki beberapa keunggulan lain juga. Perhatikan bahwa bilangan prima kecil berkontribusi sangat sedikit ke ln sum, tetapi belum mereka membutuhkan waktu yang paling saringan. Penyaringan nilai 3 membutuhkan lebih banyak waktu daripada 11, 13, 17, 19, dan 23 digabungkan . Sebagai gantinya, kita bisa melewatkan beberapa bilangan prima pertama, dan menyesuaikan ambang batas sesuai, dengan asumsi persentase tertentu dari mereka akan berlalu.

Hasil lain, adalah bahwa sejumlah nilai akan diizinkan 'lolos', yang sebagian besar mulus, tetapi mengandung satu kofaktor besar. Kami hanya bisa membuang nilai-nilai ini, tetapi anggaplah kami menemukan nilai lain yang sebagian besar mulus, dengan kofaktor yang persis sama. Kita kemudian dapat menggunakan kedua nilai ini untuk membangun y yang dapat digunakan ; karena produk mereka akan mengandung kofaktor besar ini kuadrat, tidak perlu lagi dipertimbangkan.

Menyatukan semuanya

Hal terakhir yang perlu kita lakukan adalah menggunakan nilai-nilai dari y construct x dan d yang memadai . Misalkan kita hanya mempertimbangkan faktor-faktor non-kuadrat dari y , yaitu faktor utama dari kekuatan ganjil. Kemudian, masing-masing y dapat diekspresikan dengan cara berikut:

yang dapat diekspresikan dalam bentuk matriks:

Masalahnya kemudian menjadi menemukan vektor v sedemikian rupa sehingga vM = ⦳ (mod 2) , di mana ⦳ adalah vektor nol. Artinya, untuk memecahkan ruang nol sebelah kiri M . Hal ini dapat dilakukan dengan beberapa cara, yang paling sederhana adalah dengan melakukan Gaussian Elimination pada M ^T , menggantikan operasi penambahan baris dengan xor baris . Ini akan menghasilkan sejumlah vektor berbasis ruang nol, kombinasi apa pun yang akan menghasilkan solusi yang valid.

Konstruksi x cukup mudah. Ini hanyalah produk dari Ax + B untuk masing-masing y yang digunakan. Konstruksi d sedikit lebih rumit. Jika kita mengambil produk dari semua y , kita akan berakhir dengan nilai dengan 10s dari ribuan, jika tidak 100s dari ribuan digit, yang mana kita perlu menemukan akar kuadrat. Perhitungan ini tidak praktis mahal. Sebaliknya, kita bisa melacak kekuatan bahkan bilangan prima selama proses penyaringan, dan kemudian menggunakan dan dan XOR operasi pada vektor faktor non-square untuk merekonstruksi akar kuadrat.

Saya sepertinya telah mencapai batas 30000 karakter. Ahh, kurasa itu sudah cukup.

Nah, 38! +1 Anda memecahkan skrip php saya, tidak yakin mengapa. Bahkan, setiap semi-prime lebih dari 16 digit sudah lama mematahkan skrip saya.

Namun, dengan menggunakan 8980935344490257 (86028157 * 104395301) skrip saya mengatur waktu 25,963 detik pada komputer di rumah saya (AMD Phenom 9950 2,61GHz). Jauh lebih cepat daripada komputer kerja saya yang hampir 31 detik @ 2.93GHz Core 2 Duo.

php - 757 karakter termasuk baris baru

<?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;
}
?>

Saya akan tertarik untuk melihat algoritma yang sama ini di c atau bahasa lain yang dikompilasi.