What is the fastest factorization algorithm?

Question

I wrote a program that tries to find Amicable Pairs. This requires finding the sums of the proper divisors of numbers.

Here is my current sumOfDivisors() method:

 int sumOfDivisors(int n) { int sum = 1; int bound = (int) sqrt(n); for(int i = 2; i <= 1 + bound; i++) { if (n % i == 0) sum = sum + i + n / i; } return sum; } 

So I need to do a lot of factorization, and it starts to become a real bottleneck in my application. I scored a huge amount in MAPLE, and it produced insanely fast.

What is one of the fastest factorization algorithms?

Answer 1

Pulled directly from my answer to this other question .

The method will work, but will be slow. "How big are your numbers?" defines the method of use:

• Less than 2 ^ 16 or so: Search table.
• Less than 2 ^ 70 or so: A modification of Richard Brent Pollard rho .
• Less than 10 ^ 50: Fragmentation of an elliptic Lenstra curve
• Less than 10 ^ 100: Quadratic strainer
• Over 10 ^ 100: General field sieve with numbers

Answer 2

The question in the title (and on the last line) seems to have little in common with the actual body of the question. If you are trying to find friendly pairs or calculate the sum of divisors for many numbers, then separately factoring each number (even with the fastest algorithm possible), this is an absolutely inefficient way to do this.

The sum-divisor function, σ(n) = (sum of divisors of n) , is a multiplicative function : for relatively simple m and n we have σ(mn) = σ(m)σ(n) , therefore

σ (p <sub> 1sub> ^{to <sub> 1sub>} ... p <sub> tsub> ^{to <sub> t} ) = [(p ₁ ^{k ₁ +1} -1) / (p ₁ -1)]. .. [(p <sub> tsub> ^{to <south> tsub> + 1} -1) / (p <sub> tsub> - 1)].

So, you would use any simple sieve (for example, an extended version of Sieve of Eratosthenes ) to find primes up to n , while factoring all numbers up to n. (For example, as you make your sieve, keep the smallest primary coefficient for each n. Then you can factor the number of numbers n by repeating.) This would be faster (in general) than using any particular factorization algorithm several times.

By the way: several well-known lists of friendly couples already exist (see for example here and links to MathWorld ) - are you trying to prolong the recording or just do it for fun?

Answer 3

Shore Algorithm: http://en.wikipedia.org/wiki/Shor%27s_algorithm

Of course you need a quantum computer: D

Answer 4

I would suggest starting with the same algorithm used in Maple, Quadratic Sieve .

• Choose an odd number n to factorize,
• Choose a natural number k,
• Find all p & lt = k so that k ^ 2 does not match (n mod p) to get the factor base B = p1, p2, ..., pt,
• Starting with r> floor (n), find at least t + 1 values, so that y ^ 2 = r ^ 2 - n only factors in B have,
• For each calculated y1, y2, ..., y (t + 1) you create a vector v (yi) = (e1, e2, ..., et), where ei is calculated by decreasing modulo 2 the index pi in yi,
• Use Gaussian Elimination to find some of the vectors added together, give a zero vector
• Define x as the product ri associated with yi found in the previous step and set y as p1 ^ a * p2 ^ b * p3 ^ c * .. * pt ^ z, where the exponents are half the exponents found in the factorization yi
• Calculate d = mcd (xy, n), if 1 <d <n, then d is a nontrivial factor of n, otherwise start from step 2, choosing a larger k.

The problem with these algorithms is that they really imply a lot of theory in numerical terms.

Answer 5

This is a factorization document for maple integers.

"Starting with some very simple instructions -" speed up integer factorization in Maple "- we have implemented a quadratic sieve factorization algorithm in a combination of Maple and C ..."

http://www.cecm.sfu.ca/~pborwein/MITACS/papers/percival.pdf

Answer 6

Depends on how big your numbers are. If you are looking for friendly couples, you do a lot of factorization, so the key may not be as fast as possible, but share as much work as possible between different calls. To speed up trial division, you can look at the memorandum and / or preliminary calculation of primes to the square root of the largest number that you care about. It is faster to get the basic factorization and then calculate the sum of all the factors from this than to ensure that all the loops down to sqrt (n) for each number.

If you are looking for really big friendly pairs, say more than 2 ^ 64, then on a small number of machines you cannot do this by expanding each number no matter how fast your factorization. The abbreviations you use to search for candidates can help you identify them.

Answer 7

More than 2015 C ++ version 2 ²⁷ implementation of the lookup table for 1 GB of memory:

 #include <iostream.h> // cerr, cout, and NULL #include <string.h> // memcpy() #define uint unsigned __int32 uint *factors; const uint MAX_F=134217728; // 2^27 void buildFactors(){ factors=new (nothrow) uint [(MAX_F+1)*2]; // 4 * 2 * 2^27 = 2^30 = 1GB if(factors==NULL)return; // not able to allocate enough free memory int i; for(i=0;i<(MAX_F+1)*2;i++)factors[i]=0; //Sieve of Eratosthenese factors[1*2]=1; factors[1*2+1]=1; for(i=2;i*i<=MAX_F;i++){ for(;factors[i*2] && i*i<=MAX_F;i++); factors[i*2]=1; factors[i*2+1]=i; for(int j=2;i*j<=MAX_F;j++){ factors[i*j*2]=i; factors[i*j*2+1]=j; } } for(;i<=MAX_F;i++){ for(;i<=MAX_F && factors[i*2];i++); if(i>MAX_F)return; factors[i*2]=1; factors[i*2+1]=i; } } uint * factor(uint x, int &factorCount){ if(x > MAX_F){factorCount=-1;return NULL;} uint tmp[70], at=x; int i=0; while(factors[at*2]>1){ tmp[i++]=factors[at*2]; cout<<"at:"<<at<<" tmp:"<<tmp[i-1]<<endl; at=factors[at*2+1]; } if(i==0){ cout<<"at:"<<x<<" tmp:1"<<endl; tmp[i++]=1; tmp[i++]=x; }else{ cout<<"at:"<<at<<" tmp:1"<<endl; tmp[i++]=at; } factorCount=i; uint *ret=new (nothrow) uint [factorCount]; if(ret!=NULL) memcpy(ret, tmp, sizeof(uint)*factorCount); return ret; } void main(){ cout<<"Loading factors lookup table"<<endl; buildFactors(); if(factors==NULL){cerr<<"Need 1GB block of free memory"<<endl;return;} int size; uint x=30030; cout<<"\nFactoring: "<<x<<endl; uint *f=factor(x,size); if(size<0){cerr<<x<<" is too big to factor. Choose a number between 1 and "<<MAX_F<<endl;return;} else if(f==NULL){cerr<<"ran out of memory trying to factor "<<x<<endl;return;} cout<<"\nThe factors of: "<<x<<" {"<<f[0]; for(int i=1;i<size;i++) cout<<", "<<f[i]; cout<<"}"<<endl; delete [] f; x=30637; cout<<"\nFactoring: "<<x<<endl; f=factor(x,size); cout<<"\nThe factors of: "<<x<<" {"<<f[0]; for(int i=1;i<size;i++) cout<<", "<<f[i]; cout<<"}"<<endl; delete [] f; delete [] factors; } 

Update: or sacrifice some simplicity for a slightly larger range in just 2 ²⁸

 #include <iostream.h> // cerr, cout, and NULL #include <string.h> // memcpy(), memset() //#define dbg(A) A #ifndef dbg #define dbg(A) #endif #define uint unsigned __int32 #define uint8 unsigned __int8 #define uint16 unsigned __int16 uint * factors; uint8 *factors08; uint16 *factors16; uint *factors32; const uint LIMIT_16 = 514; // First 16-bit factor, 514 = 2*257 const uint LIMIT_32 = 131074;// First 32-bit factor, 131074 = 2*65537 const uint MAX_FACTOR = 268501119; //const uint64 LIMIT_64 = 8,589,934,594; // First 64-bit factor, 2^33+1 const uint TABLE_SIZE = 268435456; // 2^28 => 4 * 2^28 = 2^30 = 1GB 32-bit table const uint o08=1, o16=257 ,o32=65665; //o64=4294934465 // TableSize = 2^37 => 8 * 2^37 = 2^40 1TB 64-bit table // => MaxFactor = 141,733,953,600 /* Layout of factors[] array * Indicies(32-bit) i Value Size AFactorOf(i) * ---------------- ------ ---------- ---------------- * factors[0..128] [1..513] 8-bit factors08[i-o08] * factors[129..65408] [514..131073] 16-bit factors16[i-o16] * factors[65409..268435455] [131074..268501119] 32-bit factors32[i-o32] * * Note: stopping at i*i causes AFactorOf(i) to not always be LargestFactor(i) */ void buildFactors(){ dbg(cout<<"Allocating RAM"<<endl;) factors=new (nothrow) uint [TABLE_SIZE]; // 4 * 2^28 = 2^30 = 1GB if(factors==NULL)return; // not able to allocate enough free memory uint i,j; factors08 = (uint8 *)factors; factors16 = (uint16 *)factors; factors32 = factors; dbg(cout<<"Zeroing RAM"<<endl;) memset(factors,0,sizeof(uint)*TABLE_SIZE); //for(i=0;i<TABLE_SIZE;i++)factors[i]=0; //Sieve of Eratosthenese //8-bit values dbg(cout<<"Setting: 8-Bit Values"<<endl;) factors08[1-o08]=1; for(i=2;i*i<LIMIT_16;i++){ for(;factors08[i-o08] && i*i<LIMIT_16;i++); dbg(cout<<"Filtering: "<<i<<endl;) factors08[i-o08]=1; for(j=2;i*j<LIMIT_16;j++)factors08[i*j-o08]=i; for(;i*j<LIMIT_32;j++)factors16[i*j-o16]=i; for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i; } for(;i<LIMIT_16;i++){ for(;i<LIMIT_16 && factors08[i-o08];i++); dbg(cout<<"Filtering: "<<i<<endl;) if(i<LIMIT_16){ factors08[i-o08]=1; j=LIMIT_16/i+(LIMIT_16%i>0); for(;i*j<LIMIT_32;j++)factors16[i*j-o16]=i; for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i; } }i--; dbg(cout<<"Setting: 16-Bit Values"<<endl;) //16-bit values for(;i*i<LIMIT_32;i++){ for(;factors16[i-o16] && i*i<LIMIT_32;i++); factors16[i-o16]=1; for(j=2;i*j<LIMIT_32;j++)factors16[i*j-o16]=i; for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i; } for(;i<LIMIT_32;i++){ for(;i<LIMIT_32 && factors16[i-o16];i++); if(i<LIMIT_32){ factors16[i-o16]=1; j=LIMIT_32/i+(LIMIT_32%i>0); for(;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i; } }i--; dbg(cout<<"Setting: 32-Bit Values"<<endl;) //32-bit values for(;i*i<=MAX_FACTOR;i++){ for(;factors32[i-o32] && i*i<=MAX_FACTOR;i++); factors32[i-o32]=1; for(j=2;i*j<=MAX_FACTOR;j++)factors32[i*j-o32]=i; } for(;i<=MAX_FACTOR;i++){ for(;i<=MAX_FACTOR && factors32[i-o32];i++); if(i>MAX_FACTOR)return; factors32[i-o32]=1; } } uint * factor(uint x, int &factorCount){ if(x > MAX_FACTOR){factorCount=-1;return NULL;} uint tmp[70], at=x; int i=0; while(at>=LIMIT_32 && factors32[at-o32]>1){ tmp[i++]=factors32[at-o32]; dbg(cout<<"at32:"<<at<<" tmp:"<<tmp[i-1]<<endl;) at/=tmp[i-1]; } if(at<LIMIT_32){ while(at>=LIMIT_16 && factors16[at-o16]>1){ tmp[i++]=factors16[at-o16]; dbg(cout<<"at16:"<<at<<" tmp:"<<tmp[i-1]<<endl;) at/=tmp[i-1]; } if(at<LIMIT_16){ while(factors08[at-o08]>1){ tmp[i++]=factors08[at-o08]; dbg(cout<<"at08:"<<at<<" tmp:"<<tmp[i-1]<<endl;) at/=tmp[i-1]; } } } if(i==0){ dbg(cout<<"at:"<<x<<" tmp:1"<<endl;) tmp[i++]=1; tmp[i++]=x; }else{ dbg(cout<<"at:"<<at<<" tmp:1"<<endl;) tmp[i++]=at; } factorCount=i; uint *ret=new (nothrow) uint [factorCount]; if(ret!=NULL) memcpy(ret, tmp, sizeof(uint)*factorCount); return ret; } uint AFactorOf(uint x){ if(x > MAX_FACTOR)return -1; if(x < LIMIT_16) return factors08[x-o08]; if(x < LIMIT_32) return factors16[x-o16]; return factors32[x-o32]; } void main(){ cout<<"Loading factors lookup table"<<endl; buildFactors(); if(factors==NULL){cerr<<"Need 1GB block of free memory"<<endl;return;} int size; uint x=13855127;//25255230;//30030; cout<<"\nFactoring: "<<x<<endl; uint *f=factor(x,size); if(size<0){cerr<<x<<" is too big to factor. Choose a number between 1 and "<<MAX_FACTOR<<endl;return;} else if(f==NULL){cerr<<"ran out of memory trying to factor "<<x<<endl;return;} cout<<"\nThe factors of: "<<x<<" {"<<f[0]; for(int i=1;i<size;i++) cout<<", "<<f[i]; cout<<"}"<<endl; delete [] f; x=30637; cout<<"\nFactoring: "<<x<<endl; f=factor(x,size); cout<<"\nThe factors of: "<<x<<" {"<<f[0]; for(int i=1;i<size;i++) cout<<", "<<f[i]; cout<<"}"<<endl; delete [] f; delete [] factors; }