mersenneforum.org  

Go Back   mersenneforum.org > Prime Search Projects > Wagstaff PRP Search

Closed Thread
 
Thread Tools
Old 2018-08-02, 15:28   #67
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

101111011112 Posts
Default

Quote:
Originally Posted by GP2 View Post
Ideally, the PRP residues should be recorded with a lot more than 64 bits. If the only goal is double-check verifiability, then 64 bits is plenty; however, larger residues might let us do quick PRP-cofactor checking for newly discovered factors. (I presume the analysis at the linked page for Mersennes applies similarly for Wagstaffs, or is "inv" different?)

So, at least 512 bits, or with an eye to future decades, who knows, maybe even 2048. Storage is very cheap.
As noted in that post this works for all odd N (like the Suyama test), for Wagstaff there is also a shorter formula w:
you should use w=(3*d*(res-s))%(2^t) instead of w=(d*(s-res))%(2^t); (what used for Mersenne numbers),
if p>t=512 (or whatever you use).

If we wouldn't want a fixed bitlength residue, then RES_t is a good choice where say t=16*floor(log2(p)), and then you have to redo the prp test only when d>2^t or when my test is showing that N/d could be a (probable) prime.
R. Gerbicz is offline  
Old 2018-08-02, 19:40   #68
GP2
 
GP2's Avatar
 
Sep 2003

5×11×47 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
As noted in that post this works for all odd N (like the Suyama test), for Wagstaff there is also a shorter formula w:
you should use w=(3*d*(res-s))%(2^t) instead of w=(d*(s-res))%(2^t); (what used for Mersenne numbers),
if p>t=512 (or whatever you use).
I have one more question:

You define
res = 3^mp mod mp
and in the formula for w we can substitute:
res_t = res mod 2^t

But for either Mersenne PRP=1,2,n,-1 or Wagstaff PRP=1,2,n,1,"3" , mprime actually calculates:
mp = 2^p ± 1
resx = 3^(mp-1) mod mp
and prints out:
resx_t = resx mod 2^t (t = 64 currently)

For instance for 2^523 − 1 it gives hex D749C83A364C1462 and for 2^523 + 1 it gives hex 388104444C578BC6


So
(3 * resx) mod mp == res
but
(3 * resx_t) mod 2^t != res_t

When t=512 the difference observed is small, only 0 or 1 or 2, but I think it would affect the calculation.

So how can we derive res_t from the available printed value resx_t ?

Last fiddled with by GP2 on 2018-08-02 at 19:42
GP2 is offline  
Old 2018-08-02, 22:30   #69
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

72×31 Posts
Default

Quote:
Originally Posted by GP2 View Post
So how can we derive res_t from the available printed value resx_t ?
Right, and there was another problem with that approach, because the exponent=(2^p+1)/3 is not a power of two
, and that would block my error check. Fortunately there is a workaround:
let N=(2^p+1)/3, if N/d is prime then with Fermat's little theorem:
3^(N/d)==3 mod (N/d), from this
3^(2^p)==3^(3*d-1) mod (N/d) and the rest is similar to the Mersenne case:
(the special case, when b=2: N=(k*2^n+c)/D is similar)

Code:
Let N=(2^p+1)/3
res=(3^(2^p) mod (2^p+1)) mod (2^t)
s=3^(3*d-1) mod (N/d)
w=3*d*(res-s) mod (2^t)
if 3*d<2^t and w>=3*d then N/d is composite,
   otherwise if 3*d>2^t then we know nothing (too large divisor!) or
   3*d<2^t and w<3*d, then N/d can be prime (we are weaker than any Fermat test, because of the "shrinked" residue)
and a sample test for Wagstaff:

Code:
prp(n)=return(Mod(3,n)^n==Mod(3,n))

prpc(maxp)={t=512;num_large_d=0;c0=c1=0;forprime(p=t,maxp,
d=1;N=(2^p+1)/3;forstep(q=2*p+1,10000*p,2*p,if(isprime(q)&&lift(Mod(2,q)^p+1)==0,d*=q));
if(d>1&&d<N,res=lift((Mod(3,2^p+1)^(2^p)))%(2^t);s=lift(Mod(3,N/d)^(3*d-1));
w=(3*d*(res-s))%(2^t);
if(3*d>2^t,num_large_d+=1);
if(3*d<2^t&&w>=3*d,c0++);
if(3*d<2^t&&w<3*d,c1++);
if(prp(N/d),print(["Found prp ",p,d]));
if(w<3*d&&2^t>3*d,print(["Candidate! ",p,d]))));
print("Number of large divisors=",num_large_d);
print(c0);print(c1)}
prpc(10000)
R. Gerbicz is offline  
Old 2018-08-03, 02:58   #70
axn
 
axn's Avatar
 
Jun 2003

520310 Posts
Default

Quote:
Originally Posted by diep View Post
Possible. Until 10M you won't find anything new.
There is a possibility, however remote, of a hidden FFT bug that could miss a prime. For mersennes, we have doublecheck with different rotation, that can detect the issue even if there is FFT bug. But I don't believe Wagstaff numbers can do the bit rotation trick (corrections welcome!). So your confidence in the integrity of the result is predicated upon your confidence about the software.

Having said that, a single pass with GEC PRP should be sufficient to rule out such issues.
axn is offline  
Old 2018-08-03, 03:45   #71
GP2
 
GP2's Avatar
 
Sep 2003

50318 Posts
Default

I wrote some Python code for the Gerbicz quick-test method which determines in many cases that a PRP cofactor is composite without actually doing a PRP cofactor test.

I find Python a little more comfortable to use than PARI/GP.

The code is attached (gerbicz_quickprp.tar.gz)
Attached Files
File Type: gz gerbicz_quickprp.tar.gz (70.1 KB, 111 views)

Last fiddled with by GP2 on 2018-08-03 at 04:03
GP2 is offline  
Old 2018-08-03, 12:01   #72
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

72×31 Posts
Default

Quote:
Originally Posted by axn View Post
But I don't believe Wagstaff numbers can do the bit rotation trick (corrections welcome!).
As I know only the newer versions of p95 includes rotation for PRP test, the older version hasn't got this, or we can say that they used the constant shift=0.

Quote:
Originally Posted by GP2 View Post
I wrote some Python code for the Gerbicz quick-test method which determines in many cases that a PRP cofactor is composite without actually doing a PRP cofactor test.
Looks like good, have you checked that this finds all "known" (probable) primes for Mersenne/Wagstaff case in your range? (ofcourse not including the numbers where Mp,W(p) is prime). My super light sieve only finds prp numbers at p=557,1171,1627,3329 for Wagstaff case in that p=512..5000 range.
R. Gerbicz is offline  
Old 2018-08-03, 14:11   #73
GP2
 
GP2's Avatar
 
Sep 2003

5×11×47 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
As I know only the newer versions of p95 includes rotation for PRP test, the older version hasn't got this, or we can say that they used the constant shift=0.
For PRP, the JSON output has a shift-count:

Code:
{"status":"C", "exponent":78358739, "worktype":"PRP-3", "res64":"EF71CC587195371C", "residue-type":1, "fft-length":4587520, "shift-count":75805076, "error-code":"00000000", ...
For PRP-CF, the JSON output also has a shift-count:

Code:
{"status":"C", "exponent":6444727, "known-factors":"373794167,176482404169", "worktype":"PRP-3", "res64":"B3948091270BC85C", "residue-type":5, "fft-length":327680, "shift-count":3602762, "error-code":"00000000", ...
But for Wagstaff, the JSON output doesn't display a shift-count:

Code:
{"status":"C", "k":1, "b":2, "n":31873, "c":1, "known-factors":"3", "worktype":"PRP-3", "res64":"130344E20CCBF5EE", "residue-type":5, "fft-length":1536, "error-code":"00000000", ...
Maybe it's just an omission? As long as b == 2 in k × bn + c, a non-zero shift count should be possible, no?


Quote:
Looks like good, have you checked that this finds all "known" (probable) primes for Mersenne/Wagstaff case in your range? (ofcourse not including the numbers where Mp,W(p) is prime). My super light sieve only finds prp numbers at p=557,1171,1627,3329 for Wagstaff case in that p=512..5000 range.
Yes, it flags all the PRPs in that range, and also flags a few false positives: for Mersenne 1489 and 3853, and for Wagstaff 1439. It's easy to see why, those exponents have a big list of factors, and the product of the factors is larger than 512 bits.

In the code I posted, the slowest step by far is the line where I calculate the 512-bit residue myself in Python, rather than getting it from a future version of mprime that will hopefully have the capability of displaying large-bitsize residues.

Code:
 res = pow(a, mp - 1, mp) % pow2_t
However, if we use the gmpy2 library for Python (interfaces to GMP), and make one small change:

Code:
from gmpy2 import mpz
mpz_a = mpz(a)
...
...
 res = pow(mpz_a, mp - 1, mp) % pow2_t
This makes the modular exponentiation use GMP rather than Python's much slower internal C code for large numbers. So now I can calculate 512-bit residues for exponents up to 100k in maybe a day or so, and then recheck that it doesn't miss any PRPs in this larger range.

Last fiddled with by GP2 on 2018-08-03 at 14:14
GP2 is offline  
Old 2018-08-03, 17:48   #74
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

72×31 Posts
Default

Quote:
Originally Posted by GP2 View Post
But for Wagstaff, the JSON output doesn't display a shift-count:

Code:
{"status":"C", "k":1, "b":2, "n":31873, "c":1, "known-factors":"3", "worktype":"PRP-3", "res64":"130344E20CCBF5EE", "residue-type":5, "fft-length":1536, "error-code":"00000000", ...
Maybe it's just an omission? As long as b == 2 in k × bn + c, a non-zero shift count should be possible, no?
Ouch, no, or I don't see the b=2 case, but for Wagstaff it is still working:
if you start with base=3*2^shift, where 0<=shift<2*p (but it is right to restrict to 0<=shift<p).
then the final residue is:
res=(3*2^shift)^(2^p)==3^(2^p)*2^(2*shift) mod (2^p+1)
so to get the original final (non-shifted) residue:
res'=res*2^(2*p-2*shift) mod (2^p+1).
R. Gerbicz is offline  
Old 2018-08-04, 00:51   #75
GP2
 
GP2's Avatar
 
Sep 2003

5×11×47 Posts
Default

I tested up to p = 50k for both Mersenne and Wagstaff.

With 2048-bit residues there are no false positives, and obviously no false negatives.

Updated programs and datafiles with large-bitsize residues are attached.
Attached Files
File Type: gz gerbicz_quickprp_ver2.tar.gz (3.16 MB, 114 views)
GP2 is offline  
Old 2018-08-04, 02:16   #76
axn
 
axn's Avatar
 
Jun 2003

112·43 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
Ouch, no, or I don't see the b=2 case, but for Wagstaff it is still working:
if you start with base=3*2^shift, where 0<=shift<2*p (but it is right to restrict to 0<=shift<p).
then the final residue is:
res=(3*2^shift)^(2^p)==3^(2^p)*2^(2*shift) mod (2^p+1)
so to get the original final (non-shifted) residue:
res'=res*2^(2*p-2*shift) mod (2^p+1).
This looks interesting, but you need to do additional work to recover the correct residue. In mersennes, you don't need to do any work at all. Having said that, it is just shifts, so not really a big deal. I'm pretty sure this is being not done currently, so more work for George (hint!)

HOWEVER, what is the math for recovering intermediate residues? Because there is an interest in getting intermediate residues at specific milestone as well.

Anyway, if we're using GEC, then the shift trick is superfluous. It doesn't bring any extra safety to the project. It is merely a nice-to-have.

Last fiddled with by axn on 2018-08-04 at 02:16 Reason: correctly -> currently
axn is offline  
Old 2018-08-04, 10:32   #77
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

72×31 Posts
Default

Quote:
Originally Posted by axn View Post
This looks interesting, but you need to do additional work to recover the correct residue.
That is the big deal, you are worrying about O(p) cost when you already done O(p^(2+eps)) work.

Quote:
Originally Posted by axn View Post
HOWEVER, what is the math for recovering intermediate residues? Because there is an interest in getting intermediate residues at specific milestone as well.
Similar maths, but why would you store the unshifted residue, when you continue with the shifted? I'd store the shifted residue and the shift number, that is more logical.

Quote:
Originally Posted by axn View Post
Anyway, if we're using GEC, then the shift trick is superfluous. It doesn't bring any extra safety to the project. It is merely a nice-to-have.
Correct.
R. Gerbicz is offline  
Closed Thread

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Testing Mersenne Primes with Elliptic Curves a nicol Math 3 2017-11-15 20:23
New Wagstaff PRP exponents ryanp Wagstaff PRP Search 26 2013-10-18 01:33
Hot tuna! -- a p75 and a p79 by Sam Wagstaff! Batalov GMP-ECM 9 2012-08-24 10:26
Statistical odds for prime in Wagstaff vs Mersenne diep Math 27 2010-01-13 20:18
Speed of P-1 testing vs. Trial Factoring testing eepiccolo Math 6 2006-03-28 20:53

All times are UTC. The time now is 18:24.


Tue Dec 7 18:24:50 UTC 2021 up 137 days, 12:53, 2 users, load averages: 1.49, 1.62, 1.52

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, 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.