1+1 Selfridges PRP test

[paulunderwood]

This script outputs no counterexamples:

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

[Dr Sardonicus]

Originally Posted by paulunderwood

{
forstep(a=3,100,1,print(a);
<snip>
}

This may be a silly question, but did you perhaps mean forstep(a=3,100,2 ?

[carpetpool]

Originally Posted by paulunderwood

This script outputs no counterexamples:

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

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

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:

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

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

[paulunderwood]

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.

[paulunderwood]

Originally Posted by carpetpool

There appear to be a handful of pseudoprimes without the BPSW requirement:

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

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:

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

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]

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.

[Dr Sardonicus]

The reason I asked about the forstep() was that it was a forstep() rather than a for()
;-)

[LaurV]

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

[paulunderwood]

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.

[paulunderwood]

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.

[Dr Sardonicus]

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.

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

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

Yes, if "a" is even, gcd(D,E) is 4.