// 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<p.length;i++) { // the 28 known Merssenne exponents

         onep = BigInteger.valueOf(p[i]);
         twop = onep.add(onep);
         fourp = twop.add(twop);
         sixp = twop.add(fourp);
         mpp = TWO.pow(p[i]).add(ONE);

         // The trial division method for 2kp+1, skipping multiples of 3 and 5

         d = new Date(); starttime = d.getTime(); // start the clock
         divisor=TEN.pow(n[i]).divide(onep).multiply(onep).add(ONE);
         // Every 3rd trial will be a multiple of 3.
         // Every 5th trial will be a multiple of 5.
         // and neither of those can ever be a factor of mpp.
         // Speed can be approximately doubled by skipping those
         // but we need to get the the process started at the right point
         // so decrement the divisor a maximum of 29 times until we find
         // a value of divisor==0 mod 2 and 3 and 5 to start correctly.
         // Then we increment divisor once and we are ready to begin.
         while(divisor.mod(THIRTY).intValue() != 0)
            divisor = divisor.subtract(onep);
         divisor = divisor.add(onep);

         for(j=0;j<divreps;j+=15) {
            result = mpp.mod(divisor); // mod=1
            divisor = divisor.add(sixp); // mod=2,3,4,5,6 skipped
            result = mpp.mod(divisor); // mod=7
            divisor = divisor.add(fourp); // mod=8,9,10 skipped
            result = mpp.mod(divisor); // mod=11
            divisor = divisor.add(twop); //mod=12 skipped
            result = mpp.mod(divisor); // mod=13
            divisor = divisor.add(fourp); // mod=14,15,16 skipped
            result = mpp.mod(divisor); // mod=17
            divisor = divisor.add(twop); // mod=18 skipped
            result = mpp.mod(divisor); // mod=19
            divisor = divisor.add(fourp); // mod=20,21,22 skipped
            result = mpp.mod(divisor); // mod=23
            divisor = divisor.add(sixp); // mod=24,25,26,27,28 skipped
            result = mpp.mod(divisor); // mod=29
            divisor = divisor.add(twop); // mod 30 skipped
         }
         d = new Date(); finishtime = d.getTime();  // stop the clock
         duration = finishtime-starttime;
         divtime = (double)duration/divreps/MILLISECONDS_PER_MONTH*
            Math.pow(10.0,n[i]+5.0)/(2.0*p[i]);

         // if that pass took more than 20000 milliseconds do fewer next time
        if(duration > 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<modreps;j+=15) {
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=1
            divisor = divisor.add(sixp); // mod=2,3,4,5,6 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=7
            divisor = divisor.add(fourp); // mod=8,9,10 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=11
            divisor = divisor.add(twop); //mod=12 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=13
            divisor = divisor.add(fourp); // mod=14,15,16 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=17
            divisor = divisor.add(twop); // mod=18 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=19
            divisor = divisor.add(fourp); // mod=20,21,22 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=23
            divisor = divisor.add(sixp); // mod=24,25,26,27,28 skipped
            result = TWO.modPow(onep, divisor).add(ONE).mod(divisor); // mod=29
            divisor = divisor.add(twop); // mod 30 skipped
         }
         d = new Date(); finishtime = d.getTime();  // stop the clock
         duration = finishtime-starttime;
         modtime = (double)duration/modreps/MILLISECONDS_PER_MONTH*
            Math.pow(10.0,n[i]+5.0)/(2.0*p[i]);

         // if that pass took more than 20000 milliseconds do fewer next time
        if(duration > 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<rhoreps;j++) {
               result = mpp.gcd(rhox.subtract(rhoy));
               rhox = rhox.modPow(onep, mpp).add(ONE).mod(mpp);
               rhox = rhox.modPow(onep, mpp).add(ONE).mod(mpp);
               rhoy = rhoy.modPow(onep, mpp).add(ONE).mod(mpp);
            }
            d = new Date(); finishtime = d.getTime();  // stop the clock
            duration = finishtime-starttime;
            rhotime = (double)duration/rhoreps/MILLISECONDS_PER_MONTH*
               2.4*Math.sqrt(Math.pow(10.0,n[i]+5.0))/Math.sqrt(p[i]);

            // accumulate data for geometric mean calculation
            rhocoeff = (double)duration/rhoreps/((double)p[i]*p[i]*Math.log(p[i]));
            rhocoeffproduct *= rhocoeff;

            // if that pass took more than 20000 milliseconds do fewer next time
            if(duration > 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));
      }
   }
}