mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Miscellaneous Math

Reply
 
Thread Tools
Old 2020-11-17, 19:07   #1
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2×5×359 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?
paulunderwood is offline   Reply With Quote
Old 2020-11-17, 20:49   #2
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

10000110111002 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
Code:
{
forstep(a=3,100,1,print(a);
<snip>
}
This may be a silly question, but did you perhaps mean forstep(a=3,100,2 ?
Dr Sardonicus is offline   Reply With Quote
Old 2020-11-17, 20:54   #3
carpetpool
 
carpetpool's Avatar
 
"Sam"
Nov 2016

2×163 Posts
Post

Quote:
Originally Posted by paulunderwood View Post
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:
carpetpool is offline   Reply With Quote
Old 2020-11-17, 22:16   #4
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

359010 Posts
Default

Quote:
Originally Posted by Dr Sardonicus View Post
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
paulunderwood is offline   Reply With Quote
Old 2020-11-17, 22:18   #5
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2·5·359 Posts
Default

Quote:
Originally Posted by carpetpool View Post
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!
paulunderwood is offline   Reply With Quote
Old 2020-11-17, 22:43   #6
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2×5×359 Posts
Default

Quote:
Originally Posted by Dr Sardonicus View Post
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
paulunderwood is offline   Reply With Quote
Old 2020-11-17, 23:54   #7
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

10DC16 Posts
Default

The reason I asked about the forstep() was that it was a forstep() rather than a for()
;-)
Dr Sardonicus is offline   Reply With Quote
Old 2020-11-18, 04:02   #8
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

220618 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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
LaurV is offline   Reply With Quote
Old 2020-11-18, 04:12   #9
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2·5·359 Posts
Default

Quote:
Originally Posted by LaurV View Post
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
paulunderwood is offline   Reply With Quote
Old 2020-11-18, 14:48   #10
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2·5·359 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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
paulunderwood is offline   Reply With Quote
Old 2020-11-18, 15:23   #11
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

22×13×83 Posts
Default

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

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
I found the primality test, there seems to be no composite numbers that pass the test sweety439 sweety439 7 2020-02-11 19:49
Modifying the Lucas Lehmer Primality Test into a fast test of nothing Trilo Miscellaneous Math 25 2018-03-11 23:20
Double check LL test faster than first run test lidocorc Software 3 2008-12-03 15:12
Will the torture test, test ALL available memory? swinster Software 2 2007-12-01 17:54
A primality test for Fermat numbers faster than P├ępin's test ? T.Rex Math 0 2004-10-26 21:37

All times are UTC. The time now is 22:32.

Fri Feb 26 22:32:46 UTC 2021 up 85 days, 18:44, 0 users, load averages: 3.41, 2.90, 2.87

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, 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.