mersenneforum.org  

Go Back   mersenneforum.org > Prime Search Projects > And now for something completely different

Reply
 
Thread Tools
Old 2022-11-27, 19:51   #210
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

23×3×11×17 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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:
Code:
-- 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:
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.

Last fiddled with by paulunderwood on 2022-11-27 at 20:59
paulunderwood is offline   Reply With Quote
Old 2022-11-29, 01:04   #211
Batalov
 
Batalov's Avatar
 
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

32×1,117 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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?
Code:
...
       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
Batalov is offline   Reply With Quote
Old 2022-11-29, 01:57   #212
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

32·179 Posts
Default

Quote:
Originally Posted by Batalov View Post
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?
Code:
...
       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.
R. Gerbicz is offline   Reply With Quote
Old 2022-11-29, 14:25   #213
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

23·3·11·17 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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.

Last fiddled with by paulunderwood on 2022-11-29 at 14:26
paulunderwood is offline   Reply With Quote
Old 2022-11-29, 14:38   #214
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

161110 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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.

Last fiddled with by R. Gerbicz on 2022-11-29 at 14:40
R. Gerbicz is offline   Reply With Quote
Old 2022-11-29, 20:20   #215
frmky
 
frmky's Avatar
 
Jul 2003
So Cal

72·53 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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?
frmky is offline   Reply With Quote
Old 2022-11-29, 21:16   #216
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

31138 Posts
Default

Quote:
Originally Posted by frmky View Post
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...
R. Gerbicz is offline   Reply With Quote
Old 2022-11-29, 21:35   #217
frmky
 
frmky's Avatar
 
Jul 2003
So Cal

72×53 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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!
frmky is offline   Reply With Quote
Old 2022-11-30, 06:28   #218
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

23·3·11·17 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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,



Code:
/*****************************************************************************/

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;
}

/*****************************************************************************/
Attached Files
File Type: txt GWNUM_hack1.txt (3.2 KB, 11 views)

Last fiddled with by paulunderwood on 2022-11-30 at 07:25
paulunderwood is offline   Reply With Quote
Old 2022-12-01, 23:40   #219
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

10001100010002 Posts
Arrow tonelli beta

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.)
Attached Files
File Type: txt tonelli_generators_beta3.txt (7.6 KB, 9 views)

Last fiddled with by paulunderwood on 2022-12-03 at 03:54
paulunderwood is offline   Reply With Quote
Old 2022-12-04, 03:17   #220
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

448810 Posts
Thumbs up New working code



(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:

Code:
#include <stdio.h>
#include "gw_utility5.c"
#include "cm-impl.h"
Code:
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.)
Attached Files
File Type: c gw_utility5.c (10.0 KB, 4 views)

Last fiddled with by paulunderwood on 2022-12-08 at 19:48
paulunderwood is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
what is the best primality software? bbb120 Proth Prime Search 75 2022-10-08 11:45
Fastest software for Mersenne primality test? JonathanM Information & Answers 25 2020-06-16 02:47
APR-CL as primality proof f1pokerspeed FactorDB 14 2014-01-09 21:06
Proof of Primality Test for Fermat Numbers princeps Math 15 2012-04-02 21:49
PRIMALITY PROOF for Wagstaff numbers! AntonVrba Math 96 2009-02-25 10:37

All times are UTC. The time now is 04:16.


Tue Jan 31 04:16:28 UTC 2023 up 166 days, 1:45, 0 users, load averages: 1.19, 0.98, 0.93

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2023, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.

≠ ± ∓ ÷ × · − √ ‰ ⊗ ⊕ ⊖ ⊘ ⊙ ≤ ≥ ≦ ≧ ≨ ≩ ≺ ≻ ≼ ≽ ⊏ ⊐ ⊑ ⊒ ² ³ °
∠ ∟ ° ≅ ~ ‖ ⟂ ⫛
≡ ≜ ≈ ∝ ∞ ≪ ≫ ⌊⌋ ⌈⌉ ∘ ∏ ∐ ∑ ∧ ∨ ∩ ∪ ⨀ ⊕ ⊗ 𝖕 𝖖 𝖗 ⊲ ⊳
∅ ∖ ∁ ↦ ↣ ∩ ∪ ⊆ ⊂ ⊄ ⊊ ⊇ ⊃ ⊅ ⊋ ⊖ ∈ ∉ ∋ ∌ ℕ ℤ ℚ ℝ ℂ ℵ ℶ ℷ ℸ 𝓟
¬ ∨ ∧ ⊕ → ← ⇒ ⇐ ⇔ ∀ ∃ ∄ ∴ ∵ ⊤ ⊥ ⊢ ⊨ ⫤ ⊣ … ⋯ ⋮ ⋰ ⋱
∫ ∬ ∭ ∮ ∯ ∰ ∇ ∆ δ ∂ ℱ ℒ ℓ
𝛢𝛼 𝛣𝛽 𝛤𝛾 𝛥𝛿 𝛦𝜀𝜖 𝛧𝜁 𝛨𝜂 𝛩𝜃𝜗 𝛪𝜄 𝛫𝜅 𝛬𝜆 𝛭𝜇 𝛮𝜈 𝛯𝜉 𝛰𝜊 𝛱𝜋 𝛲𝜌 𝛴𝜎𝜍 𝛵𝜏 𝛶𝜐 𝛷𝜙𝜑 𝛸𝜒 𝛹𝜓 𝛺𝜔