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/
>
>
> * [email protected] [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.
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.