- **Miscellaneous Math**
(*https://www.mersenneforum.org/forumdisplay.php?f=56*)

- - **1+1 Selfridges PRP test**
(*https://www.mersenneforum.org/showthread.php?t=26205*)

1+1 Selfridges PRP testThis 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. |

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

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

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

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 |

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

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

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

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

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 |

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

Powered by vBulletin® Version 3.8.11

Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.