 mersenneforum.org Playing with Fermat
 Register FAQ Search Today's Posts Mark Forums Read 2022-09-28, 15:37 #1 jnml   Feb 2012 Prague, Czech Republ 3·67 Posts Playing with Fermat Code: ? f(p)={a=Mod(2^((p+1)/2)-1,2^p-1);b=a;for(i=2,p,a=sqr(a));b==a}; ? forprime(p=5,4423,if(f(p),print(p))) 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 ? ## *** last result: cpu time 3,185 ms, real time 3,186 ms. ?   2022-09-29, 13:57 #2 Neptune   "Martin Hopf" Jul 2022 Germany 1010102 Posts Playing with Lucas-Lehmer one-liner Code: ll(p)={my(a=4,m=2^p-1,n=m+2,q=2*p);for(i=1,p-2,if(m%(i*q+1),a*=a;a=bitand(a,m)+a>>p;if(a>=m,a-=n,a-=2),break));!a}; Code: ? forprime(p=7,4423,if(ll(p),print(p))) 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 ? ## *** last result: cpu time 782 ms, real time 783 ms. Code: ? parforprime(p=7,11213,if(ll(p),print(p))) 7 13 19 17 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 9689 9941 11213 ? ## *** last result: cpu time 26,105 ms, real time 3,371 ms. Last fiddled with by Neptune on 2022-09-29 at 14:38 Reason: code shortened by 5 characters.   2022-09-29, 17:29 #3 paulunderwood   Sep 2002 Database er0rr 5·29·31 Posts Code: ? ll(p)={my(a=4,m=2^p-1,n=m+2,q=2*p);for(i=1,p-2,if(m%(i*q+1),a*=a;a=bitand(a,m)+a>>p;if(a>=m,a-=n,a-=2),break));!a}; ? gettime();forprime(p=7,7000,if(ll(p),print(p)));gettime()7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 9420 Code: ? {flt(p)=my(s=3,m=(2^p)-1,d=2*p,q=1);l=round(2*p/log(p));for(k=1,l,q+=d;if(Mod(2,q)^p==1,return(0)));for(k=2,p,s=sqr(s);s=m-bitand(s,m)-s>>p);s==3;} ? gettime();forprime(p=7,7000,if(flt(p),print(p)));gettime() 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 7572   2022-09-30, 05:49 #4 paulunderwood   Sep 2002 Database er0rr 5·29·31 Posts Slighty quicker using the bitxor function: Code: ? {flt2(p)=my(s=3,m=(2^p)-1,d=2*p,q=1);l=round(2*p/log(p));for(k=1,l,q+=d;if(Mod(2,q)^p==1,return(0)));for(k=2,p,s=sqr(s);s=bitxor(bitand(s,m),m)-s>>p);s==3;} ? gettime();forprime(p=7,7000,if(flt2(p),print(p)));gettime() 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 7523   2022-10-08, 10:30 #5 jnml   Feb 2012 Prague, Czech Republ 3×67 Posts Why it works? The OP code seems to have no false negatives and no false positives for p <= 86243 (using a compiled language but no FFT multiplication bellow 110 kbits): Code: $grep PRIME nohup.out 11:51:03 3.88µs: M5, 776ns/iteration PRIME 11:51:03 2.93µs: M7, 418ns/iteration PRIME 11:51:03 3.93µs: M13, 302ns/iteration PRIME 11:51:03 5.14µs: M17, 302ns/iteration PRIME 11:51:03 5.64µs: M19, 296ns/iteration PRIME 11:51:03 8.9µs: M31, 287ns/iteration PRIME 11:51:03 17.28µs: M61, 283ns/iteration PRIME 11:51:03 33.4µs: M89, 375ns/iteration PRIME 11:51:03 38.78µs: M107, 362ns/iteration PRIME 11:51:03 47.04µs: M127, 370ns/iteration PRIME 11:51:03 193.129µs: M521, 370ns/iteration PRIME 11:51:03 182.22µs: M607, 300ns/iteration PRIME 11:51:03 1.232338ms: M1279, 963ns/iteration PRIME 11:51:03 3.693433ms: M2203, 1.676µs/iteration PRIME 11:51:04 4.340422ms: M2281, 1.902µs/iteration PRIME 11:51:04 8.920693ms: M3217, 2.772µs/iteration PRIME 11:51:06 12.958735ms: M4253, 3.046µs/iteration PRIME 11:51:06 19.731632ms: M4423, 4.461µs/iteration PRIME 11:51:42 134.703436ms: M9689, 13.902µs/iteration PRIME 11:51:46 142.715181ms: M9941, 14.356µs/iteration PRIME 11:52:07 201.56321ms: M11213, 17.975µs/iteration PRIME 11:58:49 823.044776ms: M19937, 41.282µs/iteration PRIME 12:01:37 995.261438ms: M21701, 45.862µs/iteration PRIME 12:04:33 1.282974407s: M23209, 55.279µs/iteration PRIME 13:48:43 5.198145861s: M44497, 116.82µs/iteration PRIME 08:39:40 35.605965118s: M86243, 412.856µs/iteration PRIME$ The bit more readable form of the test is: . For example: etc. It looks much like a Fermat test, but I'm not able to derive it from . Can anybody enlighten me please? Thanks.   2022-10-08, 15:49 #6 paulunderwood   Sep 2002 Database er0rr 10001100011112 Posts Writing a=2^((p+1)/2)-1, you have the test a^(2^(p-1))==a mod 2^p-1 which is the same as a^(2^(p-1)-1)==1 mod 2^p-1; The same as a^((2^p-2)/2)==1 mod 2^p-1. All that remains is for somebody show that jacobi(a,2^p-1)==1 to show this is equivalent to a Euler-PRP test with base a.   2022-10-08, 16:16   #7
jnml

Feb 2012
Prague, Czech Republ

3·67 Posts Quote:
 Originally Posted by paulunderwood Writing a=2^((p+1)/2)-1, you have the test a^(2^(p-1))==a mod 2^p-1 which is the same as a^(2^(p-1)-1)==1 mod 2^p-1; The same as a^((2^p-2)/2)==1 mod 2^p-1. All that remains is for somebody show that jacobi(a,2^p-1)==1 to show this is equivalent to a Euler-PRP test with base a.
Thanks for the answer! Not sure I understand this part:

"a^(2^(p-1))==a mod 2^p-1 which is the same as a^(2^(p-1)-1)==1 mod 2^p-1"

Plugin in, for example p = 7 I get

"a^64==a mod 127 which is the same as a^63==1 mod 127".

But Fermat tests is about the exponent on the left being one less than the modulo value on the right, not 64 like in this case. At least
that's all I can see in the Wikipedia article. Neither is 'jacobi' mentioned on that page. Is there a better, more informative resource

Last fiddled with by jnml on 2022-10-08 at 16:30 Reason: s/that/than/ s/on/in/   2022-10-08, 17:22   #8
paulunderwood

Sep 2002
Database er0rr

5×29×31 Posts Quote:
 Originally Posted by jnml Thanks for the answer! Not sure I understand this part: "a^(2^(p-1))==a mod 2^p-1 which is the same as a^(2^(p-1)-1)==1 mod 2^p-1" Plugin in, for example p = 7 I get "a^64==a mod 127 which is the same as a^63==1 mod 127".
Yes. You have divided both side of the equivalence by a. Note that 63=(127-1)/2.
Quote:
 But Fermat tests is about the exponent on the left being one less than the modulo value on the right, not 64 like in this case. At least that's all I can see in the Wikipedia article. Neither is 'jacobi' mentioned on that page. Is there a better, more informative resource about the Fermat primality test that describes these additional details?

Last fiddled with by paulunderwood on 2022-10-08 at 17:24   2022-10-08, 18:04   #9
jnml

Feb 2012
Prague, Czech Republ

3×67 Posts Quote:
 Originally Posted by paulunderwood Yes. You have divided both side of the equivalence by a. Note that 63=(127-1)/2. See this page: https://en.wikipedia.org/wiki/Euler_pseudoprime
Thanks for the link, I think I understand now. The key thing I was missing is in the sentence "This can be factored as ...".   2022-10-08, 20:39   #10
Dr Sardonicus

Feb 2017
Nowhere

24×389 Posts Quote:
 Originally Posted by paulunderwood Writing a=2^((p+1)/2)-1, you have the test a^(2^(p-1))==a mod 2^p-1 which is the same as a^(2^(p-1)-1)==1 mod 2^p-1; The same as a^((2^p-2)/2)==1 mod 2^p-1. All that remains is for somebody show that jacobi(a,2^p-1)==1 to show this is equivalent to a Euler-PRP test with base a.
For odd p > 1 we have ( 2^((p+1)/2) - 1, 2^p - 1) = -(2^p - 1, 2^((p+1)/2) - 1).

If p ≥ 5 then (p+1)/2 ≥ 3, so that (2, 2^((p+1)/2) - 1) = +1. Then

(2^p - 1, 2^((p+1)/2) - 1) = (2^(p+1) - 2, 2^((p+1)/2) - 1) = (-1, 2^((p+1)/2) - 1) = -1.

Thus ( 2^((p+1)/2) - 1, 2^p - 1) = -(-1) = +1.

EDIT: IIRC this very question came up some time ago on the Forum, and a proof was indicated at that time.

Unfortunately I was unable to find it or recall the proof, so I would up doing it over from scratch. I think it's pretty much the same proof.

Last fiddled with by Dr Sardonicus on 2022-10-09 at 13:34 Reason: As indicated  Thread Tools Show Printable Version Email this Page Similar Threads Thread Thread Starter Forum Replies Last Post JonRussell Hardware 9 2017-04-27 11:46 Nick Puzzles 9 2013-02-13 17:17 Flatlander Miscellaneous Math 12 2012-11-29 09:56 LoveCraft Programming 7 2005-11-14 07:59 michael Puzzles 14 2004-01-17 00:15

All times are UTC. The time now is 16:20.

Sat Feb 4 16:20:16 UTC 2023 up 170 days, 13:48, 1 user, load averages: 1.32, 0.99, 0.87