快速素分解模型

Question

我正在寻找一个实现或清除算法，以便在python、伪代码或任何可读性很好的东西中获得N的主因子分解。有一些要求/事实：

• N在1到~20位之间
• 没有预先计算的查找表，但记忆功能很好。
• 不需要数学证明（例如，如果需要，可以依靠戈德巴赫猜想）
• 不需要精确，如果需要可以是概率的/确定性的

我需要一个快速的素数分解算法，不仅对它本身，而且对许多其他算法的使用，比如计算Eulerphi（n）。

我尝试过维基百科上的其他算法，但要么我不理解它们（ECM），要么我无法从算法中创建一个有效的实现（Pollard Brent）。

我对Pollard-Brent算法非常感兴趣，所以关于它的任何更多信息/实现都非常好。

谢谢！

编辑

在做了一些修改之后，我创建了一个非常快速的素数/因式分解模块。它结合了一个优化的尝试分割算法，波拉德布伦特算法，米勒拉宾素性测试和最快的素数筛我在网上找到。gcd是常规Euclid的gcd实现（二进制Euclid的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)

Answer 1

如果不想重新发明轮子，请使用库sympy

pip install sympy

使用函数^{}

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

Answer 2

Pollard Brent在Python中实现：

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

Answer 3

不需要使用primesbelow计算smallprimes，可以使用smallprimeset。

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

将primefactors分成两个函数来处理smallprimes和另一个函数来处理pollard_brent，这样可以节省一些迭代，因为smallprimes的所有幂都将从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的验证结果，可以减少使用Miller-Rabin素性检验的isprime中的重复次数。来自Wiki。

• 如果n<；1373653，则足以测试a=2和3
• 如果n<；9080191，则足以测试a=31和73
• 如果n<；4759123141，则足以测试a=2、7和61
• 如果n<；2152302898747，则足以测试a=2、3、5、7和11
• 如果n<；34749660383，则足以测试a=2、3、5、7、11和13
• 如果n<；341550071728321，则足以测试a=2、3、5、7、11、13和17。

编辑1：更正了if-else的返回调用，以将大因子附加到primefactors中的因子。