How can Wolfram Alpha compute the divisors of 2305843008139952128 in less than a second? [closed]

Question

Closed. This question needs to be more focused. It is not currently accepting answers.

---

Want to improve this question? Update the question so it focuses on one problem only by editing this post.

Closed 1 year ago.

Improve this question

I am trying to find the divisors of a huge integer I have made a question about that in Haskell but Haskell is not fast enough. I put the above number in Wolfram Alpha and the result was immediate. How this was done?

Answer 1

That's not actually a difficult factorization, since it's 2^30 * (2^31 - 1). Repeated division by two until the number is odd, then around 20k iterations of a loop attempting to divide 2147483647 by odd numbers greater than 2 but less than sqrt(2147483647)==46340. On modern processors, that's not going to take very long...

Answer 2

There are lots of algorithms for factoring, and some can be very quick for particular classes of numbers. For example, if you can use a fast prime-testing algorithm to confirm that a number is prime, you know its factors. For a composite number like 2305843008139952128 which has many small factors, something like Pollard's rho algorithm is fast. Wikipedia lists a number of factoring algorithms, many of which are special purpose.

In the general case like Wolfram Alpha, one approach would be to try many of the faster special-case algorithms first, and only resort to slower guaranteed-to-work algorithms if the faster ones don't find an answer.

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