Python (PyPy JIT v1.9 포함) ~ 1.9 초

사용하여 여러 다항식 이차 체를 . 나는 이것을 코드 도전으로 삼았으므로 외부 라이브러리 (표준 log함수 이외의 다른)를 사용하지 않기로 결정했다 . 타이밍을 설정하는 경우 PyPy JIT를 사용해야합니다 . cPython 보다 4-5 배 빠릅니다 .

업데이트 (2013-07-29) :
처음 게시 한 이후로 몇 가지 사소하지만 크게 변경하여 전체 속도를 약 2.5 배 증가시킵니다.

업데이트 (2014-08-27) :
이 게시물이 여전히 주목을 받고 my_math.py있으므로, 사용중인 모든 사람을 위해 두 가지 오류 수정을 업데이트 했습니다.

• isqrt결함이 있으며 때로는 완벽한 제곱에 매우 가까운 값에 대해 잘못된 출력을 생성합니다. 이것은 수정되었으며 훨씬 더 나은 시드를 사용하여 성능이 향상되었습니다.
• is_prime업데이트되었습니다. 완벽한 사각 2 스프링을 제거하려는 이전의 시도는 반쯤 마음에 들었습니다. 테스트 된 값에 제곱이 없는지 확인하기 위해 Mathmatica에서 사용하는 기술인 3-sprp check를 추가했습니다.

업데이트 (2014-11-24) :
계산이 끝날 때 사소한 합치점이 발견되지 않으면 프로그램은 이제 추가 다항식을 체질합니다. 이것은 이전에 코드에서로 표시되었습니다 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

샘플 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

참고 : --verbose옵션을 사용하지 않으면 타이밍이 약간 향상됩니다.

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

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

기본 개념

일반적으로, 2 차 체는 다음 관찰에 기초한다 : 임의의 홀수 복합 n 은 다음과 같이 표현 될 수있다 :

확인하기가 그리 어렵지 않습니다. 이후 , n은 홀수의 두 보조 인자 사이의 거리 , n은 짝수이어야 2D , X는 그들 사이의 중간 지점이다. 또한, n의 배수에 대해 동일한 관계가 유지됩니다.

그러한 x 와 d 를 찾을 수 있다면 , x + d 와 x-d는 모두 n 을 정의로 나누기 때문에 ( n 은 반드시 소수는 아님) n의 결과를 가져옵니다 . 사소한 합의를 허용 한 결과 다음과 같은 형태로이 관계가 더욱 약화 될 수 있습니다.

따라서 일반적으로 mod n 과 동일한 두 개의 완벽한 제곱을 찾을 수 있다면 n a la gcd (x ± d, n) 의 계수를 직접 생성 할 가능성이 높습니다 . 꽤 간단 해 보이죠?

그렇지 않으면 빼고 가능한 모든 x에 대해 철저한 검색을 수행하려면 [ √ n , √ ( 2n ) ] 에서 전체 범위를 검색해야합니다 . 이는 전체 시험 분할보다 약간 작지만 is_square각 반복에 대해 비싼 작업이 필요합니다. d 값을 확인하십시오 . n 이 √ n에 매우 가까운 인자를 가지고 있다는 것을 미리 알지 못한다면 , 시험 분할이 더 빠를 것입니다.

아마도 우리는이 관계를 더욱 약화시킬 수 있습니다. 우리가 x를 선택 했다고 가정 해 봅시다.

y 의 완전 소인수 분해 는 쉽게 알려져 있습니다. 만약 우리가 그러한 관계가 충분하다면, 우리가 그들의 제품이 완벽한 제곱이되도록 많은 y를 선택한다면 우리 는 적절한 d 를 구성 할 수 있어야한다 . 즉, 모든 주요 요소가 짝수 번 사용됩니다. 실제로, 우리가 포함하고있는 고유 한 주요 요소의 총 수보다 y 가 더 많으면 해결책이 존재합니다. 그것은 선형 방정식 시스템이됩니다. 이제 우리는 어떻게 그런 x를 선택 했는가? 그것은 체질이 시작되는 곳입니다.

체

다항식을 고려하십시오.

그런 다음 소수 p 및 정수 k 에 대해 다음이 적용됩니다.

이것은 다항식 mod p 의 근을 풀고 나면 , 즉 y (x) (0 (mod p) 와 같은 x를 찾았고 , ergo y 는 p 로 나눌 수 있다는 것을 의미합니다. 그런 x의 . 이런 식으로, x 의 범위에서 체를 걸러서 y의 작은 소인수를 식별 하고 모든 소인수가 작은 것을 찾을 수 있기를 바랍니다. k-smooth 라고하는 숫자로 , 여기서 k 는 가장 큰 소수입니다.

그러나이 방법에는 몇 가지 문제가 있습니다. x의 모든 값 이 적절 하지는 않지만 , 실제로는 √ n 주위에 거의 없습니다 . 더 작은 값은 ( -n 항 으로 인해) 크게 음수가되고 , 더 큰 값은 너무 커져서, 소인수 분해가 작은 소수만으로 구성 될 가능성은 없습니다. 이러한 x 는 여러 가지가 있지만 팩토링하는 컴포지트가 매우 작은 경우가 아니면 인수 분해를하기에 충분한 스무딩을 찾을 가능성이 거의 없습니다. 따라서 더 큰 n 의 경우 주어진 형태의 여러 다항식 을 체로 치러야 합니다.

여러 다항식

그래서 우리는 더 많은 다항식이 필요합니까? 이건 어때요:

작동합니다. 참고 것을 와 B는 말 그대로 임의의 정수 값이 될 수 있고, 수학은 여전히 보유하고 있습니다. 우리가해야 할 일은 몇 가지 임의의 값을 선택하고 다항식의 근을 풀고 0에 가까운 값을 체킹하는 것입니다. 이 시점에서 우리는 그것을 충분히 잘 부를 수 있습니다. 만약 당신이 충분한 방향으로 돌을 던지면 조만간 창문을 깰 수밖에 없습니다.

그 외에는 문제가 있습니다. 다항식의 기울기가 x 절편에서 크면 (절대 평평하지 않은 경우) 다항식 당 체에 적합한 몇 가지 값만있을 것입니다. 작동하지만 필요한 것을 얻기 전에 다항식을 많이 체질하게됩니다. 더 잘할 수 있을까요?

우리는 더 잘할 수 있습니다. 몽고메리 의 결과는 다음과 같습니다. A 와 B 를 선택 하면 만족스러운 C 가 존재합니다.

다항식 전체를 다음과 같이 다시 쓸 수 있습니다.

또한 A를 완벽한 정사각형으로 선택하면 체질하는 동안 선행 A 항을 무시할 수있어 값이 더 작고 곡선이 훨씬 평평 해집니다. 그러한 해가 존재하기 위해서는 n 은 2 차 잔차 mod √ A 이어야하며 , Legendre 기호 를 계산하여 즉시 알 수 있습니다 :
( n | √A ) = 1 . 에 대한 해결하기 위해 참고 B 의 전체 소인수 분해 √A은 (모듈 형 제곱근 취하기 위해 알려질 필요가 √n (모드 √A) 이유) √A는 일반적으로 소수로 선택됩니다.

그러면이면 x ∈ [ -M, M ] 의 모든 값에 대해 다음을 표시 할 수 있습니다 .

그리고 마지막으로, 체를 구현하는 데 필요한 모든 구성 요소가 있습니다. 아니면 우리?

요인으로서의 소수

위에서 설명한 것처럼 우리의 체에는 하나의 큰 결함이 있습니다. 이 값을 식별 할 수 X는 (A)에 발생한다 Y 로 나누어 P 하지만이 있는지의 여부를 식별 할 수있는 Y는 바이 나누어 전력 의 P . 이를 결정하기 위해 더 이상 p 로 나눌 수 없을 때까지 체질 할 값에 대해 시험 분할을 수행해야합니다 . 우리는 임 페레스에 도달 한 것 같았습니다. 체의 요점은 그렇게 할 필요가 없었습니다 . 플레이 북을 확인할 시간입니다.

꽤 유용한 것 같습니다. y 의 모든 작은 소인수 의 ln 의 합이 ln (y) 의 예상 값에 가까워 지면 y 에 다른 요인이없는 것으로 간주 됩니다. 또한 예상 값을 약간 낮추면 여러 소수의 소수를 요인으로하는 값을 부드럽게 식별 할 수도 있습니다. 이러한 방식으로, 체를 ‘사전 스크리닝’프로세스로 사용할 수 있으며 매끄럽게 될 수있는 값만 고려할 수 있습니다.

이것은 몇 가지 다른 장점도 있습니다. 작은 소수는 ln 합계에 거의 영향을 미치지 않지만 가장 많은 시간이 필요합니다. 값 3을 체질하려면 11, 13, 17, 19 및 23을 합한 것보다 더 많은 시간이 필요합니다 . 대신, 우리는 처음 몇 소수를 건너 뛰고 특정 백분율이 통과했다고 가정하여 임계 값을 적절하게 조정할 수 있습니다.

또 다른 결과는, 많은 값들이 ‘slip through’될 수 있다는 것인데, 이것은 대부분 매끄럽지 만 하나의 큰 보조 인자를 포함합니다. 우리는이 값들을 버릴 수는 있지만, 정확히 동일한 보조 인자를 가진 다른 부드러운 값을 발견했다고 가정합니다. 그런 다음이 두 값을 사용하여 사용 가능한 y 를 구성 할 수 있습니다 . 그들의 제품에는이 큰 보조 인자 제곱이 포함될 것이기 때문에 더 이상 고려할 필요가 없습니다.

함께 모아서

마지막으로해야 할 일은 y 값을 사용 하여 적절한 x 와 d를 만드는 것 입니다. y 의 비 제곱 요인 , 즉 홀수 ​​거듭 제곱의 주요 요인 만 고려한다고 가정 합니다. 그런 다음 각 y 는 다음과 같은 방식으로 표현 될 수 있습니다.

행렬 형식으로 표현할 수 있습니다.

그러면 문제 는 vM = ⦳ (mod 2) 와 같은 벡터 v 를 찾는데 , 여기서 ⦳ 는 널 벡터입니다. 즉, M 의 왼쪽 널 공간을 해결합니다 . 이것은 여러 가지 방법으로 수행 될 수 있으며, 가장 간단한 방법 은 행 추가 연산을 행 xor로 대체하여 M ^T 에서 가우시안 제거를 수행하는 것입니다 . 이로 인해 다수의 널 공간 기반 벡터가 생성되며, 이들의 조합은 유효한 솔루션을 생성합니다.

x 의 구성 은 매우 간단합니다. 사용 된 각 y에 대한 Ax + B 의 곱입니다 . d 의 구성 은 약간 더 복잡합니다. 모든 y 의 곱을 취한다면, 수십만 자릿수가 아닌 10 만 자릿수의 값으로 끝나고, 제곱근을 찾아야합니다. 이 계산은 비현실적으로 비쌉니다. 대신, 우리는 체질 과정에서 짝수의 소수를 추적 한 다음, 제곱근을 재구성하기 위해 비 제곱 요인의 벡터에 대해 and 및 xor 연산 을 사용할 수 있습니다 .

30000 자 한도에 도달 한 것 같습니다. 글쎄요, 충분하다고 생각합니다.

답변

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