mersenneforum.org 1+1 Selfridges PRP test
 User Name Remember Me? Password
 Register FAQ Search Today's Posts Mark Forums Read

 2020-11-17, 19:07 #1 paulunderwood     Sep 2002 Database er0rr 22·1,063 Posts 1+1 Selfridges PRP test This script outputs no counterexamples: Code: { forstep(a=3,100,1,print(a); D=a^2-4;E=D^2-4; forstep(n=4*E-1,50000000*E,4*E, if(kronecker(D,n)==-1&&!ispseudoprime(n), r=Mod(D,n)^((n-1)/2);s=Mod(E,n)^((n-1)/2); if(r==-1&&s==1, print([n,r,s]))))) } Can you find one?
2020-11-17, 20:49   #2
Dr Sardonicus

Feb 2017
Nowhere

5,897 Posts

Quote:
 Originally Posted by paulunderwood Code: { forstep(a=3,100,1,print(a); }
This may be a silly question, but did you perhaps mean forstep(a=3,100,2 ?

2020-11-17, 20:54   #3
carpetpool

"Sam"
Nov 2016

5×67 Posts

Quote:
 Originally Posted by paulunderwood This script outputs no counterexamples: Code: { forstep(a=3,100,1,print(a); D=a^2-4;E=D^2-4; forstep(n=4*E-1,50000000*E,4*E, if(kronecker(D,n)==-1&&!ispseudoprime(n), r=Mod(D,n)^((n-1)/2);s=Mod(E,n)^((n-1)/2); if(r==-1&&s==1, print([n,r,s]))))) } Can you find one?
There appear to be a handful of pseudoprimes without the BPSW requirement:

Code:
for(n=1,30000, if(n%2==1 & isprime(n)==0, for(a=1,n, if(kronecker(a^2-4,n)==(-1) & kronecker(a^4-8*a^2+12,n)==(-1) & Mod(a^2-4,n)^((n-1)/2)==(-1) & Mod(a^4-8*a^2-12,n)^((n-1)/2)==(-1), print([a,n])))))

Code:
703, 1387, 1891, 2071, 2743, 4187, 6943, 8227, 11359, 11773, 12403, 13019, 13747, 14383, 14701, 15251, 16441, 16531, 16589, 17261, 17767, 18721, 19093, 19951, 20567, 21667, 22411, 24727, 24929, 25531, 27217, 29341, 29891, ...

EDIT:

Code:
for(n=1,30000, if(n%2==1 & isprime(n)==0, for(a=1,n, if(kronecker(a^2-4,n)==(-1) & kronecker(a^4-8*a^2+12,n)==(-1) & Mod(a^2-4,n)^((n-1)/2)==(-1) & Mod(a^4-8*a^2-12,n)^((n-1)/2)==(1), print([a,n])))))
New list of pseudoprimes

Code:
259, 671, 703, 1649, 1891, 2059, 2701, 6001, 7471, 7957, 8911, 9211, 9881, 10963, 12403, 13213, 13747, 14111, 14491, 14701, 15203, 15251, 16531, 16589, 18631, 18721, 19951, 20263, 20567, 20591, 21349, 21667, 23377, 24461, 24727, 25351, 26599, 27613, 28939, 29341, 29539, 29891,...

Last fiddled with by carpetpool on 2020-11-17 at 21:03 Reason: s should be 1 not -1:

2020-11-17, 22:16   #4
paulunderwood

Sep 2002
Database er0rr

22·1,063 Posts

Quote:
 Originally Posted by Dr Sardonicus This may be a silly question, but did you perhaps mean forstep(a=3,100,2 ?
No. "n" is the number tested, although I do not see how to make it an argument of a function.

Last fiddled with by paulunderwood on 2020-11-17 at 22:23

2020-11-17, 22:18   #5
paulunderwood

Sep 2002
Database er0rr

10000100111002 Posts

Quote:
 Originally Posted by carpetpool There appear to be a handful of pseudoprimes without the BPSW requirement: Code: for(n=1,30000, if(n%2==1 & isprime(n)==0, for(a=1,n, if(kronecker(a^2-4,n)==(-1) & kronecker(a^4-8*a^2+12,n)==(-1) & Mod(a^2-4,n)^((n-1)/2)==(-1) & Mod(a^4-8*a^2-12,n)^((n-1)/2)==(-1), print([a,n]))))) Code: 703, 1387, 1891, 2071, 2743, 4187, 6943, 8227, 11359, 11773, 12403, 13019, 13747, 14383, 14701, 15251, 16441, 16531, 16589, 17261, 17767, 18721, 19093, 19951, 20567, 21667, 22411, 24727, 24929, 25531, 27217, 29341, 29891, ... EDIT: Code: for(n=1,30000, if(n%2==1 & isprime(n)==0, for(a=1,n, if(kronecker(a^2-4,n)==(-1) & kronecker(a^4-8*a^2+12,n)==(-1) & Mod(a^2-4,n)^((n-1)/2)==(-1) & Mod(a^4-8*a^2-12,n)^((n-1)/2)==(1), print([a,n]))))) New list of pseudoprimes Code: 259, 671, 703, 1649, 1891, 2059, 2701, 6001, 7471, 7957, 8911, 9211, 9881, 10963, 12403, 13213, 13747, 14111, 14491, 14701, 15203, 15251, 16531, 16589, 18631, 18721, 19951, 20263, 20567, 20591, 21349, 21667, 23377, 24461, 24727, 25351, 26599, 27613, 28939, 29341, 29539, 29891,...
I don't think your rearrangement of my code is the same!

2020-11-17, 22:43   #6
paulunderwood

Sep 2002
Database er0rr

22×1,063 Posts

Quote:
 Originally Posted by Dr Sardonicus This may be a silly question, but did you perhaps mean forstep(a=3,100,2 ?
On second thoughts, maybe I should! I found this counterexample for a=4: n=68368998319.

a even means gcd(D,E)>1.

Anyway I am re-running with some higher bounds.

Last fiddled with by paulunderwood on 2020-11-17 at 22:43

 2020-11-17, 23:54 #7 Dr Sardonicus     Feb 2017 Nowhere 170916 Posts The reason I asked about the forstep() was that it was a forstep() rather than a for() ;-)
2020-11-18, 04:02   #8
LaurV
Romulan Interpreter

"name field"
Jun 2011
Thailand

1001310 Posts

Quote:
 Originally Posted by paulunderwood No.
Then why the "forstep"? if the step is 1, you could just use a simple "for" which is a bit faster (internally, it uses increment instead of addition).
Edit: crosspost with Dr.S., sorry, I didn't refresh my page, from morning

Last fiddled with by LaurV on 2020-11-18 at 04:05

2020-11-18, 04:12   #9
paulunderwood

Sep 2002
Database er0rr

22·1,063 Posts

Quote:
 Originally Posted by LaurV Then why the "forstep"? if the step is 1, you could just use a simple "for" which is a bit faster (internally, it uses increment instead of addition). Edit: crosspost with Dr.S., sorry, I didn't refresh my page, from morning
It was a legacy contruction from hacking together the code, To be more efficient by far I should pre-calculate 50000000*E in the second loop!

Another speed up might be had be putting !ispseudoprime() ahead of kronecker().

Anyway, I am now running with the bound 2250000000*E and have reached a=41. The only counterexample is the aforementioned one for a=4.

Last fiddled with by paulunderwood on 2020-11-18 at 04:13

2020-11-18, 14:48   #10
paulunderwood

Sep 2002
Database er0rr

22·1,063 Posts

Quote:
 Originally Posted by paulunderwood Anyway, I am now running with the bound 2250000000*E and have reached a=41. The only counterexample is the aforementioned one for a=4.
I ran that up tp a=100, No new pseudoprimes.

I am now looping heavily over a with a shallow inner loop.

After I will hammer a=3 to 10 to see if I can find another pseudoprime.

It is remarkable that after 14 hours of computation only one pseudoprime was found.

Last fiddled with by paulunderwood on 2020-11-18 at 14:54

2020-11-18, 15:23   #11
Dr Sardonicus

Feb 2017
Nowhere

5,897 Posts

Quote:
 Originally Posted by paulunderwood I found this counterexample for a=4: n=68368998319.
I checked this example, just to see if anything interesting was going on.
Code:
? factor(68368998319)

%1 =
[106747 1]

[640477 1]

? for(i=1,#%1[,1],print(factor(%1[i,1]-1)))
[2, 1; 3, 1; 17791, 1]
[2, 2; 3, 2; 17791, 1]
Gee, that's funny...

Quote:
 a even means gcd(D,E)>1.
Yes, if "a" is even, gcd(D,E) is 4.

 Similar Threads Thread Thread Starter Forum Replies Last Post sweety439 sweety439 7 2020-02-11 19:49 Trilo Miscellaneous Math 25 2018-03-11 23:20 lidocorc Software 3 2008-12-03 15:12 swinster Software 2 2007-12-01 17:54 T.Rex Math 0 2004-10-26 21:37

All times are UTC. The time now is 15:36.

Sat Aug 13 15:36:52 UTC 2022 up 37 days, 10:24, 2 users, load averages: 1.64, 1.15, 1.00

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.

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