20220928, 15:37  #1 
Feb 2012
Prague, Czech Republ
3×67 Posts 
Playing with Fermat
Code:
? f(p)={a=Mod(2^((p+1)/2)1,2^p1);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. ? 
20220929, 13:57  #2 
"Martin Hopf"
Jul 2022
Germany
101010_{2} Posts 
Playing with LucasLehmer oneliner
Code:
ll(p)={my(a=4,m=2^p1,n=m+2,q=2*p);for(i=1,p2,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 20220929 at 14:38 Reason: code shortened by 5 characters. 
20220929, 17:29  #3 
Sep 2002
Database er0rr
5·29·31 Posts 
Code:
? ll(p)={my(a=4,m=2^p1,n=m+2,q=2*p);for(i=1,p2,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=mbitand(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 
20220930, 05:49  #4 
Sep 2002
Database er0rr
1000110001111_{2} 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 
20221008, 10:30  #5 
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 $ . 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. 
20221008, 15:49  #6 
Sep 2002
Database er0rr
10617_{8} Posts 
Writing a=2^((p+1)/2)1, you have the test a^(2^(p1))==a mod 2^p1 which is the same as a^(2^(p1)1)==1 mod 2^p1; The same as a^((2^p2)/2)==1 mod 2^p1. All that remains is for somebody show that jacobi(a,2^p1)==1 to show this is equivalent to a EulerPRP test with base a.

20221008, 16:16  #7  
Feb 2012
Prague, Czech Republ
3×67 Posts 
Quote:
"a^(2^(p1))==a mod 2^p1 which is the same as a^(2^(p1)1)==1 mod 2^p1" 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 about the Fermat primality test that describes these additional details? Last fiddled with by jnml on 20221008 at 16:30 Reason: s/that/than/ s/on/in/ 

20221008, 17:22  #8  
Sep 2002
Database er0rr
10617_{8} Posts 
Quote:
Quote:
Last fiddled with by paulunderwood on 20221008 at 17:24 

20221008, 18:04  #9  
Feb 2012
Prague, Czech Republ
3·67 Posts 
Quote:


20221008, 20:39  #10  
Feb 2017
Nowhere
1100001010000_{2} Posts 
Quote:
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 20221009 at 13:34 Reason: As indicated 

Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Just some fun playing with a Tesla P100 plus a question...  JonRussell  Hardware  9  20170427 11:46 
Playing with decimal representation  Nick  Puzzles  9  20130213 17:17 
Playing with WolframAlpha and musing.  Flatlander  Miscellaneous Math  12  20121129 09:56 
Playing with different radix  LoveCraft  Programming  7  20051114 07:59 
playing with numbers  michael  Puzzles  14  20040117 00:15 