modulear exponentiation performance options: Mod(base, modulus)^exp / gp2c / mpz_powm() / gw_powm()

[hermann@stamm-wilbrandt.de]

Most times I use GP code these days. And when it comes to very big 
numbers, I use GP code for easy prototyping first.


In a recent thread I used GP for a 3,710 decimal digit Carmichael number 
from Dubner 1989 paper.

On mersenneforum.org (you need account to view anything) there was a 
related discussion (on factoring N-1 for that Carmichael number) and I 
learned that @Neptune had determined much bigger 3-Carmichael number 
(with 3 prime factors) of 96,522 decimal digits.

I created 96522.gp gist for that number:
https://gist.github.com/Hermann-SW/2d8dbe86f105bd0b5bfdacf51ed98f3e

It computes the Carmichael number, outputs number of decimal digits and 
number of decimal digits of composite computed from N-1, P-1, Q-1 and 
R-1 (N=P*Q*R). Finally for random base a lift(Mod(a, N)^(N-1)) is 
computed and sent to output (for Carmichael numbers that is always 1).

It took 8:05min:
hermann@7950x:~$ time gp -q < 96522.gp
96522
48267
1
took 481803ms to compute lift(Mod(2+random(),n)^(n-1))

real	8m5.826s
user	8m1.829s
sys	0m3.969s
hermann@7950x:~$


Using gp2c (to compile slightly modified GP script, gp2c does not like 
"eval()") the runtime reduced to 7:52min.


Next I transpiled 96522.gp to C++ using libgmp, libgmpxx and mpz_powm():
https://github.com/Hermann-SW/prp_and_sqrtm1/blob/main/96522.cc

Using mpz_powm() runtime reduced a bit further to 6:32min().


Finally I used @paulunderwood's gw_utility (based on highly efficient 
gwnum lib from GIMPS project) gw_powm() dropin replacement and now same  
algorithm needs only 35 seconds (working on 96,522 decimal digit 
modulus):
hermann@7950x:~/prp_and_sqrtm1$ rm a.out
hermann@7950x:~/prp_and_sqrtm1$ LIBS="-lgmp -lgmpxx -Lp95v3019b14/gwnum 
-l:gwnum.a"
hermann@7950x:~/prp_and_sqrtm1$ g++ -O3 96522.cc 
paulunderwood/gw_utility.o $LIBS 2>err
hermann@7950x:~/prp_and_sqrtm1$ time ./a.out
1

real	0m35.165s
user	0m35.545s
sys	0m0.016s
hermann@7950x:~/prp_and_sqrtm1$


For the 3,710 decimal digit Carmichael number there is no real need to 
go the gw_powm() route. But for 96,522 decimal number (or bigger) it is 
good for me to have the gw_powm() option.


Regards,

Hermann.