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])))))
}
[/CODE]

Can you find one?

[QUOTE=paulunderwood;563539][CODE]
{
forstep(a=3,100,1,print(a);
<snip>
}
[/CODE][/QUOTE]
This may be a silly question, but did you perhaps mean forstep(a=3,100,[color=red]2[/color] ?

[QUOTE=paulunderwood;563539]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])))))
}
[/CODE]

Can you find one?[/QUOTE]

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]



[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, ...[/CODE]


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) & [B]Mod(a^4-8*a^2-12,n)^((n-1)/2)==(1)[/B], print([a,n])))))[/CODE]

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,...[/CODE]

[QUOTE=Dr Sardonicus;563549]This may be a silly question, but did you perhaps mean forstep(a=3,100,[color=red]2[/color] ?[/QUOTE]

No. "n" is the number tested, although I do not see how to make it an argument of a function.

[QUOTE=carpetpool;563552]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]



[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, ...[/CODE]


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) & [B]Mod(a^4-8*a^2-12,n)^((n-1)/2)==(1)[/B], print([a,n])))))[/CODE]

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,...[/CODE][/QUOTE]

I don't think your rearrangement of my code is the same!

[QUOTE=Dr Sardonicus;563549]This may be a silly question, but did you perhaps mean forstep(a=3,100,[color=red]2[/color] ?[/QUOTE]

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.

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

[QUOTE=paulunderwood;563563]No.[/QUOTE]
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 :sad:

[QUOTE=LaurV;563603]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 :sad:[/QUOTE]

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.

[QUOTE=paulunderwood;563606]
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.[/QUOTE]

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.

[QUOTE=paulunderwood;563566]I found this counterexample for a=4: n=68368998319.[/quote]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][/code]Gee, that's funny...


[quote]a even means gcd(D,E)>1.[/QUOTE]Yes, if "a" is even, gcd(D,E) is 4.