Алгоритм Брента-Полларда с бесконечным циклом
У меня есть следующая функция. Я получил его с двух сайтов и попытался адаптировать его к своему, но он не очень хорошо работал.
Когда я проверяю unsigned long int max - 2 или, но это как число 4294967293 он помещает следующий код в бесконечный цикл, где ys будет продолжать возвращаться к тому же значению.
Я очень новичок в алгоритмах факторизации и пошагово помогаю понять, почему я получаю бесконечный цикл, было бы здорово.
Следующий код — просто моя функция «RHO». У меня есть другая функция под названием gdc который такой же, как и любой другой gdc рекурсивная функция там.
unsigned long int rho(unsigned long int n)
{
if (n == 1) return n;
if (n % 2 == 0) return 2;
unsigned long int y = rand() % n;
unsigned long int x;
unsigned long long int ys = y;
unsigned long int c;
do
c = rand() % n;
while (c == 0 || c == n - 2);
unsigned long int m = 1000;
unsigned long int d = 1;
unsigned long int q = 1;
unsigned long int r = 1;
while (d == 1)
{
x = y;
for (int i = 0; i < r; i++)
{
y = y * y % n;
y += c;
if (y < c)
y += (std::numeric_limits<unsigned long>::max() - n) + 1;
y %= n;
}
int j = 0;
while (j < r && d == 1)
{
ys = y;
for (int i = 0; i < m && i < (r-j); i++)
{
y = y * y % n;
y += c;
if (y < c)
y += (std::numeric_limits<unsigned long>::max() - n) + 1;
y %= n;
q *= ((x>y) ? x - y : y - x) % n;
}
d = gcd(q, n);
j += m;
}
r *= 2;
}
if (d == n)
{
do
{
ys = ys * ys % n;
std::cout << ys << std::endl;
ys += c;
if (ys < c)
ys += (std::numeric_limits<unsigned long>::max() - n) + 1;
ys %= n;
d = gcd( ((x>ys) ? x - ys : ys - x) , n);
} while (d == 1);
}
return d;
}
Примеры, из которых я адаптировал свой код:
РЕДАКТИРОВАТЬ
Я сделал, как предложил Амд, исправил мой код и переместил повторяющиеся строки в вспомогательные функции. Однако я все еще получаю бесконечный цикл для d==n часть около дна. По какой-то причине f(ys) в конечном итоге возвращает то же самое, что и ранее, поэтому он продолжает циклически проходить через ряд значений.
uint64_t rho(uint64_t n)
{
if (n == 1) return n;
if (n % 2 == 0) return 2;
uint64_t y = rand() % n;
uint64_t x;
unsigned long long int ys = y;
uint64_t c;
do c = rand() % n; while (c == 0 || c == n - 2);
uint64_t m = 1000;
uint64_t d = 1;
uint64_t q = 1;
uint64_t r = 1;
do
{
x = y;
for (int i = 0; i <= r; i++)
y = f(y, c, n);
int j = 0;
do
{
ys = y;
for (int i = 0; i <= min(m, r-j); i++)
{
y = f(y, c, n);
q *= (abs(x,y) % n);
}
d = gcd(q, n);
j += m;
} while (j < r && d == 1);
r *= 2;
} while (d == 1);
if (d == n)
{
do
{
ys = f(ys, c, n);
d = gcd(abs(x, ys), n);
} while (d == 1);
}
return d;
}
Решение
Это должно всегда заканчиваться. К тому времени, когда вы доберетесь до этой точки в алгоритме, x будет в цикле (так как y началось в x и прошел весь путь, чтобы обнаружить цикл, используя обнаружение цикла Брента). Значение ys началось в y который начался в x, так что это также будет продолжаться вокруг цикла и в конечном итоге встретиться с x снова (см. Обнаружение цикла Флойд или Черепаха и Заяц). На этом этапе вы будете иметь gcd(ys,x)==xи этот последний цикл завершится.
В опубликованной реализации есть несколько ошибок, которые, по моему мнению, могут вызывать проблему:
x = y;
for (int i = 0; i < r; i++) // should be strictly less than
...
ys = y;
for (int i = 0; i < min(m, r-j); i++) // again, strictly less than
{
y = f(y, c, n);
q = (q*abs(x,y)) % n; // needs "mod" operator AFTER multiplication
}
...
Вы также можете заменить инициализацию c с
uint64_t c = (rand() % (n-3))+1
если вы хотите число в диапазоне [1, n-3].
Вот оригинальная статья Ричарда П. Брента: Усовершенствованный алгоритм факторизации Монте-Карло
Другие решения
Редактировать:
это работает для меня, а не в ловушке внутри бесконечного цикла:
#include<iostream>
#include<stdint.h>
#define min(a,b) (a<b?a:b)
#define abs(x,y) (x > y? x - y : y - x)
uint64_t gcd(uint64_t m, uint64_t n) {
while (true) {
int r = m % n;
if (r == 0) {
return n;
}
m = n;
n = r;
}
}
uint64_t f(uint64_t y, uint64_t c, uint64_t n) {
y = (y * y) % n;
y += c;
if (y < c)
y += (std::numeric_limits<uint32_t>::max() - n) + 1;
y %= n;
return y;
}
uint64_t rho(uint64_t n)
{
if (n == 1) return n;
if (n % 2 == 0) return 2;
uint64_t y = rand() % n;
uint64_t x;
uint64_t ys = y;
uint64_t c;
do c = rand() % n; while (c == 0 || c == n - 2);
uint64_t m = 1000;
uint64_t d = 1;
uint64_t q = 1;
uint64_t r = 1;
do
{
x = y;
for (int i = 0; i <= r; i++)
y = f(y, c, n);
int j = 0;
do
{
ys = y;
for (int i = 0; i <= min(m, r - j); i++)
{
y = f(y, c, n);
q *= (abs(x, y) % n);
}
d = gcd(q, n);
j += m;
} while (j < r && d == 1);
r *= 2;
} while (d == 1);
if (d == n)
{
do
{
ys = f(ys, c, n);
d = gcd(abs(x, ys), n);
} while (d == 1);
}
return d;
}
int main() {
std::cout << rho(std::numeric_limits<uint32_t>::max() - 2) << "\n"; //9241
}
Старый:
в соответствии с улучшенным алгоритмом факторизации Монте-Карло: http://wwwmaths.anu.edu.au/~brent/pd/rpb051i.pdf
С.: 182-183
ошибка: ys = x * x% n;
правильно: ys = ys * ys% n; // ys = f (ys)
пожалуйста используйте uint32_t или uint64_t … из stdint.h для хорошего и чистого стиля кодирования:
#include<iostream>
#include<stdint.h>
int gcd(int m, int n) {
while (true) {
int r = m % n;
if (r == 0) {
return n;
}
m = n;
n = r;
}
}
#define min(a,b) (a<b?a:b)
#define abs(x,y) (x > y? x - y : y - x)
uint32_t f(uint32_t y, uint32_t c, uint32_t n) {
y = (y * y) % n;
y += c;
if (y < c)
y += (std::numeric_limits<uint32_t>::max() - n) + 1;
y %= n;
return y;
}
//http://wwwmaths.anu.edu.au/~brent/pd/rpb051i.pdf
//pp:182-183
uint32_t rho(uint32_t n) {
if (n == 1) return n;
if (n % 2 == 0) return 2;
uint32_t y = rand() % n; // y = x0
uint32_t x;
uint64_t ys = y;
uint32_t c;
do { c = rand() % n; } while (c == 0 || c == n - 2);
uint32_t m = 1000;
uint32_t d = 1;
uint32_t q = 1;
uint32_t r = 1;
do
{
x = y;
for (int i = 1; i <= r; i++) y = f(y, c, n);
int j = 0;
do
{
ys = y;
for (int i = 1; i <= min(m, r - j); i++)
{
y = f(y, c, n);
q = q * abs(x, y) % n;
}
d = gcd(q, n);
j += m;
} while (j < r && d == 1);
r *= 2;
} while (d == 1);
if (d == n)
{
do
{
ys = f(ys, c, n);//not: ys = f(x,c,n);
d = gcd(abs(x, ys), n);
} while (d == 1);
}
return d;
}