Playing with Fermat

jnml · 2022-09-28, 15:37

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.
?

Neptune · 2022-09-29, 13:57

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.

paulunderwood · 2022-09-29, 17:29

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

paulunderwood · 2022-09-30, 05:49

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

jnml · 2022-10-08, 10:30

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:

$\left( 2^{\frac{p+1}{2}}-1 \right) ^{2^{p-1}} \equiv 2^{p+1 \over 2}-1 \ (mod \ 2^p-1)$ .

For example:

$7^{16} \equiv 7 \ (mod \ 31)$

$15^{64} \equiv 15 \ (mod \ 127)$

$63^{1024} \equiv 1202 \ (mod \ 2047)$

$127^{4096} \equiv 127 \ (mod \ 8191)$

etc.

It looks much like a Fermat test, but I'm not able to derive it from $a^{x-1} \equiv 1 \ (mod \ x)$ . Can anybody enlighten me please? Thanks.

paulunderwood · 2022-10-08, 15:49

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.

jnml · 2022-10-08, 16:16

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
about the Fermat primality test that describes these additional details?

paulunderwood · 2022-10-08, 17:22

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?

See this page: https://en.wikipedia.org/wiki/Euler_pseudoprime

jnml · 2022-10-08, 18:04

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 ...".

Dr Sardonicus · 2022-10-08, 20:39

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.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Just some fun playing with a Tesla P100 plus a question...	JonRussell	Hardware	9	2017-04-27 11:46
Playing with decimal representation	Nick	Puzzles	9	2013-02-13 17:17
Playing with WolframAlpha and musing.	Flatlander	Miscellaneous Math	12	2012-11-29 09:56
Playing with different radix	LoveCraft	Programming	7	2005-11-14 07:59
playing with numbers	michael	Puzzles	14	2004-01-17 00:15

2022-10-08, 15:49	#6
paulunderwood Sep 2002 Database er0rr 2³·3·11·17 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.