Reputation: 6692
I am trying to implement Pollard Rho integer factorization in C/C++.Google gives me a Java implementation of the problem here.
I don't know Java that well,so what I came up with this.My implemenation in C++ works for most cases but not in few like the one "9999", I used there.
I am aware that C++ didn't have Biginteger class so I can't have the full functionality as it gives in JAVA but I want to factorize 15 digits numbers that's sufficient for unsigned long long
Please point out what wrong in my implementation.
Upvotes: 3
Views: 2675
Reputation: 339
You can try this ~100 lines C implementation of Pollard Rho :
There is some helpers :
#include <stdlib.h>
#include <stdint.h>
typedef uint_fast64_t num ;
static inline num mul_mod(num a, num b, const num mod) {
// Return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
num res = 0, tmp;
for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (tmp = b) >= mod - b ? tmp -= mod : 0, b += tmp);
return res % mod;
}
static inline num square_root(num n) {
// Return the number that was multiplied by itself to reach N.
num a = 0, b, c;
for (b = 1ULL << 62; b; c = a + b, n -= c &= -(n >= c), a = (a >> 1) | (c & b), b >>= 2);
// Variable n contains the remainder.
return a;
}
There is a required prime checker :
static int is_prime(const num n, size_t iterations) {
// Perform a Miller-Rabin (strong probable prime) test.
num a = 0, b, c, d, e, f; int h, i;
if ((n == 1) == (n & 1)) return n == 2;
for (b = c = n - 1, h = 0; !(b & 1); b >>= 1, ++h);
for (; iterations--;) {
for (size_t g = 0; g < sizeof(a); ((char*)&a)[g++] = rand()); // random input.
do for (d = e = 1 + a % c, f = n; (d %= f) && (f %= d););
while (d > 1 && f > 1);
for (d = f = 1; f <= b; f <<= 1);
for (; f >>= 1; d = mul_mod(d, d, n), f & b && (d = mul_mod(e, d, n)));
if (d == 1) continue;
for (i = h; i-- && d != c; d = mul_mod(d, d, n));
if (d != c) return 0;
}
return 1;
}
There is two factorization functions :
static inline num factor_worker_2(const num n, size_t limit) {
// Perform a Pollard's Rho probabilistic test.
size_t a = -1, b = 2;
num c, d = 1 + rand(), e = 1, f = 1;
for (c = d %= n; f == 1 && --limit; d = c, b <<= 1, a = -1) {
for (; f |= e, f == 1 && ++a != b;) {
c = mul_mod(c, c, n);
for (++c, c *= c != n, e = n, f = c > d ? c - d : d - c; (f %= e) && (e %= f););
}
}
return f;
}
static inline num factor_worker_1(const num n) {
// Perform a trial divisions test on N.
static const num list[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 1};
size_t i;
for (i = -1; n % list[++i];);
return list[i];
}
There is a factorizarion manager :
num factor(const num n) {
// Basic factorization manager, detect primes, perfect squares, execute workers.
num res;
switch (n) {
case 0: case 1: case 2: case 3:
res = 1; break;
default:
res = factor_worker_1(n);
if (res == 1 && !is_prime(n, 20)) {
res = square_root(n);
if (res * res != n)
for(;res = factor_worker_2(n, -1), res == 1 || res == n;);
}
}
return res;
}
There is a main :
#include <assert.h>
#include <stdio.h>
int main(void) {
num N;
N = 951818131364430049;
printf("factor is %zu\n", factor(N));
}
To try it :
// You can put it into a main.c file then compile + execute :
// gcc -O3 -std=c99 -Wall -pedantic main.c ; ./a.out ;
Here is the source for more informations, Thank You.
Upvotes: 0
Reputation: 6334
I've spotted one difference: the Java code assigns c
and x
to be new BigInteger(N.bitLength(), random)
, whereas the C++ code uses rand() % N
, which is a smaller random range. For the value 9999, the binary is 10011100001111, so the Java code will give c
and x
a maximum value of 16383.
Upvotes: 2
Reputation: 30561
The problem's right here:
#define abs(x) (x>0)?(x):(-x)
You're missing some parentheses in your abs
macro. Try:
#define abs(x) ((x)>0 ? (x) : -(x))
instead. (Consider what happens when abs(x-xx)
is expanded in the case x-xx <= 0
.)
Also, why does your gcd function return an int rather than a BigInteger?
You should also be aware that (assuming unsigned long long is a 64-bit integer type) this code won't work correctly for N
larger than 2**32
: if x
(or xx
) is greater than or equal to 2**32
then x*x
will wrap modulo 2**64
, giving you the wrong value for x*x % N
.
Upvotes: 10