mersenneforum.org Playing with Fermat
 User Name Remember Me? Password
 Register FAQ Search Today's Posts Mark Forums Read

 2022-09-28, 15:37 #1 jnml   Feb 2012 Prague, Czech Republ 2×32×11 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   Jul 2022 1516 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 4,373 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 4,373 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 19810 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: $\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.
 2022-10-08, 15:49 #6 paulunderwood     Sep 2002 Database er0rr 111516 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

2×32×11 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
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/

2022-10-08, 17:22   #8
paulunderwood

Sep 2002
Database er0rr

104258 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?
See this page: https://en.wikipedia.org/wiki/Euler_pseudoprime

Last fiddled with by paulunderwood on 2022-10-08 at 17:24

2022-10-08, 18:04   #9
jnml

Feb 2012
Prague, Czech Republ

2×32×11 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×3×127 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

 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 05:28.

Fri Dec 2 05:28:52 UTC 2022 up 106 days, 2:57, 0 users, load averages: 0.71, 0.90, 1.03

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.

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