快速素分解模块

我正在寻找一个实现或明确的algorithm来获得N的主要因式分解在python，伪代码或其他任何可读的。 有一些要求/事实：

• N在1到20个数字之间
• 没有预先计算的查找表，memoization是好的，但。
• 不需要mathcertificate（如果需要，可以依靠哥德巴赫猜想）
• 如果需要的话，不需要是精确的，可以是概率的/确定的

我需要一个快速素数因子分解algorithm，不仅适用于本身，还适用于许多其他algorithm，如计算Euler phi（n） 。

我尝试过维基百科等其他algorithm，但是我无法理解它们（ECM），或者我无法从algorithm（Pollard-Brent）创build工作实现。

我对Pollard-Brentalgorithm非常感兴趣，所以更多的信息/实现将会非常好。

谢谢！

编辑

经过一段时间后，我创build了一个相当快的素数/分解模块。 它结合了优化的试algorithm，Pollard-Brentalgorithm，miller-rabin素数testing和我在互联网上find的速度最快的primesieve。 gcd是一个规则的Euclid的GCD实现（二进制Euclid的GCD 比常规的要慢得多）。

赏金

哦，快乐，可以获得赏金！ 但是我怎么能赢呢？

• 在我的模块中find最优化或错误。
• 提供替代/更好的algorithm/实现。

最完整/有build设性的答案得到赏赐。

最后模块本身：

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, rk)): y = (pow(y, 2, n) + c) % n q = q * abs(xy) % 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

在Python中实现Pollard-Brent：

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

Answer 2

如果你不想重新发明轮子，可以使用库sympy

 pip install sympy 

使用函数sympy.ntheory.factorint

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

没有必要使用primesbelow计算primesbelow ，使用smallprimeset 。

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

把你的primefactors分成两个函数来处理smallprimes和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的validation结果，可以减less使用Miller-Rabin素数testing的isprime的重复。 来自Wiki 。

• 如果n <1,373,653，则testinga = 2和3就足够了;
• 如果n <9,080,191，则足以testinga = 31和73;
• 如果n <4,759,123,141，则足以testinga = 2，7和61;
• 如果n <2,152,302,898,747，则足以testinga = 2,3,5,7和11;
• 如果n <3,474,749,660,383，则足以testinga = 2,3,5,7,11和13;
• 如果n <341,550,071,728,321，则足以testinga = 2,3,5,7,11,13和17。

编辑1 ：更正了if-else返callback用，以便将主因子追加到因子中。

Answer 4

即使在目前的情况下，也有几个地方需要注意。

1. 每个循环不要做checker*checker ，使用s=ceil(sqrt(num))和checher < s
2. checher每次加2，忽略2以外的所有偶数
3. 使用divmod而不是%和//

Answer 5

你应该可以做一些素数检测，你可以看看这里， 寻找素数的快速algorithm？

您应该阅读整个博客，他列出了一些testing素数的algorithm。

Answer 6

有一个python库有一些素性testing（包括不能做的不正确的）。 这叫做pyprimes 。 认为值得一提的是后人的目的。 我不认为它包括你提到的algorithm。

Answer 7

在考虑数字2**1427 * 31时，我刚碰到这个代码中的一个bug。

  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 

此代码片段：

 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的计算，这对分解数字如2**1427 。 没有理由计算sqrt(2**1427) ， sqrt(2**1426) ， sqrt(2**1425)等等，当我们所关心的是“检查器的平方超过n “。

重写的代码在出现大数字时不会引发exception，根据时间的不同（大约是样本input2和2**718 * 31 ）的两倍。

另外注意isprime(2)返回的结果是错误的，但只要我们不依赖它就可以了。 恕我直言，你应该重写该函数的介绍

 if n <= 3: return n >= 2 ... 

Answer 8

你可以分解到一个极限，然后使用布伦特来获得更高的因子

 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,rk)): y = ((y*y)%N+c)%N q = q*(abs(xy))%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))