hermann on Sun, 18 May 2025 14:08:26 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

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

On 2025-05-18 03:24, Karim Belabas wrote:

Hi Hermann,

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


Hi Karim,
you are right, below runtimes are just for the isprime() test, not for lift(...). And those are roughly identical with the table shown previously for isprimepower().
Also there is a problem, while the result for N prime power is correct, what to do with this response that is prime but for non prime power: 

? N=3*5^2;lift(gcd(Mod(2,N)^N-2,N))
3
?


hermann@7950x:~$ gp -q < m.gp
      |  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|  0|  0|
10^ 30|  1|  1|  1|  1|  1|  1|  1|  1|  1|  1|  1|  1|  1|  0|
10^ 40|  2|  2|  1|  2|  1|  2|  2|  1|  2|  1|  2|  2|  1|  2|
10^ 50|  0|  0|  0|  0|  0|  0|  1|  0|  1|  0|  1|  0|  1|  0|
10^ 60|  5|  5|  5|  4|  5|  5|  5|  5|  5|  5|  5|  5|  5|  5|
10^ 70| 22| 23| 22| 23| 22| 23| 23| 23| 23| 22| 23| 22| 21| 21|
10^ 80| 12| 12| 11| 11| 11| 11| 11| 11| 11| 12| 12| 12| 12| 12|
10^ 90| 16| 17| 17| 16| 17| 17| 16| 17| 17| 16| 16| 17| 17| 16|
10^100| 24| 24| 24| 24| 23| 24| 24| 24| 24| 24| 24| 24| 24| 24|
10^110|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|  0|
10^120| 52| 52| 52| 52| 52| 52| 52| 52| 52| 52| 52| 52| 52| 52|
10^130|110|110|111|111|110|110|110|111|110|110|111|111|110|111|
10^140|116|117|117|117|117|116|117|117|117|116|116|117|117|117|
10^150|138|139|138|138|138|138|138|138|138|138|138|138|138|138|
10^160|220|220|220|221|220|220|220|221|221|221|221|221|221|222|
hermann@7950x:~$
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.