快速素因子分解模块 [英] Fast prime factorization module
问题描述
我要寻找一个执行或清除算法获取的的 N 的首要因子在任一蟒蛇，伪code或别的也可读。有几个要求/事实：
- 的 N 的为1至〜20位数字
- 没有pre计算的查找表，记忆化是好的，但。
- 不需要用数学证明（如可以依靠，如果需要的哥德巴赫猜想）
- 不需要是precise，允许是概率性/确定性，如果需要
我需要一个快速的质因子分解算法，不仅为自己，而是使用在许多其他类似的算法计算欧拉的岛（N）的。
我试图从维基百科和其他算法，但要么我不明白他们（ECM），或者我不能创建的算法（波拉德 - 布伦特）工作实施
我在波拉德 - 布伦特算法很感兴趣，所以任何详细信息/它的实现将是非常好的。
谢谢！
修改
附近有一座小我创建了一个pretty的快速黄金/分解模块搞乱之后。它结合了优化的审判庭算法，波拉德 - 布伦特算法，米勒 - 拉宾检验并以最快的primesieve我发现在互联网上。 GCD是一个普通的欧几里得的GCD实现（二进制欧几里德的GCD是多慢然后定期的）。
赏金
哦喜悦，赏金可以获取！但是，我怎么能赢得它？
- 找到我的模块中的optimalization或缺陷。
- 在提供替代/更好的算法/实现。
这是最完整的/建设性的答案得到赏金。
和最后的模块本身：
导入随机
高清primesbelow（N）：
＃ http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
＃输入N＆GT; = 6，返回的素数，第2版的列表; = P＆其中; N
修正= N％6' 1
N = {0：N，1：N-1，2：N + 4,3：N + 3,4：N + 2,5：N + 1} [N％6]
筛= [真] *（N // 3）
筛[0] = FALSE
因为我在范围内（INT（N ** .5）// 3 + 1）：
如果筛[我]：
K =（3 * I + 1）| 1
筛[K * K // 3 :: 2 * K] = [假] *（（N // 6 - （K * k）的// 6 - 1）// K + 1）
筛[（K * K + 4 * K - 2 * K *（我％2））// 3：2 * K] = [虚假] *（（N // 6 - （K * K + 4 * K - 2 * K *（ⅰ％2））// 6 - 1）// K + 1）
返回[2，3] + [（3 * I + 1）| 1对于i在范围（1，N- // 3 - 校正）如果筛[Ⅰ]]
smallprimeset =集（primesbelow（100000））
_smallprimeset = 100000
高清isprime（N，precision = 7）：
＃http://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time
如果n＆LT; 1：
提高ValueError错误（出界，第一个参数必须是大于0）
ELIF N'LT; = 3：
返回N'GT; = 2
ELIF N％2 == 0：
返回False
ELIF N'LT; _smallprimeset：
在smallprimeset返回否
D = N - 1
S = 0
而D％2 == 0：
ð// = 2
S + 1 =
对于重复的范围（precision）：
一个= random.randrange（2，N - 2）
X = POW（A，D，N）
如果x == 1或x == N - 1：继续
当r在范围（S - 1）：
X =战俘（X，2，N）
如果x == 1：返回False
如果x == N - 1：突破
其他：返回False
返回True
＃HTTPS：//comeon$c$con.word$p$pss.com/2010/09/18/pollard-rho-brent-integer-factorization/
高清pollard_brent（N）：
如果n％2 == 0：返回2
如果n％3 == 0：返回3
Y，C，M = random.randint（1，N-1），random.randint（1，N-1），random.randint（1，N-1）
克，R，Q = 1，1，1
一边将== 1：
X = Y
因为我在范围（R）：
Y =（POW（Y，2，N）+ C）％N
K = 0
而K＆LT; r和克== 1：
YS = Y
因为我在范围内（分钟（M，R-K））：
Y =（POW（Y，2，N）+ C）％N
Q = Q * ABS（X-Y）％N
G = GCD（Q，N）
第k + = M
R * = 2
如果g = = N：
而真正的：
YS =（POW（YS，2，N）+ C）％N
G = GCD（ABS（X - YS），N）
如果G＆GT; 1：
打破
返回摹
smallprimes = primesbelow（1000）＃似乎较低，但1000 * 1000 = 1000000，因此这将充分因素每个复合＆LT;百万
高清primefactors（N，排序= FALSE）：
因素= []
对检查中smallprimes：
而N％检查== 0：
factors.append（检查）
ñ// =检查
如果检查＆GT; N：突破
如果n＆LT; 2：收益因素
而N'GT; 1：
如果isprime（N）：
factors.append（N）
打破
系数= pollard_brent（N）＃审判庭没有完全的因素，切换到波拉德，布伦特
factors.extend（primefactors（因子））＃改乘因素由波拉德布伦特返回的不一定是主要因素
ñ// =因素
如果排序：factors.sort（）
返回因素
高清分解（N）：
因素= {}
对于P1在primefactors（N）：
尝试：
因素[P1] + = 1
除了KeyError异常：
因素[P1] = 1
返回因素
totients = {}
高清totient（N）：
如果n == 0：返回1
尝试：返回totients [N]
除了KeyError异常：通
TOT = 1
对p，实验中分解（n）的.items（）：
TOT * =（P - 1）* P **（EXP - 1）
totients [N] = TOT
返回TOT
高清GCD（A，B）：
如果== B：返回
而B＆GT; 0：A，B = B，A％B
返回
高清LCM（A，B）：
返回ABS（（一个//最大公约数（A，B））* B）
波拉德，布伦特在用Python实现：
<一个href="https://comeon$c$con.word$p$pss.com/2010/09/18/pollard-rho-brent-integer-factorization/">https://comeon$c$con.word$p$pss.com/2010/09/18/pollard-rho-brent-integer-factorization/
I am looking for an implementation or clear algorithm for getting the prime factorization of N in either python, pseudocode or anything else well-readable. There are a few demands/facts:
- N is between 1 and ~20 digits
- No pre-calculated lookup table, memoization is fine though.
- Need not to be mathematically proven (e.g. could rely on the Goldbach conjecture if needed)
- Need not to be precise, is allowed to be probabilistic/deterministic if needed
I need a fast prime factorization algorithm, not only for itself, but for usage in many other algorithms like calculating the Euler phi(n).
I have tried other algorithms from Wikipedia and such but either I couldn't understand them (ECM) or I couldn't create a working implementation from the algorithm (Pollard-Brent).
I am really interested in the Pollard-Brent algorithm, so any more information/implementations on it would be really nice.
Thanks!
EDIT
After messing around a little I have created a pretty fast prime/factorization module. It combines an optimized trial division algorithm, the Pollard-Brent algorithm, a miller-rabin primality test and the fastest primesieve I found on the internet. gcd is a regular Euclid's GCD implementation (binary Euclid's GCD is much slower then the regular one).
Bounty
Oh joy, a bounty can be acquired! But how can I win it?
- Find an optimalization or bug in my module.
- Provide alternative/better algorithms/implementations.
The answer which is the most complete/constructive gets the bounty.
And finally the module itself:
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)
Pollard-Brent in implemented in Python:
https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
这篇关于快速素因子分解模块的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
