mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Miscellaneous Math

Reply
 
Thread Tools
Old 2022-09-28, 15:37   #1
jnml
 
Feb 2012
Prague, Czech Republ

2×32×11 Posts
Default 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.
?
jnml is offline   Reply With Quote
Old 2022-09-29, 13:57   #2
Neptune
 
Jul 2022

3×7 Posts
Default 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.
Neptune is offline   Reply With Quote
Old 2022-09-29, 17:29   #3
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

54·7 Posts
Default

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 is offline   Reply With Quote
Old 2022-09-30, 05:49   #4
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

54×7 Posts
Default

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
paulunderwood is offline   Reply With Quote
Old 2022-10-08, 10:30   #5
jnml
 
Feb 2012
Prague, Czech Republ

2×32×11 Posts
Default 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:

 \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.
jnml is offline   Reply With Quote
Old 2022-10-08, 15:49   #6
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

54×7 Posts
Default

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.
paulunderwood is offline   Reply With Quote
Old 2022-10-08, 16:16   #7
jnml
 
Feb 2012
Prague, Czech Republ

3068 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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?

Last fiddled with by jnml on 2022-10-08 at 16:30 Reason: s/that/than/ s/on/in/
jnml is offline   Reply With Quote
Old 2022-10-08, 17:22   #8
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

104278 Posts
Default

Quote:
Originally Posted by jnml View Post
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

Last fiddled with by paulunderwood on 2022-10-08 at 17:24
paulunderwood is offline   Reply With Quote
Old 2022-10-08, 18:04   #9
jnml
 
Feb 2012
Prague, Czech Republ

2·32·11 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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 ...".
jnml is offline   Reply With Quote
Old 2022-10-08, 20:39   #10
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

2×3,049 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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
Dr Sardonicus is offline   Reply With Quote
Reply

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

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


Fri Dec 2 20:51:21 UTC 2022 up 106 days, 18:19, 0 users, load averages: 1.08, 1.15, 1.14

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2022, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.

≠ ± ∓ ÷ × · − √ ‰ ⊗ ⊕ ⊖ ⊘ ⊙ ≤ ≥ ≦ ≧ ≨ ≩ ≺ ≻ ≼ ≽ ⊏ ⊐ ⊑ ⊒ ² ³ °
∠ ∟ ° ≅ ~ ‖ ⟂ ⫛
≡ ≜ ≈ ∝ ∞ ≪ ≫ ⌊⌋ ⌈⌉ ∘ ∏ ∐ ∑ ∧ ∨ ∩ ∪ ⨀ ⊕ ⊗ 𝖕 𝖖 𝖗 ⊲ ⊳
∅ ∖ ∁ ↦ ↣ ∩ ∪ ⊆ ⊂ ⊄ ⊊ ⊇ ⊃ ⊅ ⊋ ⊖ ∈ ∉ ∋ ∌ ℕ ℤ ℚ ℝ ℂ ℵ ℶ ℷ ℸ 𝓟
¬ ∨ ∧ ⊕ → ← ⇒ ⇐ ⇔ ∀ ∃ ∄ ∴ ∵ ⊤ ⊥ ⊢ ⊨ ⫤ ⊣ … ⋯ ⋮ ⋰ ⋱
∫ ∬ ∭ ∮ ∯ ∰ ∇ ∆ δ ∂ ℱ ℒ ℓ
𝛢𝛼 𝛣𝛽 𝛤𝛾 𝛥𝛿 𝛦𝜀𝜖 𝛧𝜁 𝛨𝜂 𝛩𝜃𝜗 𝛪𝜄 𝛫𝜅 𝛬𝜆 𝛭𝜇 𝛮𝜈 𝛯𝜉 𝛰𝜊 𝛱𝜋 𝛲𝜌 𝛴𝜎𝜍 𝛵𝜏 𝛶𝜐 𝛷𝜙𝜑 𝛸𝜒 𝛹𝜓 𝛺𝜔