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

[hermann@stamm-wilbrandt.de]

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.