Поговорим о факторизации

Как всегда, тут я буду писать пару алгоритмов из Кнута. Второй том. Получисленные алгоритмы.

Я готовлюсь ко третьей части упаковщиков, но этот материал в статью не пойдет (мало чести гордится переписанными на сишные перепевы серьезнейшими алгоритмами).



Огромное спасибо Dmit, за то, что он обратил мое внимание на сие:

http://virlib.eunnet.net/books/numbers/

Relf'у спасибо за то, что он обратил на это внимание Дмита. С этой ссылки началось мое путешествие в теорию чисел.



Рассматривается раздел второго тома - 4.5 - "Разложение на простые множители".



Алгоритм А.

Метод последовательного деления на все простые числа поочередно. Таблица простых чисел генерируется при помощи хитрого решета Эратосфена. Алгоритм немного подкручен - не рассматриваются все четные числа, т.к. все они являются составными. Память тоже отводится по-умному (сам себя не похвалишь, кто похвалит?).

Алгоритм может факторизовать числа лишь под unsigned int максимум. Т.е. не более 4294967295 - 4-х миллионов. Работает довольно быстро.

Код приводится ниже. Код на С++, но, по сути, можно считать, что на С.

Для обработки больших чисел пришлось бы реализовать арифметику бесконечной точности (раздел 4.3 того же тома).

Код (Text):

1.  
2. #include <iostream>
3. #include <math.h>
4.  
5. using namespace std;
6.  
7. inline static unsigned pi(unsigned x)
8. {
9.     return unsigned (x / (log((double)x)-1.08366));
10. }
11.  
12. /**
13.  * Input: M - the number we are creating a sieve up to.
14.  * Output: allocated array of int's (EOF of the array is 0)
15.  *
16.  * Description is here:
17.  * [url=http://mathworld.wolfram.com/SieveofEratosthenes.html
18. ]http://mathworld.wolfram.com/SieveofEratosthenes.html
19. [/url]
20.  * The main idea is to process the odd numbers only.
21.  * Created using the following algo:
22.  * [url=http://okmij.org/ftp/Scheme/Eratosthenes-sieve.txt
23. ]http://okmij.org/ftp/Scheme/Eratosthenes-sieve.txt
24. [/url]
25.  * for efficiency #t and #f are swapped
26.  * For memory saving purposes we roughly try to calculate the
27.  * number of the primes in the given range - pi function.
28.  **
29.  */
30. static unsigned *Eratosthenes_sieve(unsigned M)
31. {
32.     char *v = 0;                // auxiliary array for storing #t or #f
33.     v = new char[M+1];
34.     unsigned i = 0;             // i - index for "v" array, 0 <= i <= M
35.  
36.     unsigned *primes = 0;       // we collect prime numbers here, the EOF marker is 0
37.     primes = new unsigned[pi(M)+1];
38.     primes[0] = 2; 
39.    
40.     unsigned k = 1;             // k - index for "primes" array, 0 <= k <= M
41.  
42.     unsigned p = 0;             // will hold the odd numbers
43.     unsigned j = 0;             // auxiliary variable to get rid of the multiples
44.    
45.     do
46.     {
47.         v[i] = 0;
48.         v[i+1] = 0;
49.         i += 2;
50.     }while(i < M);
51.        
52.     i = 0;
53.  
54.     do
55.     {
56.         p = 2*i + 3;
57.         if (p > M) break;
58.         if (!v[p]) // p is the prime indeed
59.         {
60.             //cout << p << endl;
61.             primes[k] = p;
62.             k++;
63.             // now we remove all the odd multiples
64.             for(j = 1; j <= M; j++)
65.             {
66.                 unsigned tmp = (2*j + 1)*(2*i + 3);
67.  
68.                 if (tmp <= M)
69.                     (v[tmp] = 1);
70.                 else
71.                     break;
72.             }
73.         }
74.  
75.         i++;
76.     }while(i < M);
77.     primes[k] = 0;
78.  
79.     delete v;
80.     return primes;
81. }
82.  
83. /**
84.  * Simple factoring algorithm.
85.  * Donald Knuth. Algorithm 4.5.4A.
86.  **
87.  */
88. void fact_A(unsigned N)
89. {
90.     //t and p are NOT used, we are just printing out the results
91.     unsigned k = 0, n = N;
92.     unsigned *d = Eratosthenes_sieve((unsigned)ceil(sqrt((double)N+1)));
93.    
94.     do
95.     {
96.         unsigned q = n/d[k];
97.         unsigned r = n%d[k];
98.        
99.         if(r)
100.         {
101.             if(q > d[k])
102.             {
103.                 k++;
104.                 continue;
105.             }
106.             else
107.             {
108.                 cout << n << endl;
109.                 break;
110.             }
111.         }
112.         else
113.         {
114.             cout << d[k] << endl;
115.             n = q;
116.         }
117.     } while(n != 1);
118.  
119.     delete d;
120. }
121.  
122.  
123.  
124. void main()
125. {
126.     fact_A(4294967293);
127. }
128.  
129.  

 

Есть еще такая любопытная функция - функция Эйлера. Подробнее функция описана все по той же ссылке:

"§ 3. Важнейшие функции в теории чисел", раздел "Пункт 14. Примеры мультипликативных функций.", "Пример 4. Функция Эйлера.".



Для вычисления этой функции требуется факторизовать число (читай - привести в каноническую форму). Алгоритм fact_A вполне справляется с этой задачей. Сама функция реализуется достаточно просто (только я категорически не согласен с тем, что написано тут:



http://algolist.manual.ru/maths/count_fast/phi_n.php).



Ребята предлагают считать степени, а степени считать для вычислений - это порнография в данном случае, т.к. есть другая формула, приведенная в уже упомянутой книжке.



Прошу прощения за использование С++ - в этот раз действительно используется С++ и STL, т.к. мне было лениво писать собственную реализацию для того, чтобы получить уникальные значения простых чисел (если непонятно, то разложите число 50 на простые множители - это будет 5*5*2 - так вот, функция Эйлера требует только 5 и 2, т.е. вторая пятерка должна быть исключена).



Поэтому имеем:

Код (Text):

1.  
2. double Euler_totient(unsigned x)
3. {
4.     vector <unsigned> factored;
5.    
6.     fact_A(x, factored);
7.    
8.     vector <unsigned>::iterator new_end = unique(factored.begin(), factored.end());
9.     factored.erase(new_end, factored.end());
10.  
11.     double n = x;
12.  
13.     for(vector <unsigned>::iterator i = factored.begin(); i != factored.end(); i++)
14.     {
15.         n *= (1 - 1.0/(*i));
16.     }
17.    
18.     return n;
19. }
20.  

Например, для 1000 такая функция даст 400. Это означает, что есть 400 простых чисел в интервале от 0 до 1000.

Также fact_A-алгоритм немного модифицирован. Теперь он принимает ссылку на вектор.

 

volodya, Вы немного не правы. Ф-ция Эйлера возвращает количество натуральных чисел, меньших ее аргумента и взаимно простых с ним. Действительно, EulerPhi[1000] == 400. Но количество простых чисел в интервале [1, 1000] равно pi(1000) == 168.



P.S. Спасибо за приведенную ссылку.

P.P.S. Очень рекомендую пакет Mathematica - IMHO весьма полезный для реверсинга инструмент

 

volodya, Вы немного не правы.



Ох, да Спасибо, что поправил!



"Функция Эйлера phi( a ) есть количество чисел из ряда 0, 1, 2,..., a - 1, взаимно простых с a ."



Да, Mathematica, штука, конечно, хорошая, но мне нравится Maple

 

[offtop]

volodya

Maple? Честно говоря, я боюсь её использовать. В ней стоолько багов, причём таких откровенных, что зачастую приходится пересчитывать в других пакетах…

 

Ладно, продолжаем. Сегодня у нас в гостях метод Полларда (Pollard Rho method).



Линки на теоретический ликбез:

http://algolist.manual.ru/maths/teornum/factor/p-1.php

http://mathworld.wolfram.com/PollardRhoFactorizationMethod.html





На

http://algolist.manual.ru/maths/teornum/factor/

есть линк на "A Survey of Modern Integer Factorization Methods" - подберите.



Линки на понажимать кнопочки:

http://linguistlist.org/~zheng/courseware/pollard.html

http://www.louisville.edu/~dawill03/crypto/PollardRho.html



Для Поллардовского метода нам понадобятся некоторые дополнительные функции.

1) gcd - алгоритм Евклида:

http://www.wasm.ru/forum/index.php?action=vthread&forum=17&topic=6839



2) вероятностный тест Миллера-Рабина на простоту числа - взят из Кнута и дополнительно требует:

2.1) быстрое возведение в степень

2.2) представление числа в форме n - 1 = 2^kq:

http://www.wasm.ru/forum/index.php?action=vthread&forum=17&topic=6841



Теперь мы готовы. Код:

Код (Text):

1.  
2. static unsigned powmod(unsigned a, unsigned k, unsigned n)
3. {
4.   unsigned b=1;
5.  
6.   while (k)
7.   {
8.       if (!(k%2))
9.       {
10.           k /= 2;
11.           a = (a*a)%n;
12.       }
13.       else
14.       {
15.           k--;
16.           b = (b*a)%n;
17.       }
18.   }
19.   return b;
20. }
21.  
22. /* Donald Knuth, vol 2., algorithm 4.5.4P - probability test for prime*/
23. static bool isprime(unsigned n)
24. {
25.     srand(time(0));
26.  
27.     unsigned x = rand() % (n - 1);
28.         unsigned k, q;
29.     do_N(n, k, q);
30.  
31.     unsigned j = 0;
32.     unsigned y = powmod(x, q, n);
33.    
34.     while(j < k)
35.     {
36.             if(((!j) && (y == 1)) || (y == n - 1))
37.             return true;
38.         if((j > 0) && (y == 1))
39.             return false;
40.  
41.         j++;
42.         y = (y*y)%n;
43.         }
44.  
45.     return false;
46. }
47.  
48. static unsigned pollard_rho(unsigned n)
49. {
50.     if(isprime(n))
51.         return 0;
52.  
53.     unsigned factor = 1, x = 5, y = 2;
54.  
55.     while(gcd(x-y, n) == 1)
56.     {
57.         x = (x*x + 1)%n;
58.         x = (x*x + 1)%n;
59.         y = (y*y + 1)%n;
60.     }
61.    
62.     factor = gcd(x-y, n);
63.     if(factor < n)
64.         return factor;
65.  
66.     return 0;
67. }
68.  
69.  
70. void main()
71. {
72.     cout << pollard_rho(561);
73. }
74.  
75.  

 

Еще хорошая ссылочка:



http://www.frenchfries.net/paul/factoring/theory/index.html

 

Метод разложения на множители Ферма лучше всего описан тут:

http://www.wasm.ru/forum/index.php?action=vthread&forum=18&topic=6661



Эту книгу стоит купить. В виде алгоритма, что предлагает Кнут, это может выглядеть так:

Код (Text):

1. #include <iostream>
2. #include <math.h>
3.  
4. using namespace std;
5.  
6. /* 4.5.4C */
7. void fermat(unsigned n)
8. {
9.     double x = 2*floor(sqrt((double)n)) + 1;
10.     unsigned y = 1;
11.     double r = pow((floor(sqrt((double)n))), 2) - n;
12.  
13.     for(;;)
14.     {
15.         if(!r)
16.         {
17.             cout << (x-y)/2 << endl;
18.             return;
19.         }
20.  
21.         r += x;
22.         x += 2;
23.  
24.         do
25.         {
26.             r -= y;
27.             y += 2;
28.         }while(r > 0);
29.     }
30. }
31.  
32.  
33. void main()
34. {
35.     fermat(377);
36. }

 

В этом разделе ранее был приведен популярный вероятностный тест Рабина-Миллера на простоту числа (был взят из Кнута). Не так давно тройка веселых индусов открыла куда более крутой детерминистический тест. Тест очень быстр. Хм. Но и не слишком-то прост

Реализацию теста можно, например, найти в Miracl.

Сама статья находится по адресу:



http://www.cse.iitk.ac.in/news/primality.html.



Чтобы написать этот алгоритм на С/С++ самому, вам потребуется большинство из процедур, приведенных в этом топике выше, а также процедура проверки числа на степень. Процедуру можно найти тут:



http://www.wasm.ru/forum/index.php?action=vthread&forum=17&topic=6888

 

http://www.cse.iitk.ac.in/news/primality.html Error 404

 

Loger

Как ни странно, после поиска в Гугле и закачки

старого варианта алгоритма

http://www.cse.iitk.ac.in/primality.pdf

страничка опять появилась. У них там все через базу данных организовано - видимо свой кеш почистили....

 

В ссылке от volodya лишняя точка (после html), посему у Loger и Error 404 появлялось )

 

Код (Text):

1. static unsigned pollard_rho(unsigned n)
2. {
3. ...
4. return 0;
5. ...
6. }

Ноль возвращается, когда алгоритм не смог найти делителя (возможно, число простое). Не подскажете, что делать в таких случаях? Или, может быть, есть более продвинутые алгоритмы (мне нужно порядка корня - корня 4-й степени из N).

 

Перед факторизацией надо прогонять через тест на простоту.

 

maxdiver

Попробовать с другими параметрами и/или другим начальным значением. На простоту число всегда нужно проверять предварительно и пытаться разложить, только если известно, что оно составное.

Есть, но они значительно сложнее в реализации. Правда и скорость больше. Смотри различные варианты quadratic sieve(стандарт для не очень больших чисел) и number field sieve.

 

Если, кто читал последнюю версию моего мана http://xproject-all.narod.ru/factorize_numbers_rus_descrip_ver6.pdf там мной была допущена ошибка в рассуждениях по MOD CATCHER.
а, вообще, замечу, что официально известные алгосы факторизации смехотворны для нумберов на 1024 и более бит.

 

да, ну maple кульная штука, а вот mathematica имеет забавную особенность сильно грузить проц при простое)

 

Ну раз уж подняли тему...
Ссылки передохли, но вот то, шо надо:
Алгоритм называется AKS. Он далеко не самый быстрый, да и сами авторы утверждают: "Certain algorithms need to generate prime numbers in order to construct cryptographic keys, but algorithms to accomplish this which can be executed very efficiently already existed ... These algorithms that are commonly used in practice are actually faster than the ones [AKS] proposed"
тут:
http://web.archive.org/web/20060528185603/http://crypto.cs.mcgill.ca/~stiglic/PRIMES_P_FAQ.html

 

CreatorCray

Желательно пользоваться поиском, AKS как раз недавно обсуждали.

http://www.wasm.ru/forum/viewtopic.php?pid=155502#p155502
http://www.wasm.ru/forum/viewtopic.php?pid=155593#p155593

 

Я попытался сделать так: рандомная функция это f(x)=x*x+C (mod) N. И если произошла неудача, то увеличиваем C на единицу. Но это же очень коряво и всё равно тормознуто. Ведь мы так можем перебирать значения C прямо до корня из N, т.е. до бесконечности.
Возможно, нужно улучшить метод. Скажем, на algolist'е предлагается навороченный метод p-1 из двух шагов. Правда, я не смог его реализовать

А какие быстрые методы проверки на простоту? Я имею в виду не вероятностные, как isprime() в посте выше, а четкие тесты на простоту (в Wikipedia нашёл какие-то, но проблемы с реализацией).