Skip to content

Commit e796f16

Browse files
committed
pollard refactor
1 parent 8b936c5 commit e796f16

File tree

1 file changed

+18
-15
lines changed

1 file changed

+18
-15
lines changed

content/english/hpc/algorithms/factorization.md

Lines changed: 18 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -271,12 +271,13 @@ u64 diff(u64 a, u64 b) {
271271
return a > b ? a - b : b - a;
272272
}
273273

274+
const u64 SEED = 42;
275+
274276
u64 find_factor(u64 n) {
275-
u64 x = x0, y = x0, g = 1;
277+
u64 x = SEED, y = SEED, g = 1;
276278
while (g == 1) {
277-
x = f(x, a, n);
278-
y = f(y, a, n);
279-
y = f(y, a, n);
279+
x = f(f(x, n), n); // advance x twice
280+
y = f(y, n); // advance y once
280281
g = gcd(diff(x, y));
281282
}
282283
return g;
@@ -290,13 +291,13 @@ While it processes 25k 30-bit numbers — almost 15 times slower than the fastes
290291
Floyd's cycle-finding algorithm has a problem in that it does more iterator increments than necessary. One way to solve it is to memorize the values that the faster iterator visits and compute the gcd using the difference of $x_i$ and $x_{\lfloor i / 2 \rfloor}$, but it can also be done without extra memory using this trick:
291292
292293
```c++
293-
u64 find_factor(u64 n, u64 x0 = 2, u64 a = 1) {
294-
u64 x = x0, y = x0;
294+
u64 find_factor(u64 n) {
295+
u64 x = SEED;
295296
296297
for (int l = 256; l < (1 << 20); l *= 2) {
297-
x = y;
298+
u64 y = x;
298299
for (int i = 0; i < l; i++) {
299-
y = f(y, a, n);
300+
x = f(x, n);
300301
if (u64 g = gcd(diff(x, y), n); g != 1)
301302
return g;
302303
}
@@ -313,14 +314,14 @@ We can remove the logarithm from the asymptotic using the fact that if one of $a
313314
```c++
314315
const int M = 1024;
315316

316-
u64 find_factor(u64 n, u64 x0 = 2, u64 a = 1) {
317-
u64 x = x0, y = x0, p = 1;
317+
u64 find_factor(u64 n) {
318+
u64 x = SEED;
318319

319320
for (int l = M; l < (1 << 20); l *= 2) {
320-
x = y;
321+
u64 y = x, p = 1;
321322
for (int i = 0; i < l; i += M) {
322323
for (int j = 0; j < M; j++) {
323-
y = f(y, a, n);
324+
y = f(y, n);
324325
p = (u128) p * diff(x, y) % n;
325326
}
326327
if (u64 g = gcd(p, n); g != 1)
@@ -340,6 +341,8 @@ The next step is to actually apply [Montgomery Multiplication](/hpc/number-theor
340341
341342
This is exactly the type of problem when we need specific knowledge, because we have 64-bit modulo by not-compile-constants, and compiler can't really do much to optimize it.
342343
344+
We do not need to convert numbers out of Montgomery representation before computing the GCD.
345+
343346
```c++
344347
struct Montgomery {
345348
u64 n, nr;
@@ -369,13 +372,13 @@ const int M = 1024;
369372
370373
u64 find_factor(u64 n, u64 x0 = 2, u64 a = 1) {
371374
Montgomery m(n);
372-
u64 y = x0;
375+
u64 x = SEED;
373376
374377
for (int l = M; l < (1 << 20); l *= 2) {
375-
u64 x = y, p = 1;
378+
u64 y = x, p = 1;
376379
for (int i = 0; i < l; i += M) {
377380
for (int j = 0; j < M; j++) {
378-
y = f(y, a, m);
381+
x = f(x, m);
379382
p = m.multiply(p, diff(x, y));
380383
}
381384
if (u64 g = gcd(p, n); g != 1)

0 commit comments

Comments
 (0)