我试图找到一个巨大整数的除数,我在Haskell中提出了一个问题,但Haskell不够快。我把上面的数字放在 Wolfram Alpha 中,结果立竿见影。这是怎么做到的?
这实际上不是一个困难的分解,因为它是 2^30 * (2^31 - 1)。重复除以 2,直到数字为奇数,然后尝试将2147483647除以大于 2 但小于 sqrt(2147483647)==46340 的奇数进行大约 20k 次迭代。在现代处理器上,这不会花费很长时间...
有很多算法,有些算法对于特定的数字类可以非常快。 例如,如果可以使用快速素数测试算法来确认某个数字是素数,则您就知道其因子。 对于像2305843008139952128这样具有许多小因子的复合数,像Pollard的rho算法这样的东西是快速的。 维基百科列出了许多分解算法,其中许多是特殊用途的。
在像 Wolfram Alpha 这样的一般情况下,一种方法是先尝试许多更快的特殊情况算法,只有在更快的算法找不到答案时,才求助于较慢的保证工作算法。
from random import randint
from fractions import gcd
from math import floor,sqrt
from itertools import combinations
import pdb
def congruent_modulo_n(a,b,n):
return a % n == b % n
def isprimeA(n,k=5):
i=0
while i<k:
a=randint(1,n-1)
if congruent_modulo_n(a**(n-1),1,n) == False:
return False
i=i+1
return True
def powerof2(n):
if n==0: return 1
return 2<<(n-1)
def factoringby2(n):
s=1
d=1
while True:
d=n//powerof2(s)
if isodd(d): break
s=s+1
return (s,d)
def modof2(n):
a0=n>>1
a1=a0<<1
return n-a1
def iseven(m):
return modof2(m)==0
def isodd(m):
return not iseven(m)
class miller_rabin_exception(Exception):
def __init__(self,message,odd=True,lessthan3=False):
self.message=message
self.odd=odd
self.lessthan3=lessthan3
def __str__(self):
return self.message
def miller_rabin_prime_test(n,k=5):
if iseven(n): raise miller_rabin_exception("n must be odd",False)
if n<=3: raise miller_rabin_exception("n must be greater than 3",lessthan3=True)
i=0
s,d=factoringby2(n-1)
z=True
while i<k:
a=randint(2,n-2)
for j in range(0,s):
u=powerof2(j)*d
#x=a**u % n
x=pow(a,u,n)
if x==1 or x==n-1:
z=True
break
else:z=False
i=i+1
return z
def f(x,n):
return pow(x,2,n)+1
def isprime(N):
if N==2 or N==3:
return True
elif N<2:
return False
elif iseven(N):
return False
else:
return miller_rabin_prime_test(N)
def algorithmB(N,outputf):
if N>=2:
#B1
x=5
xx=2
k=1
l=1
n=N
#B2
while(True):
if isprime(n):
outputf(n)
return
while(True):
#B3
g=gcd(xx-x,n)
if g==1:
#B4
k=k-1
if k==0:
xx=x
l=2*l
k=l
x=pow(x,2,n)+1
else:
outputf(g)
if g==n:
return
else:
n=n//g
x=x % n
xx=xx % n
break
def algorithmA(N):
p={}
t=0
k=0
n=N
while(True):
if n==1: return p
for dk in primes_gen():
q,r=divmod(n,dk)
if r!=0:
if q>dk:
continue
else:
t=t+1
p[t]=n
return p
else:
t=t+1
p[t]=dk
n=q
break
def primes_gen():
add=2
yield 2
yield 3
yield 5
p=5
while(True):
p+=add
yield p
if add==2:
add=4
else:
add=2
def algorithmM(a,m,visit):
n=len(a)
m.insert(0,2)
a.insert(0,0)
j=n
while(j!=0):
visit(a[1:])
j=n
while a[j]==m[j]-1:
a[j]=0
j=j-1
a[j]=a[j]+1
def factorization(N):
s=[]
algorithmB(N,s.append)
s1=filter(lambda x:not isprime(x),s)
d=map(algorithmA,s1)
f=map(lambda x:x.values(),d)
r=reduce(lambda x,y:x+y,f,[])+filter(isprime,s)
distinct_factors=set(r)
m=map(r.count,distinct_factors)
return zip(distinct_factors,m)
def divisors(N):
prime_factors=factorization(N)
a=[0]*len(prime_factors)
m=[]
p=[]
for x,y in prime_factors:
m.append(y+1)
p.append(x)
l=[]
algorithmM(a,m,l.append)
result=[]
for x in l:
result.append(reduce(lambda x,y:x*y,map(lambda x,y:x**y,p,x)))
return result