Najszybsze rozkładanie na czynniki pierwsze

Napisz program, który w jak najkrótszym czasie będzie rozkładał liczbę na półpierwszą.

Do celów testowych użyj tego: 38! +1 (523022617466601111760007224100074291200000001)

Jest równa: 14029308060317546154181 × 37280713718589679646221

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

Korzystanie z wielokrotnego wielomianowego sita kwadratowego . Uznałem to za wyzwanie dla kodu, więc zdecydowałem się nie używać żadnych bibliotek zewnętrznych (jak logsądzę, poza funkcją standardową ). Podczas pomiaru czasu należy użyć PyPy JIT , ponieważ powoduje on taktowanie 4-5 razy szybsze niż cPython .

Aktualizacja (2013-07-29):
Od czasu pierwszego opublikowania wprowadziłem kilka drobnych, ale znaczących zmian, które zwiększają ogólną prędkość około 2,5x.

Aktualizacja (2014-08-27):
Ponieważ ten post wciąż jest przedmiotem zainteresowania, zaktualizowałem my_math.pypoprawianie dwóch błędów dla każdego, kto może go używać:

• isqrtbył wadliwy, czasami wytwarzając niepoprawny wynik dla wartości bardzo zbliżonych do idealnego kwadratu. Zostało to poprawione, a wydajność wzrosła dzięki zastosowaniu znacznie lepszego seedu.
• is_primezostał zaktualizowany. Moja poprzednia próba usunięcia idealnego kwadratowego 2-sprpsa była w najlepszym wypadku bez serca. Dodałem kontrolę 3-sprp - technikę stosowaną przez Mathmatica - aby upewnić się, że testowana wartość jest wolna od kwadratów.

Aktualizacja (24.11.2014):
Jeśli pod koniec obliczeń nie znaleziono nietrywialnych kongruencji, program wyświetla teraz dodatkowe wielomiany. Zostało to wcześniej oznaczone w kodzie jako 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)

moja_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

Przykładowe 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

Uwaga: niestosowanie tej --verboseopcji da nieco lepsze czasy:

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

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

Podstawowe koncepcje

Zasadniczo sito kwadratowe opiera się na następującej obserwacji: każdy nieparzysty kompozyt n może być reprezentowany jako:

Nie jest to trudne do potwierdzenia. Ponieważ n jest nieparzyste, odległość między dowolnymi dwoma kofaktorami n musi wynosić nawet 2d , gdzie x jest punktem środkowym między nimi. Co więcej, ta sama relacja obowiązuje dla dowolnej wielokrotności n

Zauważ, że jeśli takie x i d można znaleźć, natychmiast spowoduje to (niekoniecznie pierwszą) wartość n , ponieważ x + d i x - d dzielą n z definicji. Relację tę można dodatkowo osłabić - w konsekwencji dopuszczenia potencjalnych trywialnych zgodności - do następującej postaci:

Ogólnie więc, jeśli znajdziemy dwa idealne kwadraty, które są równoważne mod n , to jest całkiem prawdopodobne, że możemy bezpośrednio wytworzyć współczynnik n a la gcd (x ± d, n) . Wydaje się to dość proste, prawda?

Tyle że nie. Jeśli zamierzamy przeprowadzić wyczerpujące przeszukanie wszystkich możliwych x , musielibyśmy przeszukać cały zakres od [ √ n , √ ( 2n ) ], który jest nieznacznie mniejszy niż pełny podział próbny, ale wymaga również kosztownej is_squareoperacji przy każdej iteracji do potwierdź wartość d . O ile nie jest znany wcześniej, że n ma czynniki bardzo blisko √ n , podział próbny może być szybciej.

Być może możemy jeszcze bardziej osłabić tę relację. Załóżmy, że wybraliśmy x , na przykład dla

pełna faktoryzacja liczb pierwszych y jest łatwo znana. Gdybyśmy mieli wystarczającą liczbę takich relacji, powinniśmy być w stanie skonstruować odpowiednie d , jeśli wybieramy liczbę y taką, że ich iloczyn jest idealnym kwadratem; to znaczy wszystkie czynniki pierwsze są wykorzystywane parzystą liczbę razy. W rzeczywistości, jeśli mamy więcej takich y niż całkowita liczba unikalnych czynników pierwszych, które one zawierają, gwarantuje się, że istnieje rozwiązanie; Staje się układem równań liniowych. Powstaje teraz pytanie, jak wybraliśmy takie x ? Tu właśnie zaczyna się przesiewanie.

Sito

Rozważ wielomian:

Następnie dla dowolnej liczby pierwszej p i liczby całkowitej k obowiązuje następująca zasada:

Oznacza to, że po rozwiązaniu pierwiastków wielomianu mod p - to znaczy, że znalazłeś x taki, że y (x) ≡ 0 (mod p) , ergo y jest podzielne przez p - to znalazłeś nieskończoną liczbę takich x . W ten sposób możesz przesiać w zakresie x , identyfikując małe czynniki pierwsze y , miejmy nadzieję, znajdując takie, dla których wszystkie czynniki pierwsze są małe. Takie liczby, znane jako k-gładkie , gdzie k jest największym zastosowanym współczynnikiem podstawowym.

Jednak z tym podejściem wiąże się kilka problemów. Nie wszystkie wartości x są odpowiednie, w rzeczywistości jest ich bardzo niewiele, skupionych wokół √ n . Mniejsze wartości staną się w dużej mierze ujemne (z powodu terminu -n ), a większe wartości staną się zbyt duże, tak że jest mało prawdopodobne, aby ich pierwsza faktoryzacja składała się tylko z małych liczb pierwszych. Będzie wiele takich x , ale chyba że faktoryzowany kompozyt jest bardzo mały, jest bardzo mało prawdopodobne, że znajdziesz wystarczająco dużo wygładzeń, aby uzyskać faktoryzację. I tak dla większego n konieczne staje się przesianie wielu wielomianów danej formy.

Wiele wielomianów

Więc potrzebujemy więcej wielomianów do przesiewania? Co powiesz na to:

To zadziała. Zauważ, że A i B mogą być dosłownie dowolnymi liczbami całkowitymi, a matematyka nadal obowiązuje. Wszystko, co musimy zrobić, to wybrać kilka losowych wartości, rozwiązać pierwiastek wielomianu i przesiać wartości bliskie zeru. W tym momencie moglibyśmy nazwać to wystarczająco dobrym: jeśli rzucisz wystarczającą liczbę kamieni w losowych kierunkach, prędzej czy później zepsujesz okno.

Tyle że z tym też jest problem. Jeśli nachylenie wielomianu jest duże na punkcie przecięcia x, co będzie, jeśli nie będzie względnie płaskie, będzie tylko kilka odpowiednich wartości do przesiewania na wielomian. To zadziała, ale skończysz przesiewać wiele wielomianów, zanim dostaniesz to, czego potrzebujesz. Czy możemy zrobić lepiej?

Możemy zrobić lepiej. Obserwacja, w wyniku Montgomery, jest następująca: jeśli A i B są wybrane w taki sposób, że istnieje pewne C spełniające

wtedy cały wielomian może zostać przepisany jako

Ponadto, jeśli A zostanie wybrany jako idealny kwadrat, wiodący składnik A może zostać pominięty podczas przesiewania, co skutkuje znacznie mniejszymi wartościami i znacznie bardziej płaską krzywą. W przypadku takiego rozwiązania istnieje, n musi być kwadratowe pozostałość mod √ , który może być znany natychmiast obliczania symbol Legendre'a :
( N | √A ) = 1 . Należy zauważyć, że aby rozwiązać B , należy znać całkowite rozkładanie na czynniki pierwsze √A (aby wziąć modułowy pierwiastek kwadratowy √n (mod √A) ), dlatego właśnie √A jest zwykle wybierana jako liczba pierwsza.

Można wówczas wykazać, że jeśli , to dla wszystkich wartości x ∈ [ -M, M ] :

A teraz wreszcie mamy wszystkie elementy niezbędne do wdrożenia naszego sita. A może my?

Potęgi liczb pierwszych jako czynniki

Nasze sito, jak opisano powyżej, ma jedną poważną wadę. Można go określić, które wartości x spowoduje y podzielna przez p , ale nie można określić, czy to y jest podzielna przez siły z p . Aby to ustalić, musielibyśmy przeprowadzić podział próbny na przesiewaną wartość, dopóki nie będzie już podzielna przez p . Wydawało się, że osiągnęliśmy impassé: cały punkt sita był taki, że nie musieliśmy tego robić. Czas sprawdzić instrukcję.

To wygląda całkiem przydatne. Jeśli suma ln wszystkich małych czynników pierwszych y jest zbliżona do oczekiwanej wartości ln (y) , to prawie pewne, że y nie ma innych czynników. Ponadto, jeśli obniżymy nieco wartość oczekiwaną, możemy również zidentyfikować wartości jako gładkie, które mają kilka mocy liczb pierwszych jako czynników. W ten sposób możemy wykorzystać sito jako proces „wstępnej kontroli” i uwzględnić tylko te wartości, które prawdopodobnie będą gładkie.

Ma to również kilka innych zalet. Należy pamiętać, że małe liczby pierwsze przyczyni się bardzo niewiele do ln sumy, ale mimo to wymaga najwięcej czasu sita. Przesiewanie wartość 3 wymaga więcej czasu niż 11, 13, 17, 19 i 23 w połączeniu . Zamiast tego możemy po prostu pominąć kilka pierwszych liczb pierwszych i odpowiednio obniżyć próg, zakładając, że pewien procent z nich minąłby.

Innym rezultatem jest to, że pewna liczba wartości będzie mogła „prześlizgnąć się”, które są w większości gładkie, ale zawierają jeden duży kofaktor. Możemy po prostu odrzucić te wartości, ale załóżmy, że znaleźliśmy inną, w większości gładką wartość, z dokładnie tym samym kofaktorem. Następnie możemy użyć tych dwóch wartości do skonstruowania użytecznego y ; ponieważ ich produkt będzie zawierał ten duży kofaktor do kwadratu, nie trzeba go już brać pod uwagę.

Kładąc wszystko razem

Ostatnią rzeczą, jaką musimy zrobić, to użyć tych wartości y skonstruować odpowiednią X i d . Załóżmy, że bierzemy pod uwagę tylko kwadratowe współczynniki y , to znaczy czynniki pierwsze mocy nieparzystej. Następnie każde y można wyrazić w następujący sposób:

które można wyrazić w postaci macierzy:

Problemem staje się znalezienie wektora v takiego, że vM = ⦳ (mod 2) , gdzie ⦳ jest wektorem zerowym. Oznacza to, że w celu rozwiązania dla lewej przestrzeni NULL M . Można tego dokonać na wiele sposobów, z których najprostszy jest wykonać eliminacji Gaussa o M ^T , zastępując operacji dodawania rząd rzędem XOR . Spowoduje to powstanie szeregu wektorów bazowych o zerowej przestrzeni, których dowolna kombinacja da prawidłowe rozwiązanie.

Konstrukcja x jest dość prosta. Jest to po prostu iloczyn Ax + B dla każdego zastosowanego y . Konstrukcja d jest nieco bardziej skomplikowana. Gdybyśmy wzięli iloczyn całego y , otrzymalibyśmy wartość z 10s tysięcy, jeśli nie 100s tysięcy cyfr, dla których musimy znaleźć pierwiastek kwadratowy. To obliczenie jest niepraktycznie drogie. Zamiast tego, możemy śledzić nawet potęg liczb pierwszych podczas procesu przesiewania, a następnie użyj a i xor operacje na wektorach czynników nie-kwadratowych zrekonstruować pierwiastek kwadratowy.

Wydaje mi się, że osiągnąłem limit 30000 znaków. Achh cóż, przypuszczam, że to wystarczy.

Cóż, twój 38! +1 złamał mój skrypt php, nie jestem pewien, dlaczego. W rzeczywistości każda półpierwsza liczba składająca się z ponad 16 cyfr łamie mój skrypt.

Jednak używając 8980935344490257 (86028157 * 104395301) mój skrypt zarządzał czasem 25,963 sekundy na moim komputerze domowym (2,61 GHz AMD Phenom 9950). Dużo szybszy niż mój komputer roboczy, który miał prawie 31 sekund przy 2,93 GHz Core 2 Duo.

php - 757 znaków w tym. nowe linie

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

Byłbym zainteresowany tym samym algorytmem w c lub innym skompilowanym języku.