You can try this simple C99 implementation of Pollard Rho algorithm :
// Integer factorization in C language.
// Decompose a composite number into a product of smaller integers.
unsigned long long pollard_rho(const unsigned long long N) {
// Require : N is a composite number, not a square.
// Ensure : you already performed trial division up to 23.
// Option : change the timeout, change the rand function.
static const int timeout = 18;
static unsigned long long rand_val = 2994439072U;
rand_val = (rand_val * 1025416097U + 286824428U) % 4294967291LLU;
unsigned long long gcd = 1, a, b, c, i = 0, j = 1, x = 1, y = 1 + rand_val % (N - 1);
for (; gcd == 1; ++i) {
if (i == j) {
if (j >> timeout)
break;
j <<= 1;
x = y; // "x" takes the previous value of "y" when "i" is a power of 2.
}
a = y, b = y; // computes y = f(y)
for (y = 0; a; a & 1 ? b >= N - y ? y -= N : 0, y += b : 0, a >>= 1, (c = b) >= N - b ? c -= N : 0, b += c);
y = (1 + y) % N; // function f performed f(y) = (y * y + 1) % N
for (a = y > x ? y - x : x - y, b = N; (a %= b) && (b %= a););
gcd = a | b; // the GCD(abs(y - x), N) was computed
// it continues until "gcd" is a non-trivial factor of N.
}
return gcd;
}
- Usually you performed some trial division before calling the algorithm
- The algorithm isn't designed to receive a prime number as input
- Two consecutive calls may not result in the same answer
- Alternately, there is a pure C quadratic sieve which factors numbers from 0 to 300-bit.
If in doubt about the primality of N you can use a C99 primality checker :
typedef unsigned long long int ulong;
ulong mul_mod(ulong a, ulong b, const ulong mod) {
ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
return res % mod;
}
ulong pow_mod(ulong n, ulong exp, const ulong mod) {
ulong res = 1; // return (n ^ exp) % mod
for (n %= mod; exp; exp & 1 ? res = mul_mod(res, n, mod) : 0, n = mul_mod(n, n, mod), exp >>= 1);
return res;
}
int is_prime(ulong N) {
// Perform a Miller-Rabin test, it should be a deterministic version.
const ulong n_primes = 9, primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
for (ulong i = 0; i < n_primes; ++i)
if (N % primes[i] == 0) return N == primes[i];
if (N < primes[n_primes - 1]) return 0;
int primality = 1, a = 0;
ulong b;
for (b = N - 1; ~b & 1; b >>= 1, ++a);
for (ulong i = 0; i < n_primes && primality; ++i) {
ulong c = pow_mod(primes[i], b, N);
if (c != 1) {
for (int j = a; j-- && (primality = c + 1 != N);)
c = mul_mod(c, c, N);
primality = !primality;
}
}
return primality;
}
To try it there is a factor function :
// return the number that was multiplied by itself to reach N.
ulong square_root(const ulong num) {
ulong res = 0, rem = num, a, b;
for (a = 1LLU << 62 ; a; a >>= 2) {
b = res + a;
res >>= 1;
if (rem >= b)
rem -= b, res += a;
}
return res;
}
ulong factor(ulong num){
const ulong root = square_root(num);
if (root * root == num) return root ;
const ulong n_primes = 9, primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
for (ulong i = 0; i < n_primes && primes[i] <= root; ++i)
if (num % primes[i] == 0) return primes[i];
if (is_prime(num))
return 1 ;
return pollard_rho(num);
}
Which is completed by the main function :
#include <assert.h>
int main(void){
for(ulong i = 2; i < 63; ++i){
ulong f = factor(i);
assert(f <= 1 || f >= i ? is_prime(i) : i % f == 0);
ulong j = (1LLU << i) - 1 ;
f = factor(j);
assert(f <= 1 || f >= j ? is_prime(j) : j % f == 0);
j = 1 | pow_mod((ulong) &main, i, -5);
f = factor(j);
assert(f <= 1 || f >= j ? is_prime(j) : j % f == 0);
}
}