// This accepts no inputs. It times calculations and prints results.
// From the table of exponents p it calculates 2^p+1
// and uses trial division to look for divisors 2kp+1
// starting at 10^n and estimates the time for trial
// division between 10^n and 10^(n+5), using the usual
// trial division and a faster powerMod method.
// Then it uses the Pollard Rho method to estimate the
// time to make a similar factorization attempt.
// Then it uses the ECM method to estimate the time to
// make a similar factorization attempt.
// Then it calculates the probability of finding a factor
// in the process of completing that pass, using Gillies'
// log(b/a) estimate from "Three New Mersenne Primes and
// a Statistical Theorey", Math. Comp. v.18 pp 93-95.
import java.math.BigInteger;
//import java.util.Date; // to try exec and not duplicate java.util.*
import java.text.*;
import java.util.*;
import java.io.*;
// Copied from http://www.javaworld.com/javaworld/jw-12-2000/jw-1229-traps.html
// Listing 4.5 GoodWindowsExec.java
// Then modified to fit the application.
// I'm STILL worried about the real precise correctness of this.
// This process is given an input stream, a flag indicating whether or not to
// just throw away the contents of the stream and the size of the buffer to
// create. Then it will block waiting for eol terminated lines of input,
// read each when ready, and if not to discard the line then print the line
// on stdout.
// Two copies of this will be started, one for stdout, one for stderr from GMP.
class readStream extends Thread {
InputStream is;
boolean discard;
int size;
readStream(InputStream is, boolean discard, int size) {
this.is = is;
this.discard = discard;
this.size = size;
}
public void run() {
// if you are curious
// System.out.println("proc starting");
try {
BufferedReader br = new BufferedReader(new InputStreamReader(is), size);
String line=null;
while((line = br.readLine()) != null)
if(!discard)
System.out.println(line);
// don't leave open streams lost out there
// I believe that closing br will close is.
br.close();
} catch (IOException ioe) {
ioe.printStackTrace();
}
// if you are curious
// System.out.println("proc finishing");
}
}
public class stats23 {
static final BigInteger ONE = new BigInteger("1");
static final BigInteger TWO = new BigInteger("2");
static final BigInteger TEN = new BigInteger("10");
static final BigInteger THIRTY = new BigInteger("30");
static final int[] p={
2203, 2281, 3217, 4253, 4423,
9689, 9941, 11213, 19937, 21701,
23209, 44497, 86243, 110503, 132049,
216091, 756839, 859433, 1257787, 1398269,
2976221, 3021377, 6972593, 13466917, 20996011,
24036583, 25964951, 30402457};
// each of these n correspond to a p above.
// others have already searched up to 10^n
static final int[] n={
40, 40, 40, 35, 35,
35, 35, 35, 35, 30,
30, 30, 30, 30, 30,
25, 25, 20, 20, 20,
20, 20, 15, 15, 15,
15, 15, 15};
// how many repetitions of ECM are needed for nn digit pass
// taken directly from the table http://www.alpertron.com.ar/ECM.HTM
// as was the code for ECM I believe, but JGW ported that.
static final int[] ecmiterations={
25, // 15 digits
90, // 20 digits
300, // 25 digits
700, // 30 digits
1800, // 35 digits
5100, // 40 digits
10600, // 45 digits
19300, // 50 digits
49000, // 55 digits
124000, // 60 digits
210000, // 65 digits
340000}; // 70 digits
// Round double with an Enn exponent to 3 digits plus Enn.
// NumberFormat doesn't have options to display numbers
// with a trailing exponent.
// So, remove the exponent using repeated division by 10,
// let NumberFormat format the number without the exponent,
// and then append the Enn to the end of the resulting string.
// WARNING! This will not correctly handle negative exponents
// negative numbers or zero.
// But for this program at the moment these aren't used.
// This should be simple to fix, but I have claimed that before.
static String roundexp (double value) {
int power_of_ten = 0;
while (value >= 9.995) { // avoid rounding to 10.00Enn
value /= 10.0;
power_of_ten++;
}
// round to two digits after the decimal point
NumberFormat fmt3 = new DecimalFormat("0.00");
return fmt3.format(value)+("E"+power_of_ten+" ").substring(0,3);
}
public static void main (String args[])
throws java.io.IOException {
int i, j;
Date d; long starttime, finishtime, duration;
final double MILLISECONDS_PER_MONTH = 1000.0*3600.0*24.0*30.0;
BigInteger mpp, divisor, onep, twop, fourp, sixp, result;
int divreps, modreps, rhoreps, ecmreps, gmpreps;
double divtime, modtime, rhotime, ecmtime, gmptime;
double rhoduration=1.0; double rhocoeff=1.0; double rhocoeffproduct=1.0;
double ecmduration=1.0; double ecmcoeff=1.0; double ecmcoeffproduct=1.0;
double gmpduration=1.0; double gmpcoeff=1.0; double gmpcoeffproduct=1.0;
ecm31 ecm31inst = new ecm31();
// Round xxxcoeff display to 4 digits
NumberFormat fmt4 = new DecimalFormat("0.000");
double probability;
// Round probability display to 3 digits
NumberFormat fmt3 = new DecimalFormat("0.00");
// switch on or off the display of xxxcoeff in output
boolean displaycoeffs = true;
// as the size of primes grows the time needed to do trial divisions
// grows, so each time a pass took more than 20 seconds the number of
// trial divisions=reps will be reduced by a factor of 2 in future
// passes. Each of these was chosen by trial and error checking.
divreps = 524288; // 2^19 so neatly divisible by two many times
modreps = 524288; // 2^19 so enough repetitions to get good timing
rhoreps = 2048; // enough repetitions to get good timing
// if you are going to change ecmreps you must understand the ecm code
// because increasing ecmreps will invoke larger digit searching, make
// slower passes, etc. I believe it would be possible to search curves
// starting at 26 and finishing before or at curve 90, or similar changes,
// but these will slow down the timing considerably and it is already slow.
ecmreps = 25; // the "first" ecm pass tests curves from 1 to 25.
gmpreps = 25; // likewise
// label the table
if(displaycoeffs)
System.out.println(" p n div mod rho c ecm c gmp c prob");
else
System.out.println(" p n div mod rho ecm gmp probability");
for(i=0;i 20000) divreps = divreps/2;
if (divreps < 15) divreps = 15;
// The trial divison method via mod for 2kp+1, skipping multiples of 3 and 5
// For large exponents mod is significantly faster than using division
// This terrific optimization courtesy JGW who was given it by someone else
d = new Date(); starttime = d.getTime(); // start the clock
divisor=TEN.pow(n[i]).divide(onep).multiply(onep).add(ONE);
while(divisor.mod(THIRTY).intValue() != 0)
divisor = divisor.subtract(onep);
divisor = divisor.add(onep);
for(j=0;j 20000) modreps = modreps/2;
if (modreps < 15) modreps = 15;
// The Pollard Rho method for 2kp+1
// This is optimised for that special case based on information in
// "Prime Numbers and Computer Methods for Factorization",
// Hans Riesel, 2nd Ed. Birkhauser. pp. 185-186.
// In this special case by using x^p instead of x^2 the algorithm is
// expected to be sqrt(p) faster.
// By using 2.4*sqrt(10^(n+5))/sqrt(p) iterations it is conjectured to
// have >90% chance of finding a factor if one exists in the range.
// I hope there aren't any errors in this, but I've said that before.
// Unfortunately, even dividing rhoreps by 2 each time the calculation
// takes too long is not enough to save us here. Finally a single
// iteration of the rho calculation will take to long and we need an
// alternate method of estimating the time to complete the full pass.
// To do this we predict calculation times for a single step when
// rhoreps has been reduced to zero, computing the coefficient rhocoeff.
// This coefficient was obtained by calculating the geometric mean of
// rhotime/rhoreps = c p^2 log(p)
// for measured values of rhotime and rhoreps, to extrapolate rhotime
// when actually measuring it would take far too long.
// Rhocoeff remains fairly stable for medium sized values of p
// but for the largest measured values of p the value begins to grow,
// perhaps because of ever increasing memory demands of the calculations.
// This should be explored to discover whether there is an error lurking
// in this or whether it really is memory demands.
// Limited testing with java options of -Xms256m -Xmx512m did not seem
// to change the timing for the smaller values of p but larger values
// were not explored and certainly should be.
// This REALLY needs to be calculated from the measured running times
// This should probably compute the geometric mean of calculated values
// rhocoeff, and then use that value when times must be estimated.
if(rhoreps==1) { // estimate time for this pass
rhoduration = (double)rhocoeff*p[i]*p[i]*Math.log(p[i]);
rhotime = (double)rhoduration/MILLISECONDS_PER_MONTH*
2.4*Math.sqrt(Math.pow(10.0,n[i]+5.0))/Math.sqrt(p[i]);
}
else { // measure time for this pass
d = new Date(); starttime = d.getTime(); // start the clock
BigInteger rhox = new BigInteger("5");
BigInteger rhoy = new BigInteger("2");
for(j=0;j 20000) rhoreps = rhoreps/2;
if (rhoreps < 1) rhoreps = 1;
if(rhoreps==1)
rhocoeff = Math.pow(rhocoeffproduct,1.0/(i+1)); // i+1 values->(i+1)th root
}
// The ECM method
// IF my understanding of how this ecm code works, questionable,
// One difficulty with this is that what we really want to know is
// how long it will take for long unproductive searches for factors
// with lots of digits. But that will take far too long to measure.
// Instead we measure the first few curves looking for small factors.
// Unfortunately this will usually find factors when it does that.
// But perhaps this will give a plausible estimate of the run times.
// Unfortunately even doing the first pass still takes far too long
// for larger values of p. So we cut corners again and don't even
// do all the curves needed for the first pass as p gets larger and
// estimate how long the pass would actually take. Even that is not
// enough as p grows and we then turn to estimating the estimate,
// based on timing of earlier passes and partial passes we calculate
// a coefficient and use that to estimate the pass for larger p.
// 32k limit on array sizes in ecm code, estimate times beyond that
if(p[i]>=32000) { // estimate time for this pass
ecmcoeff = Math.pow(ecmcoeffproduct,1.0/11); // 11 values->11th root
ecmduration = (double)(ecmcoeff*p[i]*p[i]);
}
else { // measure time for this pass
BigInteger startCurve = ONE;
BigInteger finishCurve = new BigInteger(""+ecmreps);
d = new Date(); starttime = d.getTime(); // start the clock
ecm31inst.factorize(mpp, startCurve, finishCurve);
d = new Date(); finishtime = d.getTime(); // stop the clock
duration = finishtime-starttime;
// extrapolate time this full pass would take
// if we did fewer curves to save time
ecmduration = (double)duration*ecmiterations[0]/ecmreps;
// accumulate data for geometric mean calculation
ecmcoeff = (double)ecmduration/p[i]/p[i];
ecmcoeffproduct *= ecmcoeff;
// if that pass took more than 1/2 hour do fewer next time
if(duration > 1800000) ecmreps = ecmreps/2;
if (ecmreps < 1) ecmreps = 1;
}
// estimate the time for the pass we actually need to do
// ((n[i]+5)-15)/5 gives the "number" of the that pass
// 15->pass 0, 20->pass 1, 25->pass 2, etc.
ecmtime = (double)ecmduration/ // time for this pass
MILLISECONDS_PER_MONTH* // turn that into into months
Math.pow(10,((n[i]+5)-15)/5); // 10x for each higher pass
// now try running the external GMP executable
// Version 0.2, run GMP executable, gmpreps curves, quiet mode, B1 B2 bounds
// using separate processes to handle stdout and stderr, waitFor() to finish.
if(gmpreps==1) { // estimate time for this pass
gmpduration = (double)(gmpcoeff*p[i]*p[i]);
}
else {
try {
Process proc = Runtime.getRuntime().exec("./ecm -c "+gmpreps+" -q 2000 200000");
// Read and don't discard stderr, println that if it ever happens.
// Unless something goes wrong GMP should never write to stderr.
// 8k should be MORE than enough for a stderr line, that is default
readStream stderrRead = new readStream(proc.getErrorStream(), false,8192);
// read and discard stdout, change to false if you want to see it.
// p[i]/3 SHOULD be big enough to hold the 2^p[i]+1 factorization.
readStream stdoutRead = new readStream(proc.getInputStream(), true,p[i]/3);
// start the processes running
stderrRead.start();
stdoutRead.start();
// get stream to prepare to send mpp to stdin of GMP
OutputStream stdin = proc.getOutputStream();
OutputStreamWriter osw = new OutputStreamWriter(stdin);
// apparently we don't need to set a buffer size for writing
BufferedWriter bw = new BufferedWriter(osw);
System.out.print("GMP "+'\r');
d = new Date(); starttime = d.getTime(); // start the clock
// deliver mpp to the standard input of GMP ECM
bw.write(""+mpp); bw.newLine(); bw.flush(); bw.close();
// if you are curious
// System.out.println("the waitFor begins");
int exitVal = proc.waitFor(); // wait for GMP to finish
// if you are curious
// System.out.println("the waitFor ended!");
d = new Date(); finishtime = d.getTime(); // stop the clock
duration = finishtime-starttime;
// if you are curious
// System.out.println("Duration: " + duration);
// System.out.println("Process exitValue: " + exitVal);
// extrapolate time this full pass would take
// if we did fewer curves to save time
gmpduration = (double)duration*ecmiterations[0]/gmpreps;
// accumulate data for geometric mean calculation
gmpcoeff = (double)gmpduration/p[i]/p[i];
gmpcoeffproduct *= gmpcoeff;
// if that pass took more than 1/2 hour do fewer next time
if(duration > 1800000) gmpreps = gmpreps/2;
if (gmpreps < 1) gmpreps = 1;
if(gmpreps==1)
gmpcoeff = Math.pow(gmpcoeffproduct,1.0/(i+1)); // i+1 values->(i+1)th root
}
catch (Throwable t) {
t.printStackTrace();
}
}
gmptime = (double)gmpduration/ // measured time for this run
MILLISECONDS_PER_MONTH* // turn that into into months
Math.pow(10,((n[i]+5)-15)/5); // 10x for each higher pass
// Probability of finding a factor on the actual pass, per Gillies
probability = Math.log((double)(n[i]+5)/(double)n[i]);
// Format and display the timing data and chance of finding a factor
// Elapsed time in months to do the pass using each method
// and the probability of finding a factor in that pass
if(displaycoeffs)
System.out.println((p[i]+" ").substring(0,8)+" "+n[i]+" "+
roundexp(divtime)+" "+
roundexp(modtime)+" "+
roundexp(rhotime)+" "+fmt4.format(rhocoeff*1.0E6)+" "+
roundexp(ecmtime)+" "+fmt4.format(ecmcoeff)+" "+
roundexp(gmptime)+" "+fmt4.format(gmpcoeff*1.0E2)+" "+
fmt3.format(probability));
else
System.out.println((p[i]+" ").substring(0,8)+" "+n[i]+" "+
roundexp(divtime)+" "+
roundexp(modtime)+" "+
roundexp(rhotime)+" "+
roundexp(ecmtime)+" "+
roundexp(gmptime)+" "+
fmt3.format(probability));
}
}
}