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

Question

我试图找到一个巨大整数的除数我在 Haskell 中提出了一个问题，但 Haskell 不够快。我把上面的数字放在 Wolfram Alpha 中，结果立竿见影。这是怎么做到的？

Answer 1

这实际上并不是一个困难的因式分解，因为它是 2^30 * (2^31 - 1)。重复除以 2 直到数字为奇数，然后大约 20k 次循环尝试将 2147483647 除以大于 2 但小于 sqrt(2147483647)==46340 的奇数。在现代处理器上，这不会花费很长时间......

Answer 2

有很多算法可以进行因式分解，有些算法对于特定类别的数字可能非常快。例如，如果您可以使用快速素数测试算法来确认一个数是否为素数，那么您就知道它的因数。对于像 2305843008139952128 这样有很多小因子的合数，Pollard 的 rho 算法很快。Wikipedia 列出了许多因式分解算法，其中许多是特殊用途的。

在像 Wolfram Alpha 这样的一般情况下，一种方法是首先尝试许多更快的特殊情况算法，如果更快的算法没有找到答案，则只使用较慢的保证工作算法。

Answer 3

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