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

[paulunderwood]

Originally Posted by Luminescence

fft size: 16384
avx: 8
round off 3: 0.000000
round off 2: 0.000000
prp

real	0m9,978s

Nice.

So in short: the latest code works, but the FFT size needs to be increased.

I have just put out (post 148) code that does autoincrement of FFT size on excessive round-off error. Attached is the concept test program.

[frmky]

*** Maximum error of 0.500000 excessive at iteration 108154. Increasing FFT size.
prp

real    0m41.429s

Being a bit more conservative with a limit of 0.35 gives

*** Maximum error of 0.375000 excessive at iteration 108660. Increasing FFT size.
prp

real    0m41.536s

which switches a bit more than 500 iterations earlier.

[paulunderwood]

Attached is a version that uses a quick gwswap() rather than a slow gwcopy(). I have added the corresponding Ver_5 code to post 148.

[paulunderwood]

The ultimate Ver_6 has been uploaded to post 148. It uses functions and a neat piece of code to calculate odd d for n-1=d*2^r thanks to frmky. It might be a little bit quicker. As I state there, comment out the printf(); function messaging about FFT increases if you find them annoying at runtime.

frmky is working on 3-SPRP GMP code to be used above a certain threshold based this UTM page. This is useful for ARM based computers.

[frmky]

Originally Posted by paulunderwood

frmky is working on 3-SPRP GMP code to be used above a certain threshold based this UTM page. This is useful for ARM based computers.

Below is the code I settled on. I thank Paul for his suggestions and comments! I have tested it on many known primes and base-3 strong pseudoprimes.

This code is slower than Paul's version using gwnum but is faster than mpz_probab_prime_p(). It's useful where gwnum is not available.

int cm_nt_is_prime (mpz_t n)

{
   unsigned long j, r;
   int prime = 0;
   mpz_t d, a, x, n_sub_1;

   if (mpz_sizeinbase(n, 2) < 4096) return (mpz_probab_prime_p(n, 0) > 0);
   if (mpz_even_p(n)) return 0;

   mpz_inits(d, a, x, n_sub_1, NULL);

   mpz_sub_ui(n_sub_1, n, 1);
   r = mpz_scan1(n_sub_1, 0);
   mpz_fdiv_q_2exp(d, n_sub_1, r);

   mpz_set_ui(a, 3);
   mpz_powm(x, a, d, n);

   if ((mpz_cmp_ui(x, 1) == 0) || (mpz_cmp(x, n_sub_1) == 0)) prime = 1;
   else {
      for (j = 1; j < r; j++) {
         mpz_powm_ui(x, x, 2, n);
         if (mpz_cmp(x, n_sub_1) == 0) {
            prime = 1;
            break;
         }
      }
   }
    
   mpz_clears(d, a, x, n_sub_1, NULL);
   return prime;
}