mersenneforum.org  

Go Back   mersenneforum.org > Great Internet Mersenne Prime Search > Math

Reply
 
Thread Tools
Old 2018-06-22, 20:59   #1
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

1,409 Posts
Default 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 http://mersenneforum.org/showthread....732#post471732
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:
https://www.mersenne.org/report_expo...1344943&full=1
https://www.mersenne.org/report_expo...2389403&full=1
https://www.mersenne.org/report_expo...2462281&full=1
https://www.mersenne.org/report_expo...3009703&full=1
https://www.mersenne.org/report_expo...4373911&full=1

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.
[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).
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
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)
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)
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)
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
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.
R. Gerbicz is offline   Reply With Quote
Old 2018-06-23, 16:30   #2
kriesel
 
kriesel's Avatar
 
"TF79LL86GIMPS96gpu17"
Mar 2017
US midwest

22×32×127 Posts
Default 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...;)

Last fiddled with by kriesel on 2018-06-23 at 16:33
kriesel is offline   Reply With Quote
Old 2018-06-23, 19:22   #3
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

1,409 Posts
Default

Quote:
Originally Posted by kriesel View Post
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!
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:
Originally Posted by kriesel View Post
A) step one is the PRP test if no factors are found yet. Always save res512 from that?
It is the same test even if a factor is known, we need 3^mp mod mp.
Quote:
Originally Posted by kriesel View Post
B) does the first cofactor run generate the res512 and send it to the (enhanced to accept it) PrimeNet server?
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:
Originally Posted by kriesel View Post
C) what's the language of your code sections for PRP-c and ll-c?
Used Pari-Gp, it's available both on Linux and Windows.

Quote:
Originally Posted by kriesel View Post
D) For confidence in accurate transmission of res512 values, do we use some agreed standard like CRC32 from the start?
Yes, you can do that.

Quote:
Originally Posted by kriesel View Post
F) What's inv there? Inverse element in some group?
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:
Originally Posted by kriesel View Post
"Not get confused" (Too late...;)
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].
R. Gerbicz is offline   Reply With Quote
Old 2018-06-23, 19:31   #4
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

1,409 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
It is the same test even if a factor is known, we need 3^mp mod mp.
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)
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:
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)

Last fiddled with by R. Gerbicz on 2018-06-23 at 19:38 Reason: adds+grammar
R. Gerbicz is offline   Reply With Quote
Old 2018-06-23, 19:48   #5
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

202618 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
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].
2\equiv 2^2-2\bmod{M_n}\equiv(-2)^2-2\bmod{M_n} 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 M_n once you fall to the -2 path backwards you hit 0 mod M_n as a square and sqrt.

Last fiddled with by science_man_88 on 2018-06-23 at 19:49
science_man_88 is offline   Reply With Quote
Old 2018-06-24, 14:43   #6
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

26018 Posts
Default

Quote:
Originally Posted by science_man_88 View Post
2\equiv 2^2-2\bmod{M_n}\equiv(-2)^2-2\bmod{M_n} 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 M_n once you fall to the -2 path backwards you hit 0 mod M_n as a square and sqrt.
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.
R. Gerbicz is offline   Reply With Quote
Old 2018-06-27, 05:43   #7
preda
 
preda's Avatar
 
"Mihai Preda"
Apr 2015

22·32·5·7 Posts
Default

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?
preda is online now   Reply With Quote
Old 2018-06-27, 08:29   #8
axn
 
axn's Avatar
 
Jun 2003

22·11·107 Posts
Default

Quote:
Originally Posted by preda View Post
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?
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 ).
axn is offline   Reply With Quote
Old 2018-06-27, 09:22   #9
preda
 
preda's Avatar
 
"Mihai Preda"
Apr 2015

22×32×5×7 Posts
Default

Quote:
Originally Posted by axn View Post
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 ).
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).

Last fiddled with by preda on 2018-06-27 at 09:22
preda is online now   Reply With Quote
Old 2018-06-27, 09:28   #10
preda
 
preda's Avatar
 
"Mihai Preda"
Apr 2015

22·32·5·7 Posts
Default

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.

Last fiddled with by preda on 2018-06-27 at 09:35
preda is online now   Reply With Quote
Old 2018-06-27, 10:26   #11
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

101100000012 Posts
Default

Quote:
Originally Posted by axn View Post
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 ).
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:
Originally Posted by preda View Post
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.
See,
Code:
let f(p,e)=mp=2^p-1;return(lift(Mod(3,mp)^e)%(2^512))
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
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)].

Last fiddled with by R. Gerbicz on 2018-06-27 at 10:36
R. Gerbicz is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Daylight Saving Time davieddy Lounge 27 2015-02-26 18:45
Prime95 NOT Saving!!! Primeinator Information & Answers 9 2008-09-22 19:00
Not Saving BioRules Information & Answers 9 2008-05-31 13:52
Saving computation in ECM dave_dm Factoring 8 2004-06-12 14:18
saving over a network crash893 Software 11 2004-05-06 14:15

All times are UTC. The time now is 21:00.

Tue Oct 20 21:00:18 UTC 2020 up 40 days, 18:11, 1 user, load averages: 1.94, 1.98, 1.91

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.