快速素数分解模块
Posted
技术标签:
【中文标题】快速素数分解模块【英文标题】：Fast prime factorization module 【发布时间】：2011-01-10 04:19:07 【问题描述】：我正在寻找一种实现或清晰的算法，以便在 python、伪代码或其他任何东西中获得 N 的主要因素 -可读。有一些要求/限制：
N 介于 1 到 ~20 位之间 没有预先计算的查找表，但是记忆很好 不需要经过数学证明（例如，如果需要，可以依赖哥德巴赫猜想） 不需要精确，如果需要，可以是概率/确定性我需要一个快速的素数分解算法，不仅是为了它本身，而且是为了在许多其他算法中使用，比如计算欧拉 phi(n)。
我尝试过来自 Wikipedia 等的其他算法，但要么我无法理解它们 (ECM)，要么我无法从该算法 (Pollard-Brent) 创建一个有效的实现。
我对 Pollard-Brent 算法真的很感兴趣，所以如果有更多关于它的信息/实现会非常好。
谢谢！
编辑
经过一番折腾后，我创建了一个非常快速的素数/因式分解模块。它结合了优化的试除法算法、Pollard-Brent 算法、米勒-拉宾素数测试和我在互联网上找到的最快的素筛。 gcd 是常规 Euclid 的 GCD 实现（二进制 Euclid 的 GCD 比常规的慢得多）。
赏金
哦，快乐，可以获得赏金！但是我怎样才能赢呢？
在我的模块中查找优化或错误。 提供替代/更好的算法/实现。最完整/最有建设性的答案将获得赏金。
最后是模块本身：
import random
def primesbelow(N):
# http://***.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)
【问题讨论】：
@wheaties - 这就是while checker*checker <= num 的用途。
你可能会发现这个帖子很有用：***.com/questions/1024640/calculating-phik-for-1kn/…
为什么标准库中没有这样的东西？当我搜索时，我只找到一百万个 Project Euler 解决方案提案，以及其他人指出其中的缺陷。这不就是库和错误报告的用途吗？
@endolith 除了像 Prject Euler 这样的东西之外，它并没有太多用处。当然不足以将其放入标准库中。
@nightcracker：分解数字没有实际用途？？
【参考方案1】：
如果您不想重新发明***，请使用库sympy
pip install sympy
使用函数sympy.ntheory.factorint
给定一个正整数
n，factorint(n)返回一个包含n的主要因子作为键及其各自的多重性 作为价值观。例如：
例子：
>>> 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
【讨论】：
感谢分享这个，@Colonel Panic。这正是我一直在寻找的：一个维护良好的库中的整数分解模块，而不是代码的 sn-ps。 是的，sympy 很简单！【参考方案2】：无需使用primesbelow 计算smallprimes，使用smallprimeset 即可。
smallprimes = (2,) + tuple(n for n in xrange(3,1000,2) if n in smallprimeset)
将你的primefactors 分成两个函数来处理smallprimes 和另一个用于处理pollard_brent，这样可以节省几次迭代，因为小素数的所有幂都将从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。
编辑 1：更正了 if-else 的返回调用以将大因素附加到 primefactors 中的因素。
【讨论】：
享受你的+100（自赏金以来，你是唯一一个回答的人）。不过，您的bigfactors 效率非常低，因为factors.extend(bigfactors(factor)) 递归回大因子，这完全是错误的（如果 pollard-brent 找到因子 234892，您不想再用 pollard-brent 分解它）。如果您将factors.extend(bigfactors(factor)) 更改为factors.extend(primefactors(factor, sort)) 就可以了。
一个素数因子称为大因子，那么很明显，pollard-brent 获得的下一个因子中没有小素数的幂。
如果效率低下，我不会回答这个问题。一旦调用从primefactors转到bigfactors，n中将没有小于1000的因子，因此pollard-brent无法返回因子小于1000的数字。如果不清楚，请回复以便我将用更多解释编辑我的答案
该死的 NVM，当然。如果 N 不包含小素数找到的任何因子，则 N 的因子 F 也不会 >.
您的primefactors 确实效率低下，因为bigfactors 被称为所有具有小因数且最大因数大于2 的数字。这是因为最后一个因数因n > int(n ** .5) + 1 而被跳过。在通过if n in smallprimes:`factors.append(n)`else:`...bigfactors.... The second bug is that return any_list.extend(...)` 调用 bigfactorc 之前可以很容易地修复它返回 None。它也很容易修复。【参考方案3】：
即使是现在的，也有几个地方需要注意。
-
不要在每个循环中都使用
checker*checker，使用s=ceil(sqrt(num)) 和checher < s
checher 每次都要加 2，忽略除 2 以外的所有偶数
使用divmod 代替% 和//
【讨论】：
我需要做 checker*checker 因为 num 不断减少。我将实现偶数跳过。 divmod 大大减少了函数（它将在每个循环中计算 //，而不是仅在检查器除以 n 时）。 @night，你只需要在你改变num时重新计算s然后
是的，我认为在搞砸的同时 :) 重新计算 sqrt 似乎比 checker*checker 更快。
@nightcracker: 让N=n*n+1、ceil(sqrt(N)) 的成本是n*n 的2~4 倍，num 不会那么频繁地变化。
乘法比平方根快，Python 无论如何也很难取 bignums 的平方根。当前版本的代码出现此错误：limit = int(n ** .5) + 1 ... OverflowError: long int too large to convert to float【参考方案4】：
你可能应该做一些你可以看这里的素数检测， Fast algorithm for finding prime numbers?
不过，您应该阅读整个博客，他列出了一些用于测试素性的算法。
【讨论】：
【参考方案5】：有一个 Python 库，其中包含一系列素性测试（包括不正确的测试）。它被称为pyprimes。认为为了后代的目的值得一提。我不认为它包括你提到的算法。
【讨论】：
【参考方案6】：你可以分解到一个极限，然后使用布伦特得到更高的因子
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))
【讨论】：
【参考方案7】：我在分解数字 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
这段代码sn-p：
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
应该改写成
for checker in smallprimes:
while n % checker == 0:
factors.append(checker)
n //= checker
if checker > n: break
无论如何，这可能会在实际输入上执行得更快。平方根很慢——基本上相当于许多乘法——smallprimes 只有几十个成员，这样我们就可以从紧密的内部循环中删除n ** .5 的计算，这在分解像@987654327 这样的数字时肯定很有帮助@。根本没有理由计算 sqrt(2**1427)、sqrt(2**1426)、sqrt(2**1425) 等等，而我们只关心“检查器的 [square of the] 是否超过 n”。
根据timeit（样本输入2 和2**718 * 31），重写后的代码在呈现大数字时不会引发异常，和大约是两倍。 p>
还要注意isprime(2) 返回了错误的结果，但是只要我们不依赖它就可以了。恕我直言，您应该将该功能的介绍重写为
if n <= 3:
return n >= 2
...
【讨论】：
以上是关于快速素数分解模块的主要内容，如果未能解决你的问题，请参考以下文章