|
| Посл.отвђт | Сообщенiе |
|
|
Дата: Авг 13, 2004 22:25:55 Как всегда, тут я буду писать пару алгоритмов из Кнута. Второй том. Получисленные алгоритмы. Я готовлюсь ко третьей части упаковщиков, но этот материал в статью не пойдет (мало чести гордится переписанными на сишные перепевы серьезнейшими алгоритмами). Огромное спасибо Dmit, за то, что он обратил мое внимание на сие: http://virlib.eunnet.net/books/numbers/ Relf'у спасибо за то, что он обратил на это внимание Дмита. С этой ссылки началось мое путешествие в теорию чисел. Рассматривается раздел второго тома - 4.5 - "Разложение на простые множители". Алгоритм А. Метод последовательного деления на все простые числа поочередно. Таблица простых чисел генерируется при помощи хитрого решета Эратосфена. Алгоритм немного подкручен - не рассматриваются все четные числа, т.к. все они являются составными. Память тоже отводится по-умному (сам себя не похвалишь, кто похвалит?). Алгоритм может факторизовать числа лишь под unsigned int максимум. Т.е. не более 4294967295 - 4-х миллионов. Работает довольно быстро. Код приводится ниже. Код на С++, но, по сути, можно считать, что на С. Для обработки больших чисел пришлось бы реализовать арифметику бесконечной точности (раздел 4.3 того же тома).
#include <iostream>
#include <math.h>
using namespace std;
inline static unsigned pi(unsigned x)
{
return unsigned (x / (log((double)x)-1.08366));
}
/**
* Input: M - the number we are creating a sieve up to.
* Output: allocated array of int's (EOF of the array is 0)
*
* Description is here:
* http://mathworld.wolfram.com/SieveofEratosthenes.html
* The main idea is to process the odd numbers only.
* Created using the following algo:
* http://okmij.org/ftp/Scheme/Eratosthenes-sieve.txt
* for efficiency #t and #f are swapped
* For memory saving purposes we roughly try to calculate the
* number of the primes in the given range - pi function.
**
*/
static unsigned *Eratosthenes_sieve(unsigned M)
{
char *v = 0; // auxiliary array for storing #t or #f
v = new char[M+1];
unsigned i = 0; // i - index for "v" array, 0 <= i <= M
unsigned *primes = 0; // we collect prime numbers here, the EOF marker is 0
primes = new unsigned[pi(M)+1];
primes[0] = 2;
unsigned k = 1; // k - index for "primes" array, 0 <= k <= M
unsigned p = 0; // will hold the odd numbers
unsigned j = 0; // auxiliary variable to get rid of the multiples
do
{
v[i] = 0;
v[i+1] = 0;
i += 2;
}while(i < M);
i = 0;
do
{
p = 2*i + 3;
if (p > M) break;
if (!v[p]) // p is the prime indeed
{
//cout << p << endl;
primes[k] = p;
k++;
// now we remove all the odd multiples
for(j = 1; j <= M; j++)
{
unsigned tmp = (2*j + 1)*(2*i + 3);
if (tmp <= M)
(v[tmp] = 1);
else
break;
}
}
i++;
}while(i < M);
primes[k] = 0;
delete v;
return primes;
}
/**
* Simple factoring algorithm.
* Donald Knuth. Algorithm 4.5.4A.
**
*/
void fact_A(unsigned N)
{
//t and p are NOT used, we are just printing out the results
unsigned k = 0, n = N;
unsigned *d = Eratosthenes_sieve((unsigned)ceil(sqrt((double)N+1)));
do
{
unsigned q = n/d[k];
unsigned r = n%d[k];
if(r)
{
if(q > d[k])
{
k++;
continue;
}
else
{
cout << n << endl;
break;
}
}
else
{
cout << d[k] << endl;
n = q;
}
} while(n != 1);
delete d;
}
void main()
{
fact_A(4294967293);
}
|
|
|
Дата: Авг 14, 2004 01:20:24 · Поправил: volodya Есть еще такая любопытная функция - функция Эйлера. Подробнее функция описана все по той же ссылке: "§ 3. Важнейшие функции в теории чисел", раздел "Пункт 14. Примеры мультипликативных функций.", "Пример 4. Функция Эйлера.". Для вычисления этой функции требуется факторизовать число (читай - привести в каноническую форму). Алгоритм fact_A вполне справляется с этой задачей. Сама функция реализуется достаточно просто (только я категорически не согласен с тем, что написано тут: http://algolist.manual.ru/maths/count_fast/phi_n.php). Ребята предлагают считать степени, а степени считать для вычислений - это порнография в данном случае, т.к. есть другая формула, приведенная в уже упомянутой книжке. Прошу прощения за использование С++ - в этот раз действительно используется С++ и STL, т.к. мне было лениво писать собственную реализацию для того, чтобы получить уникальные значения простых чисел (если непонятно, то разложите число 50 на простые множители - это будет 5*5*2 - так вот, функция Эйлера требует только 5 и 2, т.е. вторая пятерка должна быть исключена). Поэтому имеем:
double Euler_totient(unsigned x)
{
vector <unsigned> factored;
fact_A(x, factored);
vector <unsigned>::iterator new_end = unique(factored.begin(), factored.end());
factored.erase(new_end, factored.end());
double n = x;
for(vector <unsigned>::iterator i = factored.begin(); i != factored.end(); i++)
{
n *= (1 - 1.0/(*i));
}
return n;
}
Например, для 1000 такая функция даст 400. Это означает, что есть 400 простых чисел в интервале от 0 до 1000. Также fact_A-алгоритм немного модифицирован. Теперь он принимает ссылку на вектор. |
|
|
Дата: Авг 14, 2004 01:46:25 „Например, для 1000 такая функция даст 400. Это означает, что есть 400 простых чисел в интервале от 0 до 1000.“ volodya, Вы немного не правы. Ф-ция Эйлера возвращает количество натуральных чисел, меньших ее аргумента и взаимно простых с ним. Действительно, EulerPhi[1000] == 400. Но количество простых чисел в интервале [1, 1000] равно pi(1000) == 168. P.S. Спасибо за приведенную ссылку. P.P.S. Очень рекомендую пакет Mathematica - IMHO весьма полезный для реверсинга инструмент ;-) |
|
|
Дата: Авг 14, 2004 21:38:23 · Поправил: volodya volodya, Вы немного не правы. Ох, да :( Спасибо, что поправил! "Функция Эйлера phi( a ) есть количество чисел из ряда 0, 1, 2,..., a - 1, взаимно простых с a ." Да, Mathematica, штука, конечно, хорошая, но мне нравится Maple :) |
|
|
Дата: Авг 15, 2004 10:22:36 [offtop] volodya Maple? Честно говоря, я боюсь её использовать. В ней стоолько багов, причём таких откровенных, что зачастую приходится пересчитывать в других пакетах… |
|
|
Дата: Авг 17, 2004 22:37:48 · Поправил: volodya Ладно, продолжаем. Сегодня у нас в гостях метод Полларда (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 = 2kq: http://www.wasm.ru/forum/index.php?action=vthread&forum=17&topic=6841 Теперь мы готовы. Код:
static unsigned powmod(unsigned a, unsigned k, unsigned n)
{
unsigned b=1;
while (k)
{
if (!(k%2))
{
k /= 2;
a = (a*a)%n;
}
else
{
k--;
b = (b*a)%n;
}
}
return b;
}
/* Donald Knuth, vol 2., algorithm 4.5.4P - probability test for prime*/
static bool isprime(unsigned n)
{
srand(time(0));
unsigned x = rand() % (n - 1);
unsigned k, q;
do_N(n, k, q);
unsigned j = 0;
unsigned y = powmod(x, q, n);
while(j < k)
{
if(((!j) && (y == 1)) || (y == n - 1))
return true;
if((j > 0) && (y == 1))
return false;
j++;
y = (y*y)%n;
}
return false;
}
static unsigned pollard_rho(unsigned n)
{
if(isprime(n))
return 0;
unsigned factor = 1, x = 5, y = 2;
while(gcd(x-y, n) == 1)
{
x = (x*x + 1)%n;
x = (x*x + 1)%n;
y = (y*y + 1)%n;
}
factor = gcd(x-y, n);
if(factor < n)
return factor;
return 0;
}
void main()
{
cout << pollard_rho(561);
}
|
|
|
Дата: Авг 17, 2004 23:26:04 |
|
|
Дата: Авг 18, 2004 00:20:39 Метод разложения на множители Ферма лучше всего описан тут: http://www.wasm.ru/forum/index.php?action=vthread&forum=18&topic=6661 Эту книгу стоит купить. В виде алгоритма, что предлагает Кнут, это может выглядеть так: #include <iostream>
#include <math.h>
using namespace std;
/* 4.5.4C */
void fermat(unsigned n)
{
double x = 2*floor(sqrt((double)n)) + 1;
unsigned y = 1;
double r = pow((floor(sqrt((double)n))), 2) - n;
for(;;)
{
if(!r)
{
cout << (x-y)/2 << endl;
return;
}
r += x;
x += 2;
do
{
r -= y;
y += 2;
}while(r > 0);
}
}
void main()
{
fermat(377);
} |
|
|
Дата: Авг 20, 2004 23:49:09 В этом разделе ранее был приведен популярный вероятностный тест Рабина-Миллера на простоту числа (был взят из Кнута). Не так давно тройка веселых индусов открыла куда более крутой детерминистический тест. Тест очень быстр. Хм. Но и не слишком-то прост :) Реализацию теста можно, например, найти в Miracl. Сама статья находится по адресу: http://www.cse.iitk.ac.in/news/primality.html. Чтобы написать этот алгоритм на С/С++ самому, вам потребуется большинство из процедур, приведенных в этом топике выше, а также процедура проверки числа на степень. Процедуру можно найти тут: http://www.wasm.ru/forum/index.php?action=vthread&forum=17&topic=6888 |
|
|
Дата: Авг 21, 2004 14:13:25 |
|
|
Дата: Авг 21, 2004 14:55:23 Loger „http://www.cse.iitk.ac.in/news/primality.html Error 404“ Как ни странно, после поиска в Гугле и закачки старого варианта алгоритма http://www.cse.iitk.ac.in/primality.pdf страничка опять появилась. У них там все через базу данных организовано - видимо свой кеш почистили.... |
|
|
Дата: Ноя 1, 2004 11:43:59 В ссылке от volodya лишняя точка (после html), посему у Loger и Error 404 появлялось :)) |