Re: faster for very big primes and small exponents than "isprimepower()"

[Karim Belabas <Karim.Belabas@u-bordeaux.fr>]

Hi Hermann,

  for a fair comparison between your procedure and isprimepower, then
you must include a primality proof (isprime) in the former.

If the goal is to try and produce a base n given N = n^k without checking
whether n is prime, then benchmark your procedure against ispower.

Cheers,

    K.B.
-- 
Pr. Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/


* hermann@stamm-wilbrandt.de [2025-05-18 01:30]:
> On mersenne forum I learned about a method to determine prime of prime
> power:
> https://www.mersenneforum.org/node/1055690?p=1077586#post1077586
> 
> After I simplified modulus and exponent, it can be proven with Euler's
> theorem and induction that for odd prime powers N=p^n
> 2^N ≡ 2 (mod N)
> 
> This determines the prime for most odd prime powers:
> lift(gcd(Mod(2,N)^N-2,N))
> 
> The exceptions are the two Wiefrich primes below 18*10^18
> (1093 and 3511), for which the square of the prime is returned.
> 
> 
> PARI/GP "isprimepower(N,&r);return(r)" is nearly always much faster.
> I created script to demonstrate that the gcd is faster for very big primes
> and low exponents:
> https://www.mersenneforum.org/node/1055690?p=1077586#post1077586
> 
> hermann@7950x:~$ gp -q < <(curl --no-progress-meter https://gist.githubusercontent.com/Hermann-SW/e4c2ccccd918e59766dc82b01993d41b/raw/e7480819b8a57331d84cd30f6a9f89650a44388e/measure.gp)
> runtimes in milliseconds, with N=nextprime(base)^expo
> 
> isprimepower(N,&r);return(r)
>       |  2|  3|  4|  5|  6|  7|  8|  9| 10| 11| 12| 13| 14| 15|
>       |---|---|---|---|---|---|---|---|---|---|---|---|---|---|
> 10^ 10|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|
> 10^ 20|  0|  0|  0|  0|  0|  0|  0|  1|  0|  0|  0|  0|  0|  0|
> 10^ 30|  1|  1|  1|  1|  1|  1|  1|  1|  1|  1|  0|  0|  1|  1|
> 10^ 40|  2|  2|  2|  2|  2|  2|  1|  2|  2|  2|  2|  2|  1|  1|
> 10^ 50|  1|  0|  0|  1|  0|  0|  1|  0|  0|  1|  0|  0|  0|  0|
> 10^ 60|  5|  5|  5|  5|  5|  5|  5|  4|  4|  5|  5|  5|  5|  5|
> 10^ 70| 23| 23| 23| 23| 23| 23| 23| 23| 23| 23| 23| 23| 23| 23|
> 10^ 80| 13| 13| 13| 13| 12| 12| 12| 12| 12| 12| 12| 12| 12| 12|
> 10^ 90| 18| 18| 19| 19| 19| 19| 18| 18| 18| 17| 16| 16| 17| 16|
> 10^100| 23| 23| 24| 24| 24| 24| 23| 23| 23| 23| 23| 23| 23| 24|
> 10^110|  0|  1|  0|  1|  0|  1|  0|  1|  0|  1|  0|  1|  0|  1|
> 10^120| 51| 50| 51| 51| 51| 51| 51| 51| 51| 51| 51| 51| 51| 51|
> 10^130|109|109|108|108|109|109|108|109|108|109|108|108|109|109|
> 10^140|115|115|115|115|115|115|115|115|114|115|115|115|115|115|
> 10^150|136|136|136|136|136|136|136|136|136|136|136|136|136|136|
> 10^160|217|216|216|216|216|216|216|216|216|216|216|217|216|216|
> 
> lift(gcd(Mod(2,N)^N-2,N))
>       |  2|  3|  4|  5|  6|  7|  8|  9| 10| 11| 12| 13| 14| 15|
>       |---|---|---|---|---|---|---|---|---|---|---|---|---|---|
> 10^ 10|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|
> 10^ 20|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  1|  0|
> 10^ 30|  0|  0|  0|  0|  0|  0|  1|  0|  0|  1|  0|  1|  0|  1|
> 10^ 40|  0|  0|  0|  0|  0|  0|  1|  0|  1|  1|  1|  1|  2|  1|
> 10^ 50|  0|  0|  0|  0|  1|  0|  1|  1|  1|  2|  1|  2|  3|  3|
> 10^ 60|  0|  0|  0|  1|  0|  1|  1|  1|  2|  3|  3|  3|  5|  5|
> 10^ 70|  0|  0|  0|  1|  1|  1|  1|  2|  2|  4|  4|  6|  6|  8|
> 10^ 80|  0|  0|  0|  0|  1|  1|  3|  3|  4|  5|  6|  7|  9| 11|
> 10^ 90|  0|  0|  0|  1|  2|  2|  3|  4|  5|  7|  8| 10| 12| 15|
> 10^100|  0|  0|  1|  1|  1|  3|  4|  5|  6|  9| 11| 13| 16| 20|
> 10^110|  0|  0|  0|  2|  2|  4|  5|  7|  8| 10| 14| 18| 21| 26|
> 10^120|  0|  0|  1|  2|  3|  4|  6|  8| 11| 14| 18| 22| 28| 32|
> 10^130|  0|  1|  1|  2|  4|  5|  8| 10| 13| 17| 22| 28| 33| 41|
> 10^140|  0|  1|  2|  3|  4|  7|  9| 13| 16| 21| 27| 33| 40| 49|
> 10^150|  1|  1|  2|  3|  5|  8| 11| 15| 20| 26| 32| 39| 48| 55|
> 10^160|  0|  1|  2|  4|  6|  9| 13| 17| 24| 31| 38| 47| 54| 58|
> hermann@7950x:~$
> 
> 
> Regards,
> 
> Hermann.