빠른 소인수 분해 모듈

빠른 소인수 분해 모듈

---

파이썬, 의사 코드 또는 기타 잘 읽을 수있는 N의 소인수 분해를 얻기위한 구현 또는 명확한 알고리즘 을 찾고 있습니다. 몇 가지 요구 / 사실이 있습니다.

• N 은 1 ~ 20 자리 사이입니다.
• 미리 계산 된 조회 테이블은 없지만 메모는 괜찮습니다.
• 수학적으로 증명할 필요가 없습니다 (예 : 필요한 경우 Goldbach 추측에 의존 할 수 있음).
• 정확할 필요는 없으며 필요한 경우 확률 적 / 결정적 일 수 있음

자체뿐만 아니라 Euler phi (n) 계산과 같은 다른 많은 알고리즘에서 사용하려면 빠른 소인수 분해 알고리즘이 필요합니다 .

Wikipedia 등에서 다른 알고리즘을 시도했지만 이해하지 못하거나 (ECM) 알고리즘에서 작동하는 구현을 만들 수 없습니다 (Pollard-Brent).

저는 Pollard-Brent 알고리즘에 정말 관심이 있으므로 더 많은 정보 / 구현이 정말 좋을 것입니다.

감사!

편집하다

조금 엉망으로 만든 후 꽤 빠른 프라임 / 팩 토라이 제이션 모듈을 만들었습니다. 최적화 된 시험 분할 알고리즘, Pollard-Brent 알고리즘, miller-rabin 소수성 테스트 및 인터넷에서 찾은 가장 빠른 소수 체를 결합합니다. gcd는 일반 유클리드의 GCD 구현입니다 (바이너리 유클리드의 GCD는 일반 유클리드의 GCD 보다 훨씬 느립니다).

하사품

오 기쁨, 현상금을 얻을 수 있습니다! 하지만 어떻게 이길 수 있습니까?

• 내 모듈에서 최적화 또는 버그를 찾으십시오.
• 대안 / 더 나은 알고리즘 / 구현을 제공합니다.

가장 완전하고 건설적인 대답이 현상금을받습니다.

마지막으로 모듈 자체 :

import random

def primesbelow(N):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    #""" Input N>=6, Returns a list of primes, 2 <= p < N """
    correction = N % 6 > 1
    N = {0:N, 1:N-1, 2:N+4, 3:N+3, 4:N+2, 5:N+1}[N%6]
    sieve = [True] * (N // 3)
    sieve[0] = False
    for i in range(int(N ** .5) // 3 + 1):
        if sieve[i]:
            k = (3 * i + 1) | 1
            sieve[k*k // 3::2*k] = [False] * ((N//6 - (k*k)//6 - 1)//k + 1)
            sieve[(k*k + 4*k - 2*k*(i%2)) // 3::2*k] = [False] * ((N // 6 - (k*k + 4*k - 2*k*(i%2))//6 - 1) // k + 1)
    return [2, 3] + [(3 * i + 1) | 1 for i in range(1, N//3 - correction) if sieve[i]]

smallprimeset = set(primesbelow(100000))
_smallprimeset = 100000
def isprime(n, precision=7):
    # http://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time
    if n < 1:
        raise ValueError("Out of bounds, first argument must be > 0")
    elif n <= 3:
        return n >= 2
    elif n % 2 == 0:
        return False
    elif n < _smallprimeset:
        return n in smallprimeset


    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1

    for repeat in range(precision):
        a = random.randrange(2, n - 2)
        x = pow(a, d, n)

        if x == 1 or x == n - 1: continue

        for r in range(s - 1):
            x = pow(x, 2, n)
            if x == 1: return False
            if x == n - 1: break
        else: return False

    return True

# https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
def pollard_brent(n):
    if n % 2 == 0: return 2
    if n % 3 == 0: return 3

    y, c, m = random.randint(1, n-1), random.randint(1, n-1), random.randint(1, n-1)
    g, r, q = 1, 1, 1
    while g == 1:
        x = y
        for i in range(r):
            y = (pow(y, 2, n) + c) % n

        k = 0
        while k < r and g==1:
            ys = y
            for i in range(min(m, r-k)):
                y = (pow(y, 2, n) + c) % n
                q = q * abs(x-y) % n
            g = gcd(q, n)
            k += m
        r *= 2
    if g == n:
        while True:
            ys = (pow(ys, 2, n) + c) % n
            g = gcd(abs(x - ys), n)
            if g > 1:
                break

    return g

smallprimes = primesbelow(1000) # might seem low, but 1000*1000 = 1000000, so this will fully factor every composite < 1000000
def primefactors(n, sort=False):
    factors = []

    for checker in smallprimes:
        while n % checker == 0:
            factors.append(checker)
            n //= checker
        if checker > n: break

    if n < 2: return factors

    while n > 1:
        if isprime(n):
            factors.append(n)
            break
        factor = pollard_brent(n) # trial division did not fully factor, switch to pollard-brent
        factors.extend(primefactors(factor)) # recurse to factor the not necessarily prime factor returned by pollard-brent
        n //= factor

    if sort: factors.sort()

    return factors

def factorization(n):
    factors = {}
    for p1 in primefactors(n):
        try:
            factors[p1] += 1
        except KeyError:
            factors[p1] = 1
    return factors

totients = {}
def totient(n):
    if n == 0: return 1

    try: return totients[n]
    except KeyError: pass

    tot = 1
    for p, exp in factorization(n).items():
        tot *= (p - 1)  *  p ** (exp - 1)

    totients[n] = tot
    return tot

def gcd(a, b):
    if a == b: return a
    while b > 0: a, b = b, a % b
    return a

def lcm(a, b):
    return abs((a // gcd(a, b)) * b)

---

Python으로 구현 된 Pollard-Brent :

https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/

---

바퀴를 재발 명하고 싶지 않다면 라이브러리 sympy를 사용하십시오.

pip install sympy

기능 사용 sympy.ntheory.factorint

>>> from sympy.ntheory import factorint
>>> factorint(10**20+1)
{73: 1, 5964848081: 1, 1676321: 1, 137: 1}

매우 큰 숫자를 고려할 수 있습니다.

>>> factorint(10**100+1)
{401: 1, 5964848081: 1, 1676321: 1, 1601: 1, 1201: 1, 137: 1, 73: 1, 129694419029057750551385771184564274499075700947656757821537291527196801: 1}

---

계산할 필요가 없다 smallprimes사용하여 primesbelow, 사용을 smallprimeset그것을 위해가.

smallprimes = (2,) + tuple(n for n in xrange(3,1000,2) if n in smallprimeset)

primefactors처리를 위해 두 개의 함수로 나누고를 위해 smallprimes다른 함수를 나누면 pollard_brentsmallprimes의 모든 거듭 제곱이 n에서 나뉘므로 몇 번의 반복을 절약 할 수 있습니다.

def primefactors(n, sort=False):
    factors = []

    limit = int(n ** .5) + 1
    for checker in smallprimes:
        print smallprimes[-1]
        if checker > limit: break
        while n % checker == 0:
            factors.append(checker)
            n //= checker


    if n < 2: return factors
    else : 
        factors.extend(bigfactors(n,sort))
        return factors

def bigfactors(n, sort = False):
    factors = []
    while n > 1:
        if isprime(n):
            factors.append(n)
            break
        factor = pollard_brent(n) 
        factors.extend(bigfactors(factor,sort)) # recurse to factor the not necessarily prime factor returned by pollard-brent
        n //= factor

    if sort: factors.sort()    
    return factors

Pomerance, Selfridge, Wagstaff 및 Jaeschke의 검증 된 결과를 고려 isprime하여 Miller-Rabin 소수성 검정을 사용 하는 반복을 줄일 수 있습니다 . 에서 위키 .

• n <1,373,653이면 a = 2와 3을 테스트하는 것으로 충분합니다.
• n <9,080,191이면 a = 31 및 73을 테스트하는 것으로 충분합니다.
• n <4,759,123,141이면 a = 2, 7, 및 61을 테스트하는 것으로 충분합니다.
• n <2,152,302,898,747이면 a = 2, 3, 5, 7, 및 11을 테스트하는 것으로 충분합니다.
• n <3,474,749,660,383이면 a = 2, 3, 5, 7, 11 및 13을 테스트하는 것으로 충분합니다.
• n <341,550,071,728,321이면 a = 2, 3, 5, 7, 11, 13 및 17을 테스트하는 것으로 충분합니다.

편집 1 :의 if-else요인에 큰 요인을 추가하기 위해 의 반환 호출을 수정 했습니다 primefactors.

---

현재에도 주목해야 할 점이 여러 개 있습니다.

1. Don't do checker*checker every loop, use s=ceil(sqrt(num)) and checher < s
2. checher should plus 2 each time, ignore all even numbers except 2
3. Use divmod instead of % and //

---

You should probably do some prime detection which you could look here, Fast algorithm for finding prime numbers?

You should read that entire blog though, there is a few algorithms that he lists for testing primality.

---

There's a python library with a collection of primality tests (including incorrect ones for what not to do). It's called pyprimes. Figured it's worth mentioning for posterity's purpose. I don't think it includes the algorithms you mentioned.

---

You could factorize up to a limit then use brent to get higher factors

from fractions import gcd
from random import randint

def brent(N):
   if N%2==0: return 2
   y,c,m = randint(1, N-1),randint(1, N-1),randint(1, N-1)
   g,r,q = 1,1,1
   while g==1:             
       x = y
       for i in range(r):
          y = ((y*y)%N+c)%N
       k = 0
       while (k<r and g==1):
          ys = y
          for i in range(min(m,r-k)):
             y = ((y*y)%N+c)%N
             q = q*(abs(x-y))%N
          g = gcd(q,N)
          k = k + m
       r = r*2
   if g==N:
       while True:
          ys = ((ys*ys)%N+c)%N
          g = gcd(abs(x-ys),N)
          if g>1:  break
   return g

def factorize(n1):
    if n1==0: return []
    if n1==1: return [1]
    n=n1
    b=[]
    p=0
    mx=1000000
    while n % 2 ==0 : b.append(2);n//=2
    while n % 3 ==0 : b.append(3);n//=3
    i=5
    inc=2
    while i <=mx:
       while n % i ==0 : b.append(i); n//=i
       i+=inc
       inc=6-inc
    while n>mx:
      p1=n
      while p1!=p:
          p=p1
          p1=brent(p)
      b.append(p1);n//=p1 
    if n!=1:b.append(n)   
    return sorted(b)

from functools import reduce
#n= 2**1427 * 31 #
n= 67898771  * 492574361 * 10000223 *305175781* 722222227*880949 *908909
li=factorize(n)
print (li)
print (n - reduce(lambda x,y :x*y ,li))

---

I just ran into a bug in this code when factoring the number 2**1427 * 31.

  File "buckets.py", line 48, in prettyprime
    factors = primefactors.primefactors(n, sort=True)
  File "/private/tmp/primefactors.py", line 83, in primefactors
    limit = int(n ** .5) + 1
OverflowError: long int too large to convert to float

This code snippet:

limit = int(n ** .5) + 1
for checker in smallprimes:
    if checker > limit: break
    while n % checker == 0:
        factors.append(checker)
        n //= checker
        limit = int(n ** .5) + 1
        if checker > limit: break

should be rewritten into

for checker in smallprimes:
    while n % checker == 0:
        factors.append(checker)
        n //= checker
    if checker > n: break

which will likely perform faster on realistic inputs anyway. Square root is slow — basically the equivalent of many multiplications —, smallprimes only has a few dozen members, and this way we remove the computation of n ** .5 from the tight inner loop, which is certainly helpful when factoring numbers like 2**1427. There's simply no reason to compute sqrt(2**1427), sqrt(2**1426), sqrt(2**1425), etc. etc., when all we care about is "does the [square of the] checker exceed n".

The rewritten code doesn't throw exceptions when presented with big numbers, and is roughly twice as fast according to timeit (on sample inputs 2 and 2**718 * 31).

Also notice that isprime(2) returns the wrong result, but this is okay as long as we don't rely on it. IMHO you should rewrite the intro of that function as

if n <= 3:
    return n >= 2
...

참고URL : https://stackoverflow.com/questions/4643647/fast-prime-factorization-module