# HG changeset patch # User Arun Giridhar # Date 1667662316 14400 # Node ID 0bd5ea4bd31c8e174f52c9748f262d6a69218538 # Parent 7e60506a54281b12887d15badb349cec1ff1870e __isprimelarge.cc__: Minor performance tweaks. __isprimelarge.cc__: Change starting value and factorization polynomial from `x^2 - c` to `x^2 + c`. Use error_unless for rare theoretical error condition that should not happen in practice. Minor rearrangement of conditional operator. diff -r 7e60506a5428 -r 0bd5ea4bd31c libinterp/corefcn/__isprimelarge__.cc --- a/libinterp/corefcn/__isprimelarge__.cc Sat Nov 05 13:48:24 2022 +0100 +++ b/libinterp/corefcn/__isprimelarge__.cc Sat Nov 05 11:31:56 2022 -0400 @@ -216,7 +216,7 @@ pollardrho (uint64_t n, uint64_t c = 1ULL) { uint64_t i = 1ULL, j = 2ULL; // cycle index values - uint64_t x = (c+2) % n; // can also be rand () % n + uint64_t x = (c+1) % n; // can also be rand () % n uint64_t y = x; // other value in the chain uint64_t g = 0; // GCD @@ -224,12 +224,13 @@ { i++; - // Calculate x = mod (x^2 - c, n) without overflow. - x = safemultiply (x, x, n); - x = (x >= c) ? (x - c) : (x + n - c); + // Calculate x = mod (x^2 + c, n) without overflow. + x = safemultiply (x, x, n) + c; + if (x >= n) + x -= n; // Calculate GCD (abs (x-y), n). - g = (x == y) ? 0 : (x > y) ? localgcd (x - y, n) : localgcd (y - x, n); + g = (x > y) ? localgcd (x - y, n) : (x < y) ? localgcd (y - x, n) : 0; if (i == j) // cycle detected ==> double j { @@ -242,10 +243,8 @@ if (g > 1ULL) // found GCD ==> exit loop properly { - if (n % g) // nonzero remainder ==> not a factor - return 0; // rare but theoretically possible case. - else // normal case - return g; // GCD is a divisor of n by definition. + error_unless (n % g == 0); // theoretical possibility of GCD error + return g; } } }