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 2·33·67 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

24·32·31 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

2·163 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

1110001000102 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

1110001000102 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

2×33×67 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 10001011100002 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

Jun 2011
Thailand

33×347 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

361810 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

2×33×67 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

105608 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 02:22.

Fri Apr 16 02:22:02 UTC 2021 up 7 days, 21:02, 0 users, load averages: 1.26, 1.44, 1.58