Wolfram Alpha 如何在不到一秒的时间内计算出2305843008139952128的除数.



我试图找到一个巨大整数的除数,我在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

相关内容

  • 没有找到相关文章

最新更新