FastECPP software and >50000 digit primality proof (reposted from NMBRTHRY)

[paulunderwood]

Originally Posted by R. Gerbicz

See https://github.com/gerbicz/CM/blob/main/corn.c .
The original code for Cornacchia just called Pari-Gp.
At 8k digits my Cornacchia code is 5.8 times faster, and for larger N it could be even faster.

Original code on approx 8k digits input (N=10^8007+21) using 4 physical cores, without gwnum library:

-- Size [0]: 26599 bits
   Time for discriminant -1642404: 470.9 (159.5)
   largest prime of d: 97
   largest prime of h: 5
       discriminants: 1.9 (1.9)
       54 qroot:      86.3 (29.3)
    57684 Cornacchia: 184.4 (62.4)
     1644 trial div:  46.0 (11.5)
       90 is_prime:   152.4 (54.4)
-- Size [1]: 26529 bits
   Time for discriminant -6487771: 1707.2 (575.7)
   largest prime of d: 907
   largest prime of h: 13
       discriminants: 2.0 (2.0)
      156 qroot:      255.6 (86.5)
   150518 Cornacchia: 479.9 (161.7)
     3236 trial div:  90.8 (22.8)
      816 is_prime:   1349.8 (462.3)

Using my own code:

-- Size [0]: 26599 bits
   Time for discriminant -1642404: 319.8 (108.5)
   largest prime of d: 97
   largest prime of h: 5
       discriminants: 1.9 (1.9)
       54 qroot:      86.4 (29.3)
    57684 Cornacchia: 32.9 (11.2)
     1644 trial div:  46.1 (11.6)
       90 is_prime:   152.5 (54.6)
-- Size [1]: 26529 bits
   Time for discriminant -6487771: 1459.2 (493.8)
   largest prime of d: 907
   largest prime of h: 13
       discriminants: 2.0 (2.0)
      156 qroot:      254.5 (86.2)
   150518 Cornacchia: 81.7 (27.7)
     3236 trial div:  91.0 (22.9)
      816 is_prime:   1349.8 (463.5)

Use gmp's half gcd code, do constant number of Euclidean algorithm's steps on the top limbs in constant time(!),
then do several mod tricks/quadratic residue tests to decrease the further computations.
More in the code, probably this time I have over commented. Used my own Legendre symbol code, it is also faster than the built-in kronecker. (calling it only for (k/n), when n is prime, so Legendre=kronecker).

Pointless to test it below 2^512, since it falls back to the original code. Also it is worth to see
the so called second step in the ecpp [I had code that failed in that part, since not saved V in the 1st step].

To compile: maybe it would work without any problem, for me the process:

On lib folder's pari.c replace the bool cm_pari_cornacchia function() with my code.
./configure --enable-mpi
On lib folder's makefile at the 263th line add to that line -I/home/gerbicz/gmp-6.2.1 -L/home/gerbicz/gmp-6.2.1
use your installed gmp's folder name. (this line was for the CFLAGS).
make install

Thanks Robert. The only thing I had to do differently to your instructions was to use

#include "/home/paul/Downloads/gmp-6.2.1/gmp-impl.h"

I am looking forward to seeing a nice speed up in production runs.

[Batalov]

Originally Posted by paulunderwood

Thanks Robert.

Thank you, Robert.
The two hack together cut down the run time by 2/3 in two categories: Cornacchia, is_prime.
The next contender for a 2/3-cut would be probably qroot?

...
       discriminants: 99.6 (99.8)
    17174 qroot:      561090.6 (19718.9)  <-- now the largest spender
 58575671 Cornacchia: 238202.6 (8039.3)
   380884 trial div:  314247.4 (9941.6)
    98055 is_prime:   522238.2 (18328.6)
-- Size [92]: 68343 bits

[R. Gerbicz]

Originally Posted by Batalov

The two hack together cut down the run time by 2/3 in two categories: Cornacchia, is_prime.
The next contender for a 2/3-cut would be probably qroot?

...
       discriminants: 99.6 (99.8)
    17174 qroot:      561090.6 (19718.9)  <-- now the largest spender
 58575671 Cornacchia: 238202.6 (8039.3)
   380884 trial div:  314247.4 (9941.6)
    98055 is_prime:   522238.2 (18328.6)
-- Size [92]: 68343 bits

Currently interested only in clean codes, but ouch, guessing that you still not observed, for dirty cracked code for Tonelli-Shanks:
on nt.c see line 268, 300, 308, you don't need more change, the others are trivial or expected O(1) number of mulmod, just use the gwnum library for powm. Paul could make the changes.

[paulunderwood]

Originally Posted by R. Gerbicz

Currently interested only in clean codes, but ouch, guessing that you still not observed, for dirty cracked code for Tonelli-Shanks:
on nt.c see line 268, 300, 308, you don't need more change, the others are trivial or expected O(1) number of mulmod, just use the gwnum library for powm. Paul could make the changes.

This is not trivial. Can b be negative in gmp_powm(r,b,e,n)?

I am finding it difficult to debug. If someone else want to have a crack at it then please go ahead.

[R. Gerbicz]

Originally Posted by paulunderwood

This is not trivial. Can b be negative in gmp_powm(r,b,e,n)?

It is always 0<b<n.
And on line 268 the b value is quite small (fits in 32 bits), for the other two lines it is just random number from [0,n).

Btw even in general it would not be a problem because mpz_mod(b,b,n) reduces b into the [0,n) interval, and then you could do the exponentiation.

[frmky]

Originally Posted by paulunderwood

Thanks Robert. The only thing I had to do differently to your instructions was to use
#include "/home/paul/Downloads/gmp-6.2.1/gmp-impl.h"

Can this be relatively easily rewritten with mpz_size(), mpz_limbs_modify(), and mpz_limbs_finish() to remove the dependency on gmp-impl.h, or is there something I'm missing in my quick review of the code?

[R. Gerbicz]

Originally Posted by frmky

Can this be relatively easily rewritten with mpz_size(), mpz_limbs_modify(), and mpz_limbs_finish() to remove the dependency on gmp-impl.h, or is there something I'm missing in my quick review of the code?

struct hgcd_matrix is defined in gmp-impl.h, also MPN_HGCD_MATRIX_INIT_ITCH(n) etc.
OK, then we could copy those definitions, but it could be not the same in different gmp versions and could change in the future...

[frmky]

Originally Posted by R. Gerbicz

struct hgcd_matrix is defined in gmp-impl.h, also MPN_HGCD_MATRIX_INIT_ITCH(n) etc.
OK, then we could copy those definitions, but it could be not the same in different gmp versions and could change in the future...

Thanks. I agree that copying those definitions isn't a good idea. I was hoping that it could be done using the documented public GMP interface relatively easily without loss of speed. Too much to ask perhaps!

[paulunderwood]

Originally Posted by paulunderwood

I am finding it difficult to debug. If someone else want to have a crack at it then please go ahead.

I have code that gets the right answer for GWNUM powm, but get an error "malloc(): corrupted top size". I have hit a brick wall. I use M4423 as a test number. Perhaps someone can debug it,

/*****************************************************************************/

unsigned int cm_nt_mpz_tonelli_generator (mpz_ptr q, mpz_ptr z,
   mpz_srcptr p)
   /* For p an odd prime, compute and return e such that p-1 = 2^e * q with
      q odd. If e>1 (that is, p != 3 mod 4), also q and a generator z of
      the 2-Sylow subgroup of (Z/pZ)^* are returned. Otherwise q and z
      are not modified. */
{

   unsigned int e = mpz_scan1 (p, 1);

   if (e > 1) {
      mpz_tdiv_q_2exp (q, p, e);
      for (mpz_set_ui (z, 2ul); mpz_legendre (z, p) != -1;
         mpz_add_ui (z, z, 1ul));
      if ( mpz_sizeinbase ( q, 2 ) < 120 )
         mpz_powm (z, z, q, p);
      else {

         #define SAFE 50

         char string [110000];
         int i, length;
         unsigned fft_size = 0;
         giant gz, gb, gp;
         gwnum wr, wb, wc;
         gwhandle *gwdata;

         length = ( strlen ( string ) >> 2 ) + 8;
         gp = newgiant ( length );
         mpz_get_str ( string, 10, p ); ctog ( string, gp );

         gwdata = (gwhandle*) malloc ( sizeof ( gwhandle ) );
         gwinit ( gwdata );
         gwset_larger_fftlen_count ( gwdata , fft_size );
         gwsetup_general_mod_giant ( gwdata, gp );
         wr = gwalloc ( gwdata );
         wb = gwalloc ( gwdata );
         wc = gwalloc ( gwdata );

         mpz_get_str ( string, 10, z );
         gb = newgiant ( length );
         ctog ( string, gb );
         gianttogw ( gwdata, gb, wb );

         int j = mpz_sizeinbase ( q, 2 ) - 2;
         int k = j - SAFE;
         gwcopy ( gwdata, wb, wr );

         for ( i = j; i > k; i-- ) {
            gwsquare2_carefully ( gwdata, wr, wr );
            if ( mpz_tstbit ( q, i ) ) gwmul ( gwdata, wb, wr );
         }

         for ( i = k; i > SAFE; i-- ) {
            gwcopy ( gwdata, wr, wc );
            gwmul ( gwdata, wr, wr );
            if ( mpz_tstbit ( q, i ) ) gwmul ( gwdata, wb, wr );
      /*
            if ( gw_get_maxerr ( gwdata ) < 0.44 ) { gwswap ( wr, wc );
            } else {
               printf ( "*** Maximum error %lf excessive at bit %d. Increasing FFT size.\n", gw_get_maxerr( gwdata ), i );
               gwtogiant ( gwdata, wr, gz );
               gwdone ( gwdata );
               ++fft_size;
               gwdata = (gwhandle*) malloc ( sizeof ( gwhandle ) );
               gwinit ( gwdata );
               gwset_larger_fftlen_count ( gwdata, fft_size );
               gwsetup_general_mod_giant ( gwdata, gp );
               wr = gwalloc ( gwdata );
               wb = gwalloc ( gwdata );
               wc = gwalloc ( gwdata );
               gwerror_checking ( gwdata, 1 );

               gianttogw ( gwdata, gz, wr );
               gianttogw ( gwdata, gb, wb );

               i++;
            }
*/
         }

         for ( i = SAFE; i >= 0; i-- ) {
            gwsquare2_carefully ( gwdata, wr, wr );
            if ( mpz_tstbit ( q, i ) ) gwmul ( gwdata, wb, wr );
         }

         gz = newgiant ( length );
         gwtogiant ( gwdata, wr, gz );
         gtompz( gz, z );
//mpz_out_str (NULL,10,z);printf("z\n");
         gwdone ( gwdata );
      }
   }

   return e;
}

/*****************************************************************************/

[paulunderwood]

I found the bug in my code. Attached is beta code. I think I have "length" right at 2^17 characters. Sorry for the ugliness of the code -- sheesh, two functions when one would suffice. The kick-in threshold is set at 4096 bits in 3 places. I shall inflict this code on my own production run.

(To install cm (ecpp-mpi) on Debian based systems follow EdH's instructions.)

Follow the instructions here and install Robert's code for some nice speed-ups.

Then just overwrite the two sections in lib/nt.c with the attached.

The beta2 file uses gwsmallmul rather than gwmul in the first section.

The beta3 file uses faster gwswap rather than gwcopy. It also presumes no round off error is thrown by gwsmallmul.

(beta and beta2 have been pulled because error checking was not turned on and because gwmul3_carefully was not there too -- thanks to forumite Luminescence for pointing out the errors.)

[paulunderwood]

(To install cm (ecpp-mpi) on Debian based systems follow EdH's instructions.)

After the instructions given here to get GWNUM working, install Robert's code for some nice speed-ups.

To cut along story short, attached is a C file that can be dropped into CM's lib directory. Then make the following changes to lib/nt.c:

#include <stdio.h>
#include "gw_utility5.c"
#include "cm-impl.h"

int cm_nt_is_prime (mpz_t a) {

   return ( gw_prp ( a ) );

}

and rename the three mpz_powm to gw_powm.

Also in nt/curve.c add #include gw_utility5.c and change mpz_powm to gw_powm.

Let me know if you have any trouble with the C code.

(Please re-download gw_utility5.c because I have added a mpz_clear.)