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

[email protected]
Newsgroups gmane.comp.mathematics.pari.user
Message-ID <[email protected]>
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.
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.