什么是最快的分解算法？

Question

我写了一个试图找到Amicable Pairs的程序.这需要找到数字的适当除数的总和.

这是我目前的sumOfDivisors()方法:

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;
}

所以我需要做很多因子分解,这开始成为我应用程序的真正瓶颈.我在MAPLE中输入了一个巨大的数字,它将它快速地考虑在内.

什么是更快的分解算法？

Answer 1

直接从我对这个问题的回答中拉出来.

该方法可行,但速度很慢."你的数字有多大？" 确定使用方法:

• 小于2 ^ 16左右:查找表.
• 不到2 ^ 70左右:理查德布伦特修改了波拉德的rho算法.
• 小于10 ^ 50:Lenstra椭圆曲线分解
• 小于10 ^ 100:二次筛
• 超过10 ^ 100:General Number Field Sieve

Answer 2

标题中的问题(和最后一行)似乎与问题的实际主体没什么关系.如果你试图找到友好的对,或计算许多数的除数之和,那么单独分解每个数字(即使用最快的算法)绝对是一种低效的方法.

该加总除数函数,?(n) = (sum of divisors of n)是一个乘法功能:相对黄金m和n,我们有?(mn) = ?(m)?(n),所以

σ(p ₁ ^{k ₁} ... p _r ^{k _r})= [(p ₁ ^{k ₁ +1} -1)/(p ₁ -1)] ... [(p _r ^{k _r +1} -1)/(p _r -1 )].

所以你可以使用任何简单的筛子(例如Eratosthenes筛子的增强版本)来找到最多的素数n,并且在此过程中,所有数字的因子分解最多为n.(例如,当您进行筛选时,存储每个n 的最小素数因子.然后您可以n通过迭代对任何数字进行因式分解.)这将比使用任何单独的因子分解算法多次(整体)更快.

BTW:已经存在几个已知友好对的列表(参见例如此处和MathWorld上的链接) - 所以你是试图扩展记录,还是只是为了好玩呢？

Answer 3

Shor的算法:http: //en.wikipedia.org/wiki/Shor%27s_algorithm

当然,你需要一台量子计算机:D

Answer 4

我建议从Maple中使用的相同算法开始,Quadratic Sieve.

1. 选择奇数n进行分解,
2. 选择一个自然数k,
3. 搜索所有p <= k,使k ^ 2与(n mod p)不一致,得到因子基数B = p1,p2,...,pt,
4. 从r > floor(n)开始搜索至少t + 1个值,使y ^ 2 = r ^ 2 - n都只有B中的因子,
5. 对于刚刚计算的每个y1,y2,...,y(t + 1),你生成一个向量v(yi)=(e1,e2,...,et),其中ei是通过减去模2的指数pi来计算的.在y,
6. 使用高斯消元法找到一些加在一起的向量给出一个空向量
7. 将x设置为与上一步中找到的yi相关的ri的乘积,并将y设置为p1 ^ a*p2 ^ b*p3 ^ c*..*pt ^ z其中指数是在分解中找到的指数的一半义
8. 计算 d = mcd(xy,n),如果1 <d <n则d是非平凡因子n,否则从步骤2开始选择更大的k.

这些算法的问题在于它们在数值计算中真正意味着很多理论.

Answer 5

这是Maple中整数分解的一篇论文.

"从一些非常简单的指令开始 - "使Maple中的整数分解更快" - 我们在Maple和C ......的组合中实现了Quadratic Sieve因子分解算法......"

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

Answer 6

取决于你的数字有多大。如果您正在寻找友好的配对，您将进行大量分解，因此关键可能不是尽快分解，而是在不同调用之间共享尽可能多的工作。为了加速试验除法，您可以查看记忆，和/或预先计算素数直到您关心的最大数字的平方根。获得质因数分解，然后从中计算所有因数的总和比对每个数字一直循环到 sqrt(n) 更快。

如果您正在寻找非常大的友好对，比如大于 2^64，那么在少数机器上，无论您的因式分解有多快，您都无法通过对每个数字进行因式分解来做到这一点。您用来寻找候选人的捷径可能会帮助您将他们考虑在内。

Answer 7

1GB 内存的更多 2015 C++ 版本 2 ²⁷查找表实现：

#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;
}

更新：或牺牲一些简单性以在 2 ^{28 之后获得}更大的范围

#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;
}