我正在寻找一种实现或清晰的算法,以便在 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://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)
最佳答案
如果您不想重新发明轮子,请使用库 sympy
pip install sympy
Given a positive integer
n
,factorint(n)
returns a dict containing the prime factors ofn
as keys and their respective multiplicities as values. For example:
例子:
>>> 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}
关于python - 快速素数分解模块,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/4643647/