wolframalpha - Wolfram Alphaはどのようにして2305843008139952128の約数を1秒以内に計算できますか?

Question

Haskellでそれについて質問した巨大な整数の約数を見つけようとしていますが、Haskellは十分に高速ではありません. 上記の数値を Wolfram Alpha に入力すると、すぐに結果が得られました。これはどのように行われたのですか？

Answer 1

これは 2^30 * (2^31 - 1) なので、実際には難しい因数分解ではありません。数が奇数になるまで 2 による除算を繰り返し、次に 2147483647 を 2 より大きく sqrt(2147483647)==46340 未満の奇数で除算しようとするループの約 20k 反復。最新のプロセッサでは、それほど長くはかかりません...

Answer 2

因数分解には多くのアルゴリズムがあり、特定のクラスの数に対して非常に高速なアルゴリズムもあります。たとえば、高速な素数テスト アルゴリズムを使用して数値が素数であることを確認できれば、その因数がわかります。2305843008139952128 のような多数の小さな因数を持つ合成数の場合、Pollard の rho アルゴリズムのようなものが高速です。ウィキペディアには多くのファクタリング アルゴリズムがリストされており、その多くは特殊な目的のものです。

Wolfram Alpha のような一般的なケースでは、最初に多くの高速な特殊ケース アルゴリズムを試して、高速なアルゴリズムで答えが見つからない場合にのみ、動作が保証された低速なアルゴリズムに頼るという 1 つの方法があります。

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