Here is the code for a factor-like program that uses Pollard's Rho
algorithm. It doesn't do any error-checking or anything fancy, but its
output is identical to that of factor, for testing purposes.
It is slower for small inputs, so we should probably fall back to the
wheel method for these. However, I haven't yet found an input for which
it takes more than a second or so.
--Trevor
"Mathematics is like checkers in being suitable for the young, not too
difficult, amusing, and without peril to the state." --Plato
On Mon, 5 Jan 2004, Jim Meyering wrote:
> > Are there any plans to implement faster algorithms in factor? For 64-bit
> > integers Pollard's rho method would be a good choice. I have an
> > implementation that is hundreds of times faster than factor for
> > large inputs. Is anyone interested in this?
>
> Sounds interesting to me.
>
>
// Copyright (C) 2004 Trevor Wilson
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>
#include <time.h>
#include <assert.h>
#define RMTRIALS 8
#define TRIM(a, n) if (a > n) a -= n;
typedef unsigned long long int num;
typedef struct {
num factor[64];
int len;
} factorization;
const factorization FACT_NULL = {.len = 0};
void factor(num n, factorization *fac);
num rho(num n);
bool isprime(num n);
bool witness(num n, num a);
num modexp(num a, num e, num n);
num modsq(num a, num n);
num modmul(num a, num b, num n);
num gcd(num a, num n);
int main(int argc, char *argv[])
{
int i;
for(i = 1; i < argc; i++) {
num n = strtoull(argv[i], NULL, 10);
int j;
factorization f = FACT_NULL;
factor(n, &f);
printf("%llu:", n);
for(j = 0; j < f.len; j++)
printf(" %llu", f.factor[j]);
printf("\n");
}
return 0;
}
void factor(num n, factorization *fac)
{
if (isprime(n)) {
int i;
for (i = fac->len++; i > 0; i--) {
if (fac->factor[i-1] <= n)
break;
fac->factor[i] = fac->factor[i-1];
}
fac->factor[i] = n;
} else {
num d = rho(n);
factor(d, fac);
factor(n/d, fac);
}
}
num rho(num n)
// Uses Pollard's Rho algorithm to find a nontrivial factor
{
num a, b, c;
if(~n & 1)
return 2;
redo:
a = b = 1;
c = rand() % (n-1) + 1;
for (;;) {
num d;
a = modsq(a,n) + c;
TRIM(a, n);
b = modsq(b,n) + c;
TRIM(b, n);
b = modsq(b,n) + c;
TRIM(b, n);
d = gcd(a < b ? b - a : a - b, n);
if(d == n)
goto redo;
if(d != 1)
return d;
}
}
bool isprime(num n)
{
int i;
assert(n > 1);
if (~n & 1)
return n == 2;
for(i = 0; i < RMTRIALS; i++)
if (witness(n, rand() % (n-2) + 2))
return false;
return true;
}
bool witness(num n, num a)
// The Rabin-Miller primality test
{
num s, c;
int r = 0;
assert(1 < a && a < n);
for (s = n - 1; ~s & 1; s >>= 1)
r++;
c = modexp(a, s, n);
if(c == 1)
return false;
while(r--) {
if(c == n-1)
return false;
c = modsq(c,n);
}
return true;
}
num modexp(num a, num e, num n)
// Binary modular exponentiation
{
num b = 1;
for (; e; e >>= 1) {
if (e & 1)
b = modmul(b,a,n);
a = modsq(a,n);
}
return b;
}
num modsq(num a, num n)
// This seems like a lame way to do 64-bit modular multiplication,
// but I can't think of anything better.
{
int i;
num a1 = a & 0xffffffff;
num a2 = a >> 32;
num c1 = a1*a1 % n;
num c2 = (a1*a2 % n) << 1;
num c3 = a2*a2 % n;
if (c2 >= n)
c2 -= n;
for(i = 0; i < 32; i++) {
c2 <<= 1;
TRIM(c2, n);
c3 <<= 1;
TRIM(c3, n);
c3 <<= 1;
TRIM(c3, n);
}
c1 += c2;
TRIM(c1, n);
c1 += c3;
TRIM(c1, n);
return c1;
}
num modmul(num a, num b, num n)
// Same comment as above
{
int i;
num a1 = a & 0xffffffff;
num b1 = b & 0xffffffff;
num a2 = a >> 32;
num b2 = b >> 32;
num c1 = a1*b1 % n;
num c2 = a1*b2 % n + a2*b1 % n;
num c3 = a2*b2 % n;
if (c2 >= n)
c2 -= n;
for(i = 0; i < 32; i++) {
c2 <<= 1;
TRIM(c2, n);
c3 <<= 1;
TRIM(c3, n);
c3 <<= 1;
TRIM(c3, n);
}
c1 += c2;
TRIM(c1, n);
c1 += c3;
TRIM(c1, n);
return c1;
}
num gcd(num a, num n)
// Euclid's Algorithm
{
while(a) {
num t = n % a;
n = a;
a = t;
}
return n;
}
_______________________________________________
Bug-coreutils mailing list
[EMAIL PROTECTED]
http://mail.gnu.org/mailman/listinfo/bug-coreutils
