Saving tons of PRP-C and (hypotetical) LL-C tests, a new method
 

PRP-C refers to (weak) Fermat test for (2^p-1)/d.
And LL-C refers to a LL type test for (2^p-1)/d, see [url]http://mersenneforum.org/showthread.php?p=471732#post471732[/url]
ofcourse currently the main issue with LL/LL-C is that we have no strong error check for these LL (type) tests.

But first concentrating on PRP-C tests.

You can see such lines when you can check out the recent results page:
[url]https://www.mersenne.org/report_exponent/?exp_lo=1344943&full=1[/url]
[url]https://www.mersenne.org/report_exponent/?exp_lo=2389403&full=1[/url]
[url]https://www.mersenne.org/report_exponent/?exp_lo=2462281&full=1[/url]
[url]https://www.mersenne.org/report_exponent/?exp_lo=3009703&full=1[/url]
[url]https://www.mersenne.org/report_exponent/?exp_lo=4373911&full=1[/url]

What happens here? We do a PRP-C test on (2^p-1)/d1, but when a new
prime divisor appears, then we redo this, because the full residue is not available.
Ofcourse there is a folklore solution: store the full residue. But that takes lots of storage and bandwidth.
I can do much better here, and this is the trick.

Let p>3 prime and mp=2^p-1.
Compute
[CODE]
[1]: 3^mp==res mod mp, independently from d.
[/CODE]
[in practice we compute 3^(mp+1), because then mp+1 is a power of two; and the end we use an easy division by 3].
If mp/d is prime, then using Fermat's theorem we know that:
3^(mp/d)==3 mod (mp/d), so
[CODE]
[2]: 3^mp==3^d mod (mp/d) also holds (note that here we're already weaker than
Fermat, but not much weaker).
[/CODE]
From [1] and [2]:
3^mp==res==3^d mod (mp/d)
let s=3^d mod (mp/d), then
[CODE]
[3]: res=s+w*(mp/d), for some w integer, where 0<=w<d. Multiply by d, we get
d*(res-s)=w*mp, see this equation mod 2^t, where t<p, then
[/CODE]

d*(res-s)==-w mod 2^t, so
w==d*(s-res) mod 2^t, but if 2^t>d, then we have a unique candidate for w:

[CODE]
[4]: w=(d*(s-res) mod 2^t)
[/CODE]
but if this is w>=d, then mp/d is composite!!! (because we had 0<=w<d)
(not forget the only conditions: t<p and d<2^t)


It was a quite general idea, works for all numbers, just go back to [3],
suppose N is odd, and d is a non-trivial divisor of N.
[3]': res=s+w*(N/d), see this mod 2^t
w==(res-s)*inv mod 2^t, where (N/d)*inv==1 mod 2^t, here inv exists, because N/d is also odd.
and if d<2^t and d<=w then N/d is composite.
[ where w=((res-s)*inv) mod 2^t) ]. Not get confused, but in Mersenne number case it was inv=-d for p>t !


A little quick and dirty test, just get the small q prime divisors for mp where q<10000*p
, and test mp/d, where d is the product of all found prime divisors of mp.

[CODE]
prp(n)=return(Mod(3,n)^n==Mod(3,n))

prpc(maxp)={t=512;num_large_d=0;forprime(p=t,maxp,
d=1;mp=2^p-1;forstep(q=2*p+1,10000*p,2*p,if(isprime(q)&&lift(Mod(2,q)^p)==1,d*=q));
if(d>1&&d<mp,res=lift((Mod(3,mp)^mp));s=lift(Mod(3,mp)^d)%(mp/d);
w=(d*(s-res))%(2^t);
if(d>2^t,num_large_d+=1);
if(prp(mp/d),print(["Found prp ",p,d]));
if(w<d&&2^t>d,print(["Candidate! ",p,d]))));
print("Number of large divisors=",num_large_d);}
prpc(10000)
[/CODE]

The LLT variant is similar to the above test,
only s has a slightly complicated formula.

[CODE]
prp(n)=return(Mod(3,n)^n==Mod(3,n))

fmult(u,v,p)=return([(u[1]*v[1]+3*u[2]*v[2])%(2^p-1),(u[1]*v[2]+u[2]*v[1])%(2^p-1)])

fpowmod(v,e,p)={
ret=[1,0];h=v;
while(e,
if(e%2,ret=fmult(ret,h,p));
h=fmult(h,h,p);e>>=1);
return(ret)}

LLT(n,x0,it)={mp=2^n-1;x=Mod(x0,mp);for(i=1,it,x=x^2-2);return(lift(x))}

llc(maxp)={t=512;forprime(p=t,maxp,
d=1;mp=2^p-1;forstep(q=2*p+1,10000*p,2*p,if(isprime(q)&&lift(Mod(2,q)^p)==1,d*=q));
if(d>1&&d<mp,
res=LLT(p,4,p);
s=(2*(fpowmod([2,1],d+kronecker(3,mp/d),p)[1]))%(mp/d);
w=(d*(s-res))%(2^t);
if(d>2^t,num_large_d+=1);
if(prp(mp/d),print(["Found prp ",p,d]));
if(w<d&&2^t>d,print(["Candidate! ",p,d]))));
print("Number of large divisors=",num_large_d);}
llc(10000)
[/CODE]

the output in both case:
[CODE]
["Found prp ", 881, 26431]
["Candidate! ", 881, 26431]
["Found prp ", 1459, 1362463807]
["Candidate! ", 1459, 1362463807]
["Found prp ", 3041, 135391639199]
["Candidate! ", 3041, 135391639199]
["Found prp ", 3359, 6719]
["Candidate! ", 3359, 6719]
["Found prp ", 4127, 154517628583]
["Candidate! ", 4127, 154517628583]
["Found prp ", 4243, 101833]
["Candidate! ", 4243, 101833]
["Found prp ", 7673, 2563193011919]
["Candidate! ", 7673, 2563193011919]
Number of large divisors=0
[/CODE]
so found all prp, and there was no false hit.

ps.
say if 2^t>d>2^(t-1) then you'll have a good chance that w<d, but this doesn't mean that you have found a prime,
because with more than 50% chance w<d will be true.
Note that for almost all known smallish p value (p<10^9 or so) it is true that d<2^512.
So you'd only need additional PRP-C tests for those d>2^512 and for those that survived this new test.
For this we'd need to store RES512, if we'd choose t=512.

And not forget that with this mod 2^t shrink we are (still) weaker than a standard Fermat test, but not much weaker as
you can see from the example runs.

Great! And the questions begin
 

Excellent. While this does not seem to impact the primality-unknown PRP or LL, as I incompletely understand it, it should really help the PRP-C and LL-C progress. Like your previous insights, it will take a while to understand, implement, roll out, and yield advantages. And that long process starts when you introduce such things, so keep them coming!

So now some questions,

A) step one is the PRP test if no factors are found yet. Always save res512 from that?

B) does the first cofactor run generate the res512 and send it to the (enhanced to accept it) PrimeNet server?

C) what's the language of your code sections for PRP-c and ll-c?

D) For confidence in accurate transmission of res512 values, do we use some agreed standard like CRC32 from the start?

E) Is addition of res512 storage practical on the server? A question for Madpoo et al

F) What's inv there? Inverse element in some group?

"Not get confused" (Too late...;)

[QUOTE=kriesel;490360]Excellent. While this does not seem to impact the primality-unknown PRP or LL, as I incompletely understand it, it should really help the PRP-C and LL-C progress. Like your previous insights, it will take a while to understand, implement, roll out, and yield advantages. And that long process starts when you introduce such things, so keep them coming![/QUOTE]

There is no very deep Maths here, we're working in different rings, that could be the reason why this could be hard/new.

[QUOTE=kriesel;490360]
A) step one is the PRP test if no factors are found yet. Always save res512 from that?
[/QUOTE]
It is the same test even if a factor is known, we need 3^mp mod mp.
[QUOTE=kriesel;490360]
B) does the first cofactor run generate the res512 and send it to the (enhanced to accept it) PrimeNet server? [/QUOTE]

There would be only at most one prp/ll test (not counting the double checks), independently from the
number of factors of mp. And yes, we'd need to store it somewhere.

[QUOTE=kriesel;490360]
C) what's the language of your code sections for PRP-c and ll-c?[/QUOTE]
Used Pari-Gp, it's available both on Linux and Windows.

[QUOTE=kriesel;490360]
D) For confidence in accurate transmission of res512 values, do we use some agreed standard like CRC32 from the start?[/QUOTE]
Yes, you can do that.

[QUOTE=kriesel;490360]
F) What's inv there? Inverse element in some group?[/QUOTE]
It's more than a group, it is a ring, but yes ofcourse inv is the multiplicative inverse of N/d in Z[2^t].

[QUOTE=kriesel;490360]
"Not get confused" (Too late...;)[/QUOTE]

In the general case there is (res-s) factor in the "final" formula, and in the Mersenne case there was
s-res in [4], a negative sign came from inv, because inv=-d.


Returning to LL case:
since we need p iterations here and we can't get even the low 64 bits of S[p] if you know (S[p-2] mod 2^64),
we really need that two (cheap) LL iterations, and we need to keep those S[p-2] residues also, because:
1. for double checking the old tests
2. S[p]=2 is necessity if mp is prime, but I don't see the sufficiency [it could be true, but this needs a proof].

[QUOTE=R. Gerbicz;490362]It is the same test even if a factor is known, we need 3^mp mod mp.
[/QUOTE]

Just correction, we need to store res'=(3^mp mod mp) mod 2^t, where we use the smallest non-negative integer at (x mod n). And in [4] we actually used this res' value:
[CODE]
[4]: w=(d*(s-res') mod 2^t)
[/CODE]
Just seeing the PARI-Gp code they have also used this full res residue not the shrinked version, res'.


Not written, but a general trick works here if you used 3^(mp+1) or 3^(mp-1) instead of 3^mp, then there are 3 possibilities for res'.

Just for completeness, the two corrected tests, when we store the shrinked residues:
[CODE]
prp(n)=return(Mod(3,n)^n==Mod(3,n))

prpc(maxp)={t=512;num_large_d=0;forprime(p=t,maxp,
d=1;mp=2^p-1;forstep(q=2*p+1,10000*p,2*p,if(isprime(q)&&lift(Mod(2,q)^p)==1,d*=q));
if(d>1&&d<mp,res=(lift((Mod(3,mp)^mp)))%(2^t);s=lift(Mod(3,mp)^d)%(mp/d);
w=(d*(s-res))%(2^t);
if(d>2^t,num_large_d+=1);
if(prp(mp/d),print(["Found prp ",p,d]));
if(w<d&&2^t>d,print(["Candidate! ",p,d]))));
print("Number of large divisors=",num_large_d);}
prpc(10000)
[/CODE]

[CODE]
prp(n)=return(Mod(3,n)^n==Mod(3,n))

fmult(u,v,p)=return([(u[1]*v[1]+3*u[2]*v[2])%(2^p-1),(u[1]*v[2]+u[2]*v[1])%(2^p-1)])

fpowmod(v,e,p)={
ret=[1,0];h=v;
while(e,
if(e%2,ret=fmult(ret,h,p));
h=fmult(h,h,p);e>>=1);
return(ret)}

LLT(n,x0,it)={mp=2^n-1;x=Mod(x0,mp);for(i=1,it,x=x^2-2);return(lift(x))}

llc(maxp)={t=512;forprime(p=t,maxp,
d=1;mp=2^p-1;forstep(q=2*p+1,10000*p,2*p,if(isprime(q)&&lift(Mod(2,q)^p)==1,d*=q));
if(d>1&&d<mp,
res=LLT(p,4,p)%(2^t);
s=(2*(fpowmod([2,1],d+kronecker(3,mp/d),p)[1]))%(mp/d);
w=(d*(s-res))%(2^t);
if(d>2^t,num_large_d+=1);
if(prp(mp/d),print(["Found prp ",p,d]));
if(w<d&&2^t>d,print(["Candidate! ",p,d]))));
print("Number of large divisors=",num_large_d);}
llc(10000)
[/CODE]

[QUOTE=R. Gerbicz;490362]
2. S[p]=2 is necessity if mp is prime, but I don't see the sufficiency [it could be true, but this needs a proof].[/QUOTE]

[TEX]2\equiv 2^2-2\bmod{M_n}\equiv(-2)^2-2\bmod{M_n}[/TEX] it needs to fall back to the -2 case eventually because if not you won't get anything but a list of 4's. The only other possible way around this is the previous mersenne's square (congruent to the power of 2 next door mod the higher mersenne) or another pair of a and -a having their squares equal to 4 mod [TEX]M_n[/TEX] once you fall to the -2 path backwards you hit 0 mod [TEX]M_n[/TEX] as a square and sqrt.

[QUOTE=science_man_88;490366][TEX]2\equiv 2^2-2\bmod{M_n}\equiv(-2)^2-2\bmod{M_n}[/TEX] it needs to fall back to the -2 case eventually because if not you won't get anything but a list of 4's. The only other possible way around this is the previous mersenne's square (congruent to the power of 2 next door mod the higher mersenne) or another pair of a and -a having their squares equal to 4 mod [TEX]M_n[/TEX] once you fall to the -2 path backwards you hit 0 mod [TEX]M_n[/TEX] as a square and sqrt.[/QUOTE]

You proved (almost) nothing, surf the web less, and learn more Math, much more.
Ofcourse even with a proof as stated above for obvious reason(s) we'd need to store the RES64 for S[p-2]. For my algorithm this question isn't that very interesting (at database or in other level), just noted as a somewhat nice open(?) question.

Would it be useful/desirable for LL and PRP software to report to the server a larger 512-bit residue instead of the current 64-bit?

[QUOTE=preda;490668]Would it be useful/desirable for LL and PRP software to report to the server a larger 512-bit residue instead of the current 64-bit?[/QUOTE]

Yes, but the thing needed is not just a bigger version of RES64. For PRP, you need 3^Mp (mod Mp), not 3^(Mp-1) (mod Mp). For LL, you need two additional iteration of s^2-2, before you extract the RES 512.

But if you're going this route, might as well futureproof it and go with a bigger RES1024 or something ( this will fit in a varbinary(128) in sql server, so storage cost isn't high ).

[QUOTE=axn;490670]Yes, but the thing needed is not just a bigger version of RES64. For PRP, you need 3^Mp (mod Mp), not 3^(Mp-1) (mod Mp). For LL, you need two additional iteration of s^2-2, before you extract the RES 512.

But if you're going this route, might as well futureproof it and go with a bigger RES1024 or something ( this will fit in a varbinary(128) in sql server, so storage cost isn't high ).[/QUOTE]

Anyway right now I'm computing 3^(Mp + 1), and at the end do a division-by-9 to derive the 3^(Mp - 1).

Also, I can do either 512 or 1024. I'd like to know what George thinks. If he gives the green light I can do that.

(also we'd need to chose a key for reporting it in the JSON. Probably in addition to the current res64, in a first phase).

It would be nice if we had a way to derive the current res64 from the new res512 or res1024. My feeling is that it could be done with some guard high-bits ("before 0") that would sink the carry propagation from the div-3 or mul-3.

Or, could the new "big-residue" be not starting at bit-0, but e.g. centered on bit-0? i.e. shifted some amount such that it contains bits on both sides of the current res64.

[QUOTE=axn;490670]
But if you're going this route, might as well futureproof it and go with a bigger RES1024 or something ( this will fit in a varbinary(128) in sql server, so storage cost isn't high ).[/QUOTE]
Maybe there could be a dozen or so such p, where mp is not completely factored, but d>2^512.
We could switch to res1024 only for say p>2^32.
And ofcourse you'd to store it only if you have a completed LL/PRP test to save space in database.


[QUOTE=preda;490676]It would be nice if we had a way to derive the current res64 from the new res512 or res1024. My feeling is that it could be done with some guard high-bits ("before 0") that would sink the carry propagation from the div-3 or mul-3.[/QUOTE]

See,
[CODE]let f(p,e)=mp=2^p-1;return(lift(Mod(3,mp)^e)%(2^512))[/CODE]

where I need this for e=mp
then
f(p,e+1)=(3*f(p,e)-c*mp)%(2^512), but mp==-1 mod 2^512, for p>512.
so we can write:
[CODE]
f(p,e+1)=(3*f(p,e)+c)%(2^512) for some c=0,1,2 and if you solve it for f(p,e) then:
f(p,e)=lift(Mod((f(p,e+1)-c),2^512)/3) for some c=0,1,2
[/CODE]
So you can't get the exact residue, but you have "only" 3 choices for the residue,
(after the additional mod 2^64), ofcourse you would have 3 residues even for mod 2^512.

The situation is similar if you have general number (not Mersenne), or base!=3 or p<512.
(but that is really boring, since those are factored), or you have f(p,e+2) and want f(p,e);
but in this case you can also just iterate the above method, you get 9 possible residues for f(p,e).
And again you can say in the case of LL test nothing from RES512 at iteration p, due to the nonlinear squaring. [That is why we'd need to store also the "standard" RES64 at iteration (p-2)].