python 最快 因式分解_python – 快速素因子分解模块

我正在寻找一个实现或清除算法获得N的素数因子分解在Python，伪码或任何其他可读性。有几个要求/事实：

> N在1到〜20之间

>没有预先计算的查找表，记忆是好的。

>不需要在数学上证明(例如，如果需要，可以依靠Goldbach猜想)

>不需要精确，如果需要，允许是概率/确定性的

我需要一个快速素因子分解算法，不仅对于自身，而是用于许多其他算法，如计算欧拉phi(n)。

我试过维基百科等其他算法，但我不能理解他们(ECM)或我不能创建一个工作的实现从算法(波拉德布伦特)。

我对Pollard-Brent算法真的很感兴趣，所以任何更多的信息/实现将是非常好的。

谢谢！

编辑

在搞乱了一点后，我创建了一个相当快的素数/分解模块。它结合了优化的试验划分算法，Pollard-Brent算法，米勒 – 拉宾素性测试和在互联网上找到的最快的primesieve。 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)