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

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?

Hard to say, without asking them directly but it's possible their implementation is much faster than yours, or they could simply have cached the answer. — Borgleader, Sep 19 '12 at 17:18
Posing Mathematica queries about the 10^n'th for nontrivial values of n convinced me that it had indeed cached answers involving certain large primes. But the OP's question is perhaps not dependent on that, since according to the answers below, it has only one odd factor > 1 (necessarily a prime). — hardmath, Jan 24 '14 at 15:59

score 5 · Accepted Answer · answered Sep 19 '12 at 17:21

5

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...

answered Sep 19 '12 at 17:21

twalberg

59,951
11
89
84

In fact, my laptop, with only a 2GHz processor, and using a pretty naive algorithm similar to what's described above, takes only 58s to determine that 18446744073709551557 (the largest prime that fits in an unsigned 64-bit number) is prime, and less than 0.01s to factor the original number. – twalberg Sep 19 '12 at 17:54
I did not ask for factorization. I speak for finding the 62(including the number itself) divisors of this number. In my laptop with only 2GB RAM and a 1.8GHz dual core processor took me 16seconds to find all the 62 divisors. Using C++ and the pthread libray. If you try to find if a number is prime you do not have to find all the divisors. If you find one divisor then this number is not prime. By the way the number that I speak of is a perfect number that is produced using Mersenne primes. – Dragno Sep 21 '12 at 10:39
The number in the title of the question factors in <0.01s on my laptop, producing 2^30 * (2^31 - 1) - the last factor is prime. Finding all possible divisors is then a pretty simple task of recombining. All divisors either have the form 2^n (where 0<=n<=30) or (2^31 - 1)*2^n (same limits on n), giving 62 possible divisors (including 1 and the number itself). Trying to find all the possible divisors without finding the factors first is a bit brute-force-ish, and not particularly efficient, I think. – twalberg Sep 21 '12 at 14:04

score 1 · Answer 2 · answered Sep 19 '12 at 17:23

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.

score 1 · Answer 3 · answered Oct 08 '12 at 12:33

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

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

3 Answers3