mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Miscellaneous Math (https://www.mersenneforum.org/forumdisplay.php?f=56)
-   -   1+1 Selfridges PRP test (https://www.mersenneforum.org/showthread.php?t=26205)

 paulunderwood 2020-11-17 19:07

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?

 Dr Sardonicus 2020-11-17 20:49

[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] ?

 carpetpool 2020-11-17 20:54

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

 paulunderwood 2020-11-17 22:16

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

 paulunderwood 2020-11-17 22:18

[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!

 paulunderwood 2020-11-17 22:43

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

 Dr Sardonicus 2020-11-17 23:54

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

 LaurV 2020-11-18 04:02

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

 paulunderwood 2020-11-18 04:12

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

 paulunderwood 2020-11-18 14:48

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

 Dr Sardonicus 2020-11-18 15:23

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

 paulunderwood 2020-11-18 20:17

What have the romans ever done for us?

[QUOTE=paulunderwood;563661]
After I will hammer a=3 to 10 to see if I can find another pseudoprime.
[/QUOTE]

[177655186463, 3]
[286748218763, 3]
[340473667463, 3]

I will update as they are found.

 Dr Sardonicus 2020-11-19 00:21

[QUOTE=paulunderwood;563692][177655186463, 3]
[286748218763, 3]
[340473667463, 3]

I will update as they are found.[/QUOTE]

[code]? factor(177655186463)
%1 =
[99347 1]

[1788229 1]

? for(i=1,2,print(factor(%1[i,1]-1)))
[2, 1; 13, 1; 3821, 1]
[2, 2; 3, 2; 13, 1; 3821, 1]

? factor(286748218763)
%2 =
[463747 1]

[618329 1]

? for(i=1,2,print(factor(%2[i,1]-1)))
[2, 1; 3, 1; 77291, 1]
[2, 3; 77291, 1]

? factor(340473667463)
%3 =
[505327 1]

[673769 1]

? for(i=1,2,print(factor(%3[i,1]-1)))
[2, 1; 3, 1; 84221, 1]
[2, 3; 84221, 1][/code]
Hmm. That's funny...

 paulunderwood 2020-11-19 08:06

Results for a=4:

[68368998319, 4]
[416032697359, 4]
[4752200398399, 4]

All semi-primes, p*q, where gcd(p-1,q-1)=2^i*3^j

 Dr Sardonicus 2020-11-19 14:51

[QUOTE=paulunderwood;563733]Results for a=4:

[68368998319, 4]
[416032697359, 4]
[4752200398399, 4]

All semi-primes, p*q, where gcd(p-1,q-1)=2^i*3^j[/QUOTE]

No. If g= gcd(p-1,q-1) it is not g, but (p-1)/g and (q-1)/g that are 2*i*3^j; in most of the cases g = p-1.

Furthermore, in every single case, p-1 and q-1 are small multiples of a single prime.

I think there's enough numerical evidence to try to devise a script to [i]construct[/i] counterexamples p*q to your test. For example, with a = 4, D = 12, E = 140, find primes l for which (say)

p = 6*l + 1 and q = 36*l + 1 are both prime

p*q is congruent to -1 (mod 560)

kronecker(12,p*q) = - 1

[there's another condition I'm too lazy to look up]

The congruence condition mod 560 produces something like 8 possibilities for l (mod 560).

The other examples show the multipliers 6 and 36 aren't the only possibilities.

 paulunderwood 2020-11-19 15:40

For a=5:

[17310451593311, 5]

gcd(p-1,q-1) = 10*17*4079

a=6 seems to be special case.

I will continue to a=10 and if no other counterexamples turn up, concentrate on a=7

 Kebbaj 2020-11-19 22:44

for a = 3 the numbers found all end with 63

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

for a = 3 the numbers found all end with 63:
[177655186463, 3]
[286748218763, 3]
[340473667463, 3]
so I added an "if ((n-floor (n / 100) * 100) == 63" , to find only the numbers that end with 63.

[CODE]
{
forstep(a=3,3,1,print(a);
D=a^2-4;E=D^2-4;
forstep(n=177655186463,1840473667452,4*E,
if((n-floor(n/100)*100)==63,
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]

There are no other numbers up to 1840473667452.

 paulunderwood 2020-11-20 02:20

[QUOTE=Kebbaj;563804]for a = 3 the numbers found all end with 63:
[177655186463, 3]
[286748218763, 3]
[340473667463, 3]
so I added an "if ((n-floor (n / 100) * 100) == 63" , to find only the numbers that end with 63.
[/QUOTE]

You could have done n%100==63 instead.

Anyway thanks.

With the counterexample [62611454518799, 9] I need to go back to the drawing board.

I had a PM for Dr. Sardonicus, wherein he describes a method to generate counterexamples for a=4. Perhaps he will post his method here.

 Kebbaj 2020-11-20 06:32

[QUOTE=paulunderwood;563807]You could have done n%100==63 instead.

Anyway thanks.

With the counterexample [62611454518799, 9] I need to go back to the drawing board.

I had a PM for Dr. Sardonicus, wherein he describes a method to generate counterexamples for a=4. Perhaps he will post his method here.[/QUOTE]

Thank you also to have drive me into this question of this exceptions. And also to put me on the language PARI, usually I program on MMA.
At the beginning I wanted to translate your code into MMA to better measure it, but I couldn't translate the "kronecker" function
(if someone knows how to translate the "kronecker" function into MMA?).
so I downloaded PARI 2.13, and i made my first code PARI.

Occasionally, I wanted to make a calculation status in your code, to know where "n" is in the calculation, but I didn't succeed:

[CODE]
{
forstep(a=3,3,1,print(a);D=a^2-4;E=D^2-4;s=187655186463;
forstep(n=177655186463,1840473667452,4*E,
if(n>=s,s=s+10000000000;print([n,s]),
if((n-floor(n%)*100)==63,
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]
I don't know why its plant?

 paulunderwood 2020-11-20 08:13

[QUOTE=Kebbaj;563818]
At the beginning I wanted to translate your code into MMA to better measure it, but I couldn't translate the "kronecker" function
(if someone knows how to translate the "kronecker" function into MMA?).
so I downloaded PARI 2.13, and i made my first code PARI.

Occasionally, I wanted to make a calculation status in your code, to know where "n" is in the calculation, but I didn't succeed:

I don't know why its plant?[/QUOTE]

bound=100000000000*E

I don't know MMA.

 Kebbaj 2020-11-20 11:01

OK that works!:smile:
I find that the functionality of PARI GP looks a lot of thinks like MMA.
[CODE]
{
str = fileopen("C:\\pari GP\\here.log", "a");
forstep(a=3,3,1,print(a);D=a^2-4;E=D^2-4;bound=177665186463;
forstep(n=177655186463,1840473667452,4*E,if((n-floor(n/100)*100)==63,
if(n>=bound,bound=100000000+bound;print([n,bound]),
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]);
filewrite(str, ([n,r,s]))))))));fileclose(str)
}
[/CODE]
Thanks

I will try to find for a = 7

 Dr Sardonicus 2020-11-20 13:16

[QUOTE=paulunderwood;563807]I had a PM for Dr. Sardonicus, wherein he describes a method to generate counterexamples for a=4. Perhaps he will post his method here.[/QUOTE]OK, you asked for it...

I gotta admit -- I have NO idea where most of the conditions in this "1 + 1 Selfridge test" come from.

I [i]do[/i] understand that given D = a^2 - 4, if kronecker(D,n) == -1 you want Mod(D,n)^((n-1)/2) == -1 also.

I do [i]not[/i] understand why E = D^2 - 4 comes into the picture.

I do [i]not[/i] understand where the restriction Mod(n,4*E) == -1 comes from (it's implicit in the loop that generates n-values).

I do [i]not[/i] understand where the condition Mod(E,n)^((n-1)/2) == 1 comes from.

These things notwithstanding, the composite examples n = p*q produced so far, as far as I have looked at them, are all very similar in one important respect: p-1 and q-1 have a very large common factor, most of which is its largest prime factor. In many (but not all) of these, gcd(p-1,q-1) = p-1. And in many (but not all) of the examples, q-1 = 6*(p-1).

So it seemed reasonably likely that by requiring p and q = 6*p - 5 both to be prime, and also requiring that p*q == -1 (mod 4*E), I would be able to generate additional examples for a-values for which such examples had already been found.

I also wondered "why 6?" Replacing 6 with a generic multiplier A led to an answer. Taking q-1 = A*(p-1), or q = A*p - (A-1), the condition p*q == -1 (mod 4*E) becomes

A*p^2 - (A-1)*p + 1 == 0 (mod n).

The discriminant of the quadratic is A^2 - 6*A + 1 = (A - 3)^2 - 8, which is easily seen to be a square for A = 6, but no other positive integer. For A = 6 alone, then, the quadratic factors, as

(3*p - 1)*(2*p - 1).

Now for a = 4, D = 12 and E = 140. We therefore want (3*p - 1)*(2*p - 1) == 0 (mod 16*5*7).

We get lucky with the prime-power factor 16, finding just one possibility, p == 11 (mod 16). We then find p == 2 or 3 (mod 5) and p == 4 or 5 (mod 7). We then find the 4 residue classes

p == 347, 187, 523, 123 (mod 560) make p*(6*p - 5) == -1 (mod 560).

So I cobbled together a rudimentary Pari-GP script and ran it for p%560 = 347, the class of the first example found. The output is a vector of vectors [p, q, n] with n = p*q which, if my code is writ right, will "fool" the test. The fact that the first [p, q, n] in the output agrees with the first example already found is encouraging. The script took seconds, not minutes, to run. I chose 60000000 as an upper bound because it was a bit below my (then) primelimit. If I'd bumped up my primelimit, I might have found more examples.

[code]? v=[];forprime(p=3,60000000,if(p%560==347&&ispseudoprime(6*p-5),q=6*p-5;n=p*q;if(kronecker(12,n)==-1&&Mod(12,n)^((n-1)/2)==Mod(-1,n)&&Mod(140,n)^((n-1)/2)==Mod(1,n),v=concat(v,[[p,q,n]]))))[/code]

The output was
[code][[106747, 640477, 68368998319], [1719547, 10317277, 17741042713519], [2588107, 15528637, 40189774120159], [2900587, 17403517, 50480415164479], [3335707, 20014237, 66761630460559], [3355867, 20135197, 67571043150799], [3784267, 22705597, 85924041442399], [6211867, 37271197, 231523718694799], [6257227, 37543357, 234917307091039], [6672187, 40033117, 267108442816879], [6956107, 41736637, 290324512792159], [8152267, 48913597, 398756702674399], [8234587, 49407517, 406850497190479], [10052347, 60314077, 606298030988719], [10111147, 60666877, 613411711377919], [10243867, 61463197, 629620815462799], [11124187, 66745117, 742485162844879], [12469867, 74819197, 932985435636799], [12587467, 75524797, 950665889919199], [12639547, 75837277, 958548826993519], [13862587, 83175517, 1153027840682479], [14017147, 84102877, 1178882390031919], [15888667, 95331997, 1514698354777999], [17778667, 106671997, 1896485912887999], [20263387, 121580317, 2463629014953679], [21765307, 130591837, 2842371423998959], [23364667, 140187997, 3275445867301999], [24918667, 149511997, 3725639665747999], [25530187, 153181117, 3910742561878879], [25921627, 155529757, 4031584348354639], [26067787, 156406717, 4077176984125279], [26397067, 158382397, 4180830745229599], [28201387, 169208317, 4771909231335679], [30707947, 184247677, 5657867900189119], [31067467, 186404797, 5791124879439199], [31904107, 191424637, 6107232101284159], [33431227, 200587357, 6705881465197039], [33459787, 200758717, 6717343909213279], [33955387, 203732317, 6917809668141679], [36312427, 217874557, 7911553946219839], [37080187, 222481117, 8249641422328879], [37753867, 226523197, 8552126651952799], [39835387, 239012317, 9521148145461679], [44816587, 268899517, 12051158597888479], [46059787, 276358717, 12729023640613279], [46877947, 281267677, 13185251255219119], [47847307, 287083837, 13736188483676959], [49532347, 297194077, 14720720148308719], [51180427, 307082557, 15716616391511839], [51425707, 308554237, 15867619785570559], [51639067, 309834397, 15999559185587599], [52756267, 316537597, 16699341982870399], [55047787, 330286717, 18181552846345279], [56956267, 341737597, 19464097818670399], [59467867, 356807197, 21218562935838799]][/code]

I ran the same script, substituting each of the other 3 possibilities. Only one of them, p%560 == 123, produced more examples:

[code]? v=[];forprime(p=3,60000000,if(p%560==123&&ispseudoprime(6*p-5),q=6*p-5;n=p*q;if(kronecker(12,n)==-1&&Mod(12,n)^((n-1)/2)==Mod(-1,n)&&Mod(140,n)^((n-1)/2)==Mod(1,n),v=concat(v,[[p,q,n]]))))[/code]
The output was
[code][[263323, 1579933, 416032697359], [889963, 5339773, 4752200398399], [1862683, 11176093, 20817518437519], [1872763, 11236573, 21043438161199], [5898043, 35388253, 208721437888879], [9625963, 57755773, 555954933934399], [11148043, 66888253, 745673120638879], [11683963, 70103773, 819089889892399], [12092203, 72553213, 877328179898239], [12140923, 72845533, 884412007046959], [12661723, 75970333, 961915312663759], [14240923, 85445533, 1216823256146959], [14509723, 87058333, 1263192296671759], [14931403, 89588413, 1337680698633439], [16752523, 100515133, 1683882077430559], [17448043, 104688253, 1826605139938879], [17972203, 107833213, 1938000394178239], [20416603, 122499613, 2501025966274639], [20418283, 122509693, 2501437581917119], [22852603, 137115613, 3133448668990639], [23942923, 143657533, 3439581250988959], [25438123, 152628733, 3882588483388159], [27729643, 166377853, 4613598466796479], [28734283, 172405693, 4953953973473119], [29226523, 175359133, 5125137733884559], [30051403, 180308413, 5418520783353439], [32339563, 194037373, 6275083848487999], [33041803, 198250813, 6550564307735839], [33248443, 199490653, 6632753605303279], [34155643, 204933853, 6999647521682479], [34224523, 205347133, 7027907676342559], [35385403, 212312413, 7512760295907439], [38004523, 228027133, 8666062420722559], [38212843, 229277053, 8761328029791679], [38545483, 231272893, 8914525365492319], [39372043, 236232253, 9300946423102879], [40364923, 242189533, 9775961850950959], [40400203, 242401213, 9793058212646239], [43120123, 258720733, 11156069829610159], [47189083, 283134493, 13360857090339919], [47436043, 284616253, 13501068815806879], [47567083, 285402493, 13575764072937919], [49957723, 299746333, 14974644274279759], [50997643, 305985853, 15604557294344479], [51464683, 308788093, 15891681320419519], [52985083, 317910493, 16844513858175919], [54053563, 324321373, 17530725767701999], [54626443, 327758653, 17904289375861279], [55204363, 331226173, 18285129889392799], [55610923, 333665533, 18555448263416959], [55642843, 333857053, 18576755584521679], [59261563, 355569373, 21071596798909999]][/code]
I'm sure there's a simple reason the other two residue classes didn't work, but I'm too lazy to track it down.

I tried this same setup with a = 6, D = 32, E = 1020 but found no examples. My script may have been writ wrong, of course, since I was doing this while up late, but it is also possible there is a simple reason this setup didn't work for a = 6.

I started to work out a generalization p = 2*a*l + 1, q = 2*b*l + 1 (where "l" will be a "large" prime). Proceeding as before, this leads to

2*a*b*l^2 + (a+b)*l + 1 == 0 (mod n)

The discriminant of the quadratic in l is a^2 - 6*a*b + b^2, a binary quadratic form with discriminant 32, which happens to be the same as D = 6^2 - 4.

There are of course values of a and b which make this form a perfect square; for these a and b the quadratic in l will factor, possibly leading to solutions. And that's as far as I've gotten.

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