Какой самый быстрый алгоритм факторизации? - PullRequest
54 голосов
/ 15 февраля 2010

Я написал программу, которая пытается найти Дружных Паров. Для этого необходимо найти суммы правильных делителей чисел.

Вот мой текущий 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 огромное число, и оно безумно быстро его разложило.

Что является одним из самых быстрых алгоритмов факторизации?

Ответы [ 7 ]

84 голосов
/ 16 февраля 2010

Вытащил прямо из моего ответа на этот другой вопрос .

Метод будет работать, но будет медленным. "Насколько велики твои числа?" определяет метод для использования:

21 голосов
/ 15 февраля 2010

Вопрос в названии (и в последней строке), похоже, имеет мало общего с фактическим содержанием вопроса. Если вы пытаетесь найти дружественные пары или вычисляете сумму делителей для многих чисел, то раздельное разложение каждого числа (даже с самым быстрым из возможных алгоритмов) - абсолютно неэффективный способ сделать это.

Функция суммы делителей , σ(n) = (sum of divisors of n), является мультипликативной функцией : для относительно простых m и n мы имеем σ(mn) = σ(m)σ(n), поэтому

σ (р 1 к 1 ... р г * +1019 * к г * +1021 ** * тысяча двадцать-два ) = [(p 1 k 1 + 1 -1) / (p 1 -1)]… [(p * * г тысячу тридцать одна тысяча тридцать две * ** * +1033 к г + 1 * +1036 * -1) / (р * +1037 * г -1)].

Таким образом, вы можете использовать любое простое сито (например, расширенную версию Сита Эратосфена ), чтобы найти простых чисел до n, и, в процессе, разложение всех чисел до n. (Например, когда вы делаете сито, сохраняйте наименьший простой множитель каждого n. Затем вы можете позже разложить любое число n путем итерации.) Это будет быстрее (в целом), чем при использовании любого отдельного алгоритм факторизации несколько раз.

Кстати: несколько известных списков дружных пар уже существуют (см., Например, здесь и ссылки на MathWorld ) - поэтому вы пытаетесь расширить запись или делаете это только для весело

15 голосов
/ 15 февраля 2010

Алгоритм Шора: http://en.wikipedia.org/wiki/Shor%27s_algorithm

Конечно, вам нужен квантовый компьютер: D

12 голосов
/ 15 февраля 2010

Я бы предложил начать с того же алгоритма, который используется в Maple, Quadratic Sieve .

  1. Выберите нечетное число n для разложения,
  2. Выберите натуральное число k ,
  3. Поиск по всем p <= <em>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 в yi ,
  6. Используйте Устранение Гаусса , чтобы найти несколько векторов, которые сложены вместе, дают нулевой вектор
  7. Установите x как произведение ri , относящегося к yi , найденному на предыдущем шаге, и установите y как p1 ^ a * p2 ^ b * p3 ^ c * .. * pt ^ z, где показатели степени - это половина показателей, найденных при факторизации yi
  8. Рассчитать d = mcd (xy, n) , если 1 , то d является нетривиальным фактором n , иначе начните с шага 2, выбирая большее k.

Проблема этих алгоритмов состоит в том, что они действительно предполагают много теории в численном исчислении.

6 голосов
/ 15 февраля 2010

Это статья о целочисленной факторизации в Maple.

«Исходя из некоторых очень простых инструкций -« ускорить целочисленную факторизацию в Maple »- мы реализовали алгоритм факторинга Quadratic Sieve в сочетание клена и C ... "

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

4 голосов
/ 27 октября 2015

A 2015 C ++ версия 2 27 реализация таблицы поиска для 1 ГБ памяти:

#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;
}
4 голосов
/ 15 февраля 2010

Зависит от того, насколько велики ваши цифры. Если вы ищете дружественные пары, вы делаете много факторизаций, поэтому ключом может быть не в том, чтобы как можно быстрее вычислить фактор, а в том, чтобы распределить как можно больше работы между различными вызовами. Чтобы ускорить пробное деление, вы можете посмотреть на запоминание и / или предварительный расчет простых чисел до квадратного корня из наибольшего числа, которое вас волнует. Проще получить первичную факторизацию, а затем вычислить сумму всех факторов из этого, чем выполнить цикл до sqrt (n) для каждого числа.

Если вы ищете действительно большие дружные пары, скажем, больше чем 2 ^ 64, то на небольшом количестве машин вы не сможете сделать это, разложив каждое число независимо от того, насколько быстра ваша разложение. Сокращения, которые вы используете для поиска кандидатов, могут помочь вам их учесть.

...