History of factorization records by sieving methods

This graph is by Francois Morain. http://algo.inria.fr/seminars/sem00-01/morain.html
(image linked to in situ)
The relevant figure is towards the bottom of the article, in Section 4.

Penryn

It would appear that Intel's next-generation processors (Penryn) draw up to 40% less power than current CPUs, under full load. :))) Quite welcome I imagine, in the current climate (pun intended :( ).
http://www.tomshardware.com/2007/10/29/intel_penryn_4ghz_with_air_cooling/page14.html
This is no doubt due to moving to the smaller 45nm die resolution. Paul Otellini (CEO) has said that Penryn will ship 12 November 2007...

Coincidence?

Is it coincidence that Fermat's method works best with factors close to the square root of the input number, while NFS is particularly good at factoring semiprimes?

XYYXF Project milestone

Congratulations to the XYYXF project, who have just finished their original target of factoring all composites with x and y <=100, after 6 years! The search continues with x and y extended to 150...

http://xyyxf.at.tut.by/news.html#0

Daylight Savings Time

Now the clocks have gone back, all my ECMNet times are correct again ! :)

Msieve v1.29

New version of msieve (by Jason Papadopoulos) now available.
Download from:
http://www.boo.net/~jasonp/qs.html
Announcement at:
http://mersenneforum.org/showthread.php?p=117235#post117235

Fermat's Method

This 400-yr old method is the basis for many modern sieving factorization methods.
It relies on the discovery of an expression for the input number, as the difference between two (integral) squares.

http://en.wikipedia.org/wiki/Fermat%27s_factorization_method

Also, the pseudocode immediately below is taken from a PD document by Connelly Barnes of Oregon State University
http://oregonstate.edu/~barnesc/documents/factoring.pdf

function fermatFactor(N)
for x from ceil(sqrt(N)) to N
ySquared := x * x - N
if isSquare(ySquared) then
y := sqrt(ySquared)
s := (x - y)
t := (x + y)
if s <> 1 and s <> N then
return s, t
end if
end if
end for
end function

And here is its Java implementation (from superfac9) [with thanks to DT]:

BigInteger factorizefermat(BigInteger n) {
BigInteger a, bSq;
long iterations = 1L;

a = sqrt(n);
if(n.mod(a).compareTo(ZERO) == 0) return a;

while(iterations<10000000) {
bSq = a.multiply(a).subtract(n);
if(fastsquareQ(bSq))
return a.subtract(sqrt(bSq)); // bSq is a square, factorization found
a = a.add(ONE);
iterations++;
}

return n;
}

Mark Manasse

Here is a picture of Mark Manasse (reproduced by kind permission), who is notable for his involvement in the factorizations of F9, RSA-110 and RSA-120 in the early 1990's. In fact he helped develop the Number Field Sieve, which counted these numbers among its early successes.
http://www.std.org/~msm/common/f9paper.pdf
[Incidentally this paper, above, I see, also explains why WE works on prime-powers...basically Fermat's Little Theorem in combination w/ the Binomial Theorem]
http://en.wikipedia.org/wiki/RSA-110
http://en.wikipedia.org/wiki/RSA-120
http://www.std.org/~msm/common/nfspaper.pdf

"factorization"

Here's a bit of a Meta-post:
Regarding the spelling/usage of the word "factorization".
First off note that there is a synonym for "factorization", namely "factoring" - however since this also has a completely unrelated meaning (financial of some sort IIRC), I generally tend to use the longer form, unless it would sound particularly ugly.
Similarly the verb can just be "to factor", instead of "to factorize", however in this case the correct (mathematical) meaning is usually clearer, since it will be associated with an object (eg C127_113_36).
In addition (just to add to the confusion :) there is an alternative _spelling_ of "factorization" as "factorisation" (in British English).

Desktop

My desktop is pretty busy at the moment.
Here is a current screen grab.
The windows are (left-to-right):
1) Random-based WEP on (M+2)859433
2) WEP on F7
3) WEP on F8
4) WEP on (M+2)4253
5) msieve (NFS) on C127_113_36 (xyyxf)
6) ECM via ECMNet for xyyxf
7) GGNFS on rsa100
8) GGNFS on C127_113_36
All this on a 1.9GHz G5 (iMac called 'tiggatoo')!

GGNFS on PPC (update)

Basically, the problem of "not enough polynomials in mpqs", is fixed in the latest CVS code (0.77.1-20060722) available on Sourceforge (with a separate, preliminary, classical sieve stage).
This is the story of how I succeeded in obtaining, compiling and running it on my iMac G5 (with thanks to Mark Rodenkirch, who responded to my query on the yahoo GGNFS mailing list).
Step 1)
Obtain, compile and install a 64-bit GMP archive (PPC 970)
Step 2)
Download the source code of GGNFS via CVS from Sourceforge, thus:
"cvs -d:pserver:anonymous@ggnfs.cvs.sourceforge.net:/cvsroot/ggnfs
checkout -R branch_0"
The source-tree is also available for browsing at:
http://ggnfs.cvs.sourceforge.net/ggnfs/branch_0/src/
Step 3)
Necessary specifically (atm) for PPC [(probably) not for other architectures]:
tweak/hack the code to include its own implementation [supplied, but disabled] of getline() in if.c by,
a) inserting a suitable declaration for the missing getline() function at the top of file if.c
b) remming out the ifdef (and matching endif) of GGNFS_GNU_MISSING lower down, to actually activate the missing implementation of getline()
(I also investigated other methods of achieving same eg trying to define GGNFS_GNU_MISSING elsewhere to achieve the same effect, but this was the first method that worked.)
Step 4)
Compile with "make ppc_970"
Step 5)
Note that the default implementation of factLat.pl, by which factorizations are actually run, already has the correct path to the binaries, so does not need adjusting, unlike some previous versions.

Fermat number factor search

For the very latest on the search for Fermat number factors, checkout George Woltman's page at:
http://www.mersenne.org/ecmf.htm
Prime95/mprime includes an ECM factoring ability, easily controlled by a single ini file, and is particularly well-suited (and fast!) to looking for factors of larger Fermat (and Mersenneplustwo ;-) numbers. George Woltman even supplies a file of already-known factors of Fermat numbers, that mprime can read, to avoid the time spent re-discovering these (or you can start from scratch without reading-in that file, as a double-check that your set-up is finding factors successfully - this won't affect your ability to also find new factors)

Euclidean Algorithm

The Euclidean Algorithm is vital to many factorization algorithms.
It is a very fast method for finding the gcd (greatest common divisor) of two numbers, without actually performing the factorization of either.
This page has pretty much all you need to know about it:
http://en.wikipedia.org/wiki/Euclidean_algorithm

"Monte Carlo" method(s)

What is a "Monte Carlo" (factorization) method?
Essentially one where luck can play a major part - hence the rather exciting-sounding connection with gambling! This is the case whenever the process is seeded by some (pseudo-or-otherwise) random number generation (eg for the benefit of parallelization), or when there is inherent unpredictability in the algorithm itself.
Note also, that for the Rho method, for example, this Monte Carlo nature of the algorithm results in performance gains caused by a statistical effect akin to the so-called "Birthday Paradox". This latter, in turn, is the term given to the surprisingly high observed probability of at least one coincidence of two different random variables in the same range (eg two or more people's birthdays)
The following Wikipedia link talks about Monte Carlo effect in general:
http://en.wikipedia.org/wiki/Monte_Carlo_method
and here is a link to the Birthday Paradox:
http://en.wikipedia.org/wiki/Birthday_Paradox

Rho method

Today I am going to be talking a little bit about Pollard's Rho algorithm, which, for such an apparently simple algorithm, is surprisingly good at finding moderate-size factors of general numbers.

Some reference links:
http://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
[as I mentioned briefly once before]
http://mersennewiki.org/index.php/Rho_Factorization_Method

Also, the pseudocode immediately below is taken from a PD document by Connelly Barnes of Oregon State University
http://oregonstate.edu/~barnesc/documents/factoring.pdf

function pollardRho(N)
# Initial values x(i) and x(2*i) for i = 0.
xi := 2
x2i := 2
do
# Find x(i+1) and x(2*(i+1))
xiPrime := xi ^ 2 + 1
x2iPrime := (x2i ^ 2 + 1) ^ 2 + 1
# Increment i: change our running values for x(i), x(2*i).
xi := xiPrime % N
x2i := x2iPrime % N
s := gcd(xi - x2i, N)
if s <> 1 and s <> N then
return s, N/s
end if
end do
end function


The Rho method is also the main method used by GMP's demo factorization sample program (or was last time I looked) - immediately below:

void
factor_using_pollard_rho (mpz_t n, int a_int, unsigned long p)
{
mpz_t x, x1, y, P;
mpz_t a;
mpz_t g;
mpz_t t1, t2;
int k, l, c, i;

if (flag_verbose)
{
printf ("[pollard-rho (%d)] ", a_int);
fflush (stdout);
}

mpz_init (g);
mpz_init (t1);
mpz_init (t2);

mpz_init_set_si (a, a_int);
mpz_init_set_si (y, 2);
mpz_init_set_si (x, 2);
mpz_init_set_si (x1, 2);
k = 1;
l = 1;
mpz_init_set_ui (P, 1);
c = 0;

while (mpz_cmp_ui (n, 1) != 0)
{
S2:
if (p != 0)
{
mpz_powm_ui (x, x, p, n); mpz_add (x, x, a);
}
else
{
mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
}
mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
c++;
if (c == 20)
{
c = 0;
mpz_gcd (g, P, n);
if (mpz_cmp_ui (g, 1) != 0)
goto S4;
mpz_set (y, x);
}
S3:
k--;
if (k > 0)
goto S2;

mpz_gcd (g, P, n);
if (mpz_cmp_ui (g, 1) != 0)
goto S4;

mpz_set (x1, x);
k = l;
l = 2 * l;
for (i = 0; i < k; i++)
{
if (p != 0)
{
mpz_powm_ui (x, x, p, n); mpz_add (x, x, a);
}
else
{
mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
}
}
mpz_set (y, x);
c = 0;
goto S2;
S4:
do
{
if (p != 0)
{
mpz_powm_ui (y, y, p, n); mpz_add (y, y, a);
}
else
{
mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
}
mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
}
while (mpz_cmp_ui (g, 1) == 0);

if (!mpz_probab_prime_p (g, 3))
{
do
{
mp_limb_t a_limb;
mpn_random (&a_limb, (mp_size_t) 1);
a_int = (int) a_limb;
}
while (a_int == -2 || a_int == 0);

if (flag_verbose)
{
printf ("[composite factor--restarting pollard-rho] ");
fflush (stdout);
}
factor_using_pollard_rho (g, a_int, p);
break;
}
else
{
mpz_out_str (stdout, 10, g);
fflush (stdout);
fputc (' ', stdout);
}
mpz_div (n, n, g);
mpz_mod (x, x, n);
mpz_mod (x1, x1, n);
mpz_mod (y, y, n);
if (mpz_probab_prime_p (n, 3))
{
mpz_out_str (stdout, 10, n);
fflush (stdout);
fputc (' ', stdout);
break;
}
}

mpz_clear (g);
mpz_clear (P);
mpz_clear (t2);
mpz_clear (t1);
mpz_clear (a);
mpz_clear (x1);
mpz_clear (x);
mpz_clear (y);
}


This Java code is taken directly from superfac9:

BigInteger factorizerho(BigInteger n) {
BigInteger loop = new BigInteger("1");
BigInteger x = new BigInteger("5");
BigInteger y = new BigInteger("2");

while (n.gcd(x.subtract(y)).compareTo(ONE) == 0 && loop.compareTo(TENTHOUSAND) < 0) {
x = x.multiply(x).add(ONE).mod(n);
x = x.multiply(x).add(ONE).mod(n);
y = y.multiply(y).add(ONE).mod(n);
loop = loop.add(ONE);
}

return n.gcd(x.subtract(y));

}


Note that there are several variants (given different start-point values for x and y for example), but the basic "Monte-Carlo" principle remains the same.

Notable successes with the Rho method or its variants include Brent and Pollard's factorization of the eighth Fermat number in 1980. The 16-digit factor of F8 took about 2 hours to find on a UNIVAC 1100/42, apparently.

F12-WEP

And here is the WEP algorithm running (unsuccessfully) for an hour or so looking for new factors of F12:

tiggatoo:~/math/wec james$ time ./factorize3.gmp -P4096 -T1000 114689 26017793 63766529
P4096 T1000

real 59m16.287s
user 12m32.262s
sys 0m4.508s

M521^2 by WE

Here is a factorization illustrating the automatic algebraic factoring capability of superfac9, due to WE algorithm:

tiggatoo:~/math james$ java superfac9
f(2^521-1)^2

[47125446914534694131579097993419809976955095716785201420286055195012674566357244479460731079205201122720511132925006540350105785156086431086764996857554304847155991333706718342307167456986269662311038377104760933477381254100896222805785374204495333936040246318307567782851014765052850751581472024524956029996236801]
wanless...
aprtcle: 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
aprtcle: 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
duration: 115 seconds

Near-repdigit factorization(s)

And speaking of repunits (or repdigits - the factorization of these latter reduce to the former), there is also Makoto Kamada's page of factorizations of near-repdigits.
http://homepage2.nifty.com/m_kamada/math/factorizations.htm
Clearly this enlarges the options for factoring candidates significantly.

Repunit factorization

Here is a link to Yousuke Koide's page on repunit factorization(s). Repunits are numbers specifically whose decimal representation consists of all 1's.
http://www.h4.dion.ne.jp/~rep/
It would appear that the lowest repunit whose prime factorization is yet unknown, is R239 (see also Torbjorn Granlund's page at Swox).
http://swox.com/~tege/fac10m.txt
[Note that Swox is the organization behind the widely-used GMP math library
http://gmplib.org/]

GGNFS bounce!

GGNFS is cool! (and deservedly a standard reference implementation). But it would appear I have been a victim of the GGNFS "bounce", as Bob Backstrom would call it! I had almost collected enough relations to proceed to the linear algebra stage of my factorization attempt of C127_113_36 (for the xyyxf project), when this happened. Apparently, discarding heavy relations speeds up the later stage, but it's going to take me a while to replace those discarded relations. Bob assures me that this bounce will not happen again, and also suggests (as in fact I wondered) that it may well be possible to proceed to the LA immediately with msieve. However, since I'm fairly new at this (not having run even straight GGNFS before on a number this size), for now I think I've decided to wait, watch, and hopefully learn...

Benchmarks: Java vs C (Addendum)

Of course Assembler can be even faster than C. George Woltman's big number library (gwnum) - part of Prime95/mprime - for example, seems very fast, perhaps at least as fast again as GMP, particularly for extrabig numbers.
Gwnum also uses advanced FFT techniques, and because it is partly written in assembler, not surprisingly, it is not portable - only being available for the most common platform, x86.
It would be interesting to rewrite my WEP application using gwnum...
http://www.mersenne.org/freesoft.htm
http://en.wikipedia.org/wiki/Fft

Benchmarks: Java vs C (Part 4)

Here are the promised benchmarks, running random-based WEP on input numbers of, respectively, 100, 1000 and 10000 bits.
If anyone would care to (and has the patience to) extend the tables, feel free to send me results!
Result: Java is, very approximately, an order of magnitude, only, slower than C.

tiggatoo:~/math/wec james$ time java we2tpr2 < P109.in.2.txt
[109]
[1]
base#0
elapsed=0s factor=3 A=645264517869175453037177387327
base#10746
elapsed=43s factor=104124649 A=245093880973395622884386051056621
2077756847362348863128179
duration: 43 seconds


real 0m46.195s
user 0m5.185s
sys 0m0.510s

tiggatoo:~/math/wec james$ time ./factorize3.gmp -P109 -T10748 3
P109 T10748

real 0m7.176s
user 0m0.793s
sys 0m0.011s

tiggatoo:~/math/wec james$ time java we2tpr2 < P1279.in.2.txt
[1279]
[1]
base#0
elapsed=2s factor=3 A=530978962560898025744780713800683
base#11
elapsed=13s factor=706009 A=1243358263022091534365334367195221
^Cse#100
real 2m34.366s
user 0m17.542s
sys 0m0.222s

tiggatoo:~/math/wec james$ time ./factorize3.gmp -P1279 -T114 3 706009
P1279 T114

real 0m29.889s
user 0m3.344s
sys 0m0.022s

tiggatoo:~/math/wec james$ time java we2tpr2 < P11213.in.2.txt
[11213]
[1]
base#0
elapsed=1510s factor=3 A=93613194278166973511368466002199
^Cse#2
real 45m19.926s
user 5m14.923s
sys 0m2.640s

tiggatoo:~/math/wec james$ time ./factorize3.gmp -P11213 -T3 3
P11213 T3

real 5m44.733s
user 0m40.864s
sys 0m0.247s

Benchmarks: Java vs C (Part 3)

Here is the (nearly) same program (factorize3.mac.c) coded in C (also copyright JGW). It needs to be compiled with -lgmp, the GMP math library.

/* Factoring with WEP method using random base(s).
*/

#include
#include
#include

#include

#include "gmp.h"

int flag_verbose = 0;

int
factor_using_random_wep (mpz_t s, unsigned long p, unsigned long T)
{
mpz_t V, v1, v2, v3, b, bb, bbb, bbbb, wanless;
long numtrials;
int flag;

gmp_randstate_t rstate;



if (flag_verbose)
{
printf ("[wep ");
printf ("s=");
mpz_out_str (stdout, 10, s);
printf ("\tp=");
printf ("%ld", p);
printf ("\tT=");
printf ("%ld", T);
printf ("]\n\r");
fflush (stdout);
}


if (mpz_probab_prime_p (s, 3))
{
printf ("P%ld\tT0\t", p);
mpz_out_str (stdout, 10, s);
fflush (stdout);
flag=2;
return(flag);
}



gmp_randinit (rstate, GMP_RAND_ALG_LC, 128);

{
#if HAVE_GETTIMEOFDAY
struct timeval tv;
gettimeofday (&tv, NULL);
gmp_randseed_ui (rstate, tv.tv_sec + tv.tv_usec);
#else
time_t t;
time (&t);
gmp_randseed_ui (rstate, t);
#endif
}






numtrials=0;

mpz_init_set_si (b, 1);
mpz_init_set_si (bb, 1);
mpz_init_set_si (bbb, 1);
mpz_init_set_si (bbbb, 1);
mpz_init_set_si (V, 1);
mpz_init_set_si (v1, 1);
mpz_init_set_si (v2, 1);
mpz_init_set_si (v3, 1);
mpz_init_set_si (wanless, 2);

while (mpz_cmp (wanless, s) < 0)
mpz_mul_ui (wanless, wanless, 2);


while (numtrials < T && (mpz_cmp_ui (V, 1) == 0 || mpz_cmp (V, s) == 0))
{


mpz_urandomb (b, rstate, 100L);
mpz_mul_ui(bb, b, p);
mpz_mul_ui(bbb, bb, 2);
mpz_add_ui(bbbb, bbb, 1);


mpz_powm (v1, bbbb, wanless, s);
mpz_powm (v2, v1, s, s);
mpz_sub (v3, v2, v1);
mpz_gcd (V, s, v3);

if (flag_verbose)
if (numtrials%1000 == 0)
{
printf ("numtrials=%ld\tb=", numtrials);
mpz_out_str (stdout, 10, b);
printf ("\r");
fflush (stdout);
}

numtrials++;

}

if (flag_verbose)
printf("\n\r");


if (mpz_cmp_ui (V, 1) > 0 && mpz_cmp (V, s) < 0)
{
printf ("P%ld\tT%ld\t", p, numtrials);
mpz_out_str (stdout, 10, V);
printf("\tbase=");
mpz_out_str (stdout, 10, bbbb);
printf ("\n");
fflush (stdout);

flag=3;
}

else
{
printf ("P%ld\tT%ld\n", p, T);
flag=0;
}



mpz_clear (b);
mpz_clear (bb);
mpz_clear (bbb);
mpz_clear (bbbb);
mpz_clear (V);
mpz_clear (v1);
mpz_clear (v2);
mpz_clear (v3);
mpz_clear (wanless);

return (flag);

}

main (int argc, char *argv[])
{
mpz_t r, s, t, F, f;
unsigned long p, T;
int i;
int flag;

if (argc > 1 && !strcmp (argv[1], "-v"))
{
flag_verbose = 1;
argv++;
argc--;
}

mpz_init (r);
mpz_init (s);
mpz_init (t);
mpz_init (F);
mpz_init (f);

mpz_set_ui (F, 1);
mpz_set_ui (f, 1);

if (argc > 1)
{
p = 0;
for (i = 1; i < argc; i++)
{
if (!strncmp (argv[i], "-M", 2))
{
p = atoi (argv[i] + 2);
mpz_set_ui (t, 1);
mpz_mul_2exp (t, t, p);
mpz_sub_ui (t, t, 1);
}
else if (!strncmp (argv[i], "-P", 2))
{
p = atoi (argv[i] + 2);
mpz_set_ui (t, 1);
mpz_mul_2exp (t, t, p);
mpz_add_ui (t, t, 1);
}
else if (!strncmp (argv[i], "-T", 2))
{
T = atoi (argv[i] + 2);
}
else
{
mpz_set_str (f, argv[i], 0);
mpz_mul (F, F, f);
}
}


mpz_mod (r, t, F);

if (mpz_cmp_si (r, 0) != 0) {
printf ("Wrong known factors!\n");
flag =1;
}
else {

mpz_div (s, t, F);

flag = factor_using_random_wep (s, p, T);

}
}

mpz_clear(r);
mpz_clear(s);
mpz_clear(t);
mpz_clear(F);
mpz_clear(f);

exit (flag);
}

Benchmarks: Java vs C (Part 2)

Here is the program, we2tpr2.java, coded in Java (copyright JGW, all rights reserved :) ). Note how elegant, and convenient, the BigInteger facility is - this has been a standard feature of Java from its inception. Java is almost unique in that respect - most other languages rely on add-in libraries to provide big number capability.

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

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

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

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

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




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

we2tpr2 we2tpr2inst = new we2tpr2();

char c;
String sInput;

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;
p = we2tpr2inst.eval(s);

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



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;
known = we2tpr2inst.eval(s);

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



n = we2tpr2inst.mersenneplustwo(p);
if (n.remainder(known).compareTo(zero) > 0) {
System.out.println("Wrong known factors!");
return;
}
else
n = n.divide(known);

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

we2tpr2inst.factorize(n, p);

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

duration = (finishtime-starttime)/1000;

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

System.out.println();
}


public BigInteger mersenneplustwo(BigInteger p) {
BigInteger i = new BigInteger("0");

n = two;
for (i = one; i.compareTo(p) < 0; i = i.add(one))
n = n.multiply(two);

n = n.add(one);

return n;
}


public boolean factorize(BigInteger n, BigInteger p) {
boolean prime = false;
BigInteger numtested = new BigInteger("0");
BigInteger T = new BigInteger("1");
BigInteger b = new BigInteger("1");
BigInteger A = new BigInteger("2");
BigInteger wanless = new BigInteger("2");

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

// workaround - apparent java bug in modPow - JGW
if (n.compareTo(two) < 0)
return false;

if (n.remainder(two).compareTo(zero) == 0) {
System.out.println(two.toString());
return true;
}
// end workaround

while (wanless.compareTo(n) < 0)
wanless = wanless.multiply(two);

Random r = new Random();

numtested = zero;

while (T.compareTo(one) == 0 || T.compareTo(n) == 0) {
// changed JW 2005-3-23
A = new BigInteger(hundred.intValue(), r);
A = (A.multiply(two).multiply(p)).add(one);

// added JGW 2006-06-09
System.out.print("base#" + numtested + '\r');



// changed DT 2005-2-20
b = A.modPow(wanless, n);
T = n.gcd(b.modPow(n, n).subtract(b));

numtested = numtested.add(one);
}





if (T.compareTo(one) > 0 && T.compareTo(n) < 0) {

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

duration = (finishtime-starttime)/1000;


System.out.println();
System.out.println("elapsed=" + duration + "s" + '\t' + "factor=" + T.toString() + '\t' + "A=" + A.toString() + '\t');
factorize(n.divide(T), p);
}

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) {
BigInteger n = new BigInteger("0");

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

return n;
}


}

Benchmarks: Java vs C (Part 1)

Over the next few days I am going to be benchmarking a small factorization program (random-based WEP algorithm) written (approximately) identically in Java (using Java's built-in BigInteger feature) and in C (using the GMP math library), and seeing how the two compare for various size inputs.
But first some background information about Java:
http://en.wikipedia.org/wiki/Java_%28programming_language%29
and a picture (from Wikipedia) of the Java mascot, Duke:

and (also from Wikipedia) of the creator of Java (at Sun Microsystems), James Gosling:

Fermat number factorizers

Here is a nice page detailing historical (and recent) figures in the search for factors of Fermat numbers. Note that the first computer calculations are dated around 1953.
http://www.fermatsearch.org/history.php

What is a "trapdoor" function?

Term first coined by Diffie and Hellman in 1976, refers to a function that is easy to compute forwards, yet very difficult to invert, ie compute backwards. Large integer multiplication/factorization serves this purpose for RSA.
see link:
http://en.wikipedia.org/wiki/Trapdoor_function

Cryptography

I am currently listening (for free, on iTunes - search on "math 55") to Professor James Demmel's Discrete Mathematics series of lectures at University of California, Berkeley - highly recommended (so far). Anyway, Lecture 18 audio is relevant to today's post, which is on public-key cryptography, and, specifically RSA. The reason cryptography is involved is because RSA relies on the difficulty of factoring large semiprimes for keeping encoded messages secret. Some links:
http://en.wikipedia.org/wiki/Public-key_cryptography
http://en.wikipedia.org/wiki/Rsa
http://www.cs.berkeley.edu/~demmel/ma55_Fall07/LectureNotes/Lecture_15_Oct_03.txt

History of factorization records by ECM

Original at:
http://www.loria.fr/~zimmerma/records/factor.html
(http://www.loria.fr/~zimmerma/records/ecmrecord.jpg)
Reproduced by kind permission of Paul Zimmermann

222-digit SNFS completed with msieve

From the following link (by Greg Childers):
http://mersenneforum.org/showthread.php?p=115322#post115322
"The NFSNet factorization of 5^317-1 (SNFS difficulty of 222 digits) has been completed with msieve."
"...the runtime was just over 6 days"
Apparently this represents a real advance in speed of processing over the established suite for postprocessing large NFS sieve jobs. Thanks Jason (Papadopoulos)!

Msieve 1.28

New version of msieve (by Jason Papadopoulos) now available.
Download from:
http://www.boo.net/~jasonp/qs.html

Mersennewiki

...also, there is the Mersennewiki, which has a factorization section at:
http://mersennewiki.org/index.php/Factorization

mersenneforum

I'd just like to draw people's attention to the 'factoring' section of the mersenneforum (at the following link) - it has certainly made very interesting reading for me at times in the past...
http://mersenneforum.org/forumdisplay.php?f=19

Msieve v1.27

New version of msieve (by Jason Papadopoulos) now available.
Download from:
http://www.boo.net/~jasonp/qs.html
Announcement at:
http://mersenneforum.org/showthread.php?p=115491#post115491

Fermat number factorizations

Speaking of Fermat number factorization - this page:
http://www.prothsearch.net/fermat.html
lists the current state-of-play.
Note that only F5-F11 have been completely factored. The smallest Fermat number with no known factors is F14, and the smallest Fermat number whose compositeness (ie factorizability) has not been proven is F33.
[for additional/preparatory reference - Wikipedia article on Fermat numbers: http://en.wikipedia.org/wiki/Fermat_number]

Factorization of tenth Fermat number in 1995

From Richard Brent's page:
http://wwwmaths.anu.edu.au/~brent/pub/pub161.html
"We describe the complete factorization of the tenth Fermat number F10 by the elliptic curve method (ECM). The tenth Fermat number is a product of four prime factors with 8, 10, 40 and 252 decimal digits. The 40-digit factor was found after about 140 Mflop-years of computation"