Monday, 8 June 2009

Finalized(?!) (GMP-)WE algorithm code

/* Factoring general integers with WE method using random base(s).
Author: James Wanless (c) 2000-9
*/

#include
#include
#include

#include
#include
#include

gmp_randclass r (gmp_randinit_default);

mpz_class rand2(long bits)
/* returns random bits-bit integer */
{
return r.get_z_bits(bits);
}

int Rabin_Miller(mpz_class n)
/* given an integer n >= 3 returns 0 composite, 1 probably prime */
{
return mpz_probab_prime_p(n.get_mpz_t(), 10);
}

mpz_class gcd(mpz_class a, mpz_class b)
/* returns greatest common divisor of a and b */
{
mpz_t temp;
mpz_init(temp);
mpz_gcd(temp, a.get_mpz_t(), b.get_mpz_t());
mpz_class temp_class(temp);
mpz_clear(temp);
return temp_class;
}

mpz_class exp_mod(mpz_class x, mpz_class b, mpz_class n)
/* returns x ^ b mod n */
{
mpz_t temp;
mpz_init (temp);
mpz_powm(temp, x.get_mpz_t(), b.get_mpz_t(), n.get_mpz_t());
mpz_class temp_class(temp);
mpz_clear (temp);
return temp_class;
}

bool Wanless (mpz_class n)
{
mpz_class basestested = 0;
mpz_class T = 1;
mpz_class b = 1;
mpz_class A = 2;
mpz_class prodA = 2;
mpz_class wanless = 2;
long nbits = 1;

if (Rabin_Miller(n)) {
cout << n << "\n";
fflush(stdout);
return true;
}

if (n < 2)
return false;

if (n % 2 == 0) {
cout << "2\n";
fflush(stdout);
return Wanless(n / 2);
}


while (wanless < n) {
wanless *= 2;
nbits ++;
}

while (T == 1 || T == n) {
basestested++;
cout << "base#" << basestested << "\r";
fflush(stdout);

prodA = 1;
for (long i = 0; i < nbits; i++) {
A = rand2(nbits);
prodA *= A;
prodA = prodA % n;
}


b = exp_mod(prodA, wanless, n);
T = gcd(n, exp_mod(b, n, n) - b);
}

if (T > 1 && T < n) {
cout << "\n" << T << "\t[prodA=" << prodA << "]\n";
fflush(stdout);
Wanless(n / T);
}


return (false);

}

int main(void)
{
mpz_class N;

for (;;) {
cout << "number to be tested or 0 to quit:\n";
cin >> N;
if (N == 0) break;
Wanless(N);
}
return 0;
}

Thursday, 5 February 2009

Monday, 2 February 2009

CRTfac3

I wonder whether (quite simply) this was what I was aiming at...:
(copyright 2009-02-02 JGW, again, just-in-case :)

import java.math.BigInteger;
import java.util.Random;
import java.util.Date;

public class CRTfac3 {
static BigInteger zero = new BigInteger("0");
static BigInteger one = new BigInteger("1");
static BigInteger two = new BigInteger("2");
static BigInteger thousand = new BigInteger("1000");

static BigInteger n = new BigInteger("0");
static String s;

static int olds1 = 0;
static int s1 = 0;

public static void main (String args[])
throws java.io.IOException {

CRTfac3 CRTfac3inst = new CRTfac3();

char c;
String sInput;

Date d;
long starttime;
long finishtime;
long duration;

while (true) {
StringBuffer sbInput = new StringBuffer("");

while ((c = (char)System.in.read()) != '\n' && c != '\r')
sbInput.append(c);
System.in.read();
sInput = sbInput.toString().trim();


if (sInput.charAt(0) == 'f' || sInput.charAt(0) == 'F') {
s = sInput.substring(1).trim();

s1 = 0;
olds1 = 0;
n = CRTfac3inst.eval(s);

System.out.println('[' + n.toString() + ']');
}
else {
n = new BigInteger(sInput);
}


d = new Date();
starttime = d.getTime();

CRTfac3inst.factorize(n);

d = new Date();
finishtime = d.getTime();

duration = (finishtime-starttime)/1000;

System.out.println("duration: " + duration + " seconds");

System.out.println();
}
}

public boolean factorize(BigInteger n) {
boolean prime = false;
BigInteger a = new BigInteger("1");
BigInteger b = new BigInteger("1");
BigInteger T = new BigInteger("1");

if (n.isProbablePrime(1000)) {
prime = true;
System.out.println(n);
return prime;
}

while (T.compareTo(one) == 0 || T.compareTo(n) == 0) {
T = n.gcd(b);
if (a.mod(thousand).compareTo(zero) == 0)
System.out.print("a = " + a.toString() + '\r');
a = a.add(one);
b = b.multiply(a).mod(n);
}


if (T.compareTo(one) > 0 && T.compareTo(n) < 0) {
factorize(T);
factorize(n.divide(T));
}
else {
prime = false;
System.out.println("composite" +'\t' + n);
}

return prime;

}


public BigInteger evalRand(char op, BigInteger oldn) {
BigInteger n = new BigInteger("1");


switch (op) {
case 'r':
case 'R':
Random r = new Random();

n = new BigInteger(oldn.intValue(), r);
break;

default:
n = oldn;
break;
}
return n;
}

public BigInteger evalFact(BigInteger oldn, char op) {
BigInteger n = new BigInteger("1");
BigInteger i = new BigInteger("1");
BigInteger j = new BigInteger("1");
boolean prime = true;


switch (op) {
case '!':
for (i = one; i.compareTo(oldn) <= 0; i = i.add(one))
n = n.multiply(i);
break;

case '#':
for (i = one; i.compareTo(oldn) <= 0; i = i.add(one)) {
prime = true;
for (j = two; (prime == true) && (j.multiply(j).compareTo(i) <= 0); j = j.add(one))
if (i.remainder(j).compareTo(zero) == 0)
prime = false;
if (prime == true)
n = n.multiply(i);
}
break;

default:
n = oldn;
break;
}
return n;
}

public BigInteger evalPower(BigInteger oldn, BigInteger n1, char op) {
BigInteger n = new BigInteger("0");

switch (op) {
case '^':
n = oldn.pow(n1.intValue());
break;
default:
n = n1;
break;
}
return n;
}

public BigInteger evalProduct(BigInteger oldn, BigInteger n1, char op) {
BigInteger n = new BigInteger("0");

switch (op) {
case '*':
n = oldn.multiply(n1);
break;
case '/':
n = oldn.divide(n1);
break;
case '%':
n = oldn.remainder(n1);
break;
default:
n = n1;
break;
}
return n;
}

public BigInteger evalSum(BigInteger oldn, BigInteger n1, char op) {
BigInteger n = new BigInteger("0");

switch (op) {
case '+':
n = oldn.add(n1);
break;
case '-':
n = oldn.subtract(n1);
break;
default:
n = n1;
break;
}
return n;
}

public BigInteger eval(String s) {
BigInteger oldn0 = new BigInteger("0");
BigInteger oldn1 = new BigInteger("0");
BigInteger oldn2 = new BigInteger("0");
BigInteger n = new BigInteger("0");

char oldop0 = 0;
char oldop1 = 0;
char oldop2 = 0;
char op = 0;

while (s1 < s.length()) {
switch (s.charAt(s1)) {
case '(':
case '[':
case '{':
s1++;
n = eval(s);
break;

case '0':
case '1':
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
case '8':
case '9':
n = readNum(s);
break;

default:
break;
}
if (s1 < s.length()) {
switch (s.charAt(s1)) {

case ')':
case ']':
case '}':
case '!':
case '#':
case 'r':
case 'R':
case '^':
case '*':
case '/':
case '%':
case '+':
case '-':
op = s.charAt(s1);
s1++;
break;

default:
break;
}
}
else
op = 0;


switch (op) {
case 0:
case ')':
case ']':
case '}':
n = evalPower(oldn2, n, oldop2);
n = evalProduct(oldn1, n, oldop1);
n = evalSum(oldn0, n, oldop0);
return n;

case '!':
case '#':
n = evalFact(n, op);
break;

case 'r':
case 'R':
n = readNum(s);
n = evalRand(op, n);
break;

case '^':
n = evalPower(oldn2, n, oldop2);
oldn2 = n;
oldop2 = op;
break;

case '*':
case '/':
case '%':
n = evalPower(oldn2, n, oldop2);
oldop2 = 0;
n = evalProduct(oldn1, n, oldop1);
oldn1 = n;
oldop1 = op;
break;

case '+':
case '-':
n = evalPower(oldn2, n, oldop2);
oldop2 = 0;
n = evalProduct(oldn1, n, oldop1);
oldop1 = 0;
n = evalSum(oldn0, n, oldop0);
oldn0 = n;
oldop0 = op;
break;
default:
break;
}
}
return n;
}

public BigInteger readNum(String s) {
olds1 = s1;
while (s1 < s.length() && Character.isDigit(s.charAt(s1)))
s1++;
n = new BigInteger(s.substring(olds1, s1));

return n;
}


}

bear...@gmail.com wrote:
> Is it possible to use CRT to factor integers - or is this in effect
> what modern sieving algorithms already do anyway?
> J