 mersenneforum.org Status of Wagstaff testing? and testing Mersenne primes for Wagstaff-ness
 Register FAQ Search Today's Posts Mark Forums Read  2018-08-02, 15:28   #67
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

1,531 Posts Quote:
 Originally Posted by GP2 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.  2018-08-02, 19:40   #68
GP2

Sep 2003

5×11×47 Posts Quote:
 Originally Posted by R. Gerbicz 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  2018-08-02, 22:30   #69
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

1,531 Posts Quote:
 Originally Posted by GP2 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)  2018-08-03, 02:58   #70
axn

Jun 2003

122338 Posts Quote:
 Originally Posted by diep 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.  2018-08-03, 03:45   #71
GP2

Sep 2003

1010000110012 Posts 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 gerbicz_quickprp.tar.gz (70.1 KB, 119 views)

Last fiddled with by GP2 on 2018-08-03 at 04:03  2018-08-03, 12:01   #72
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

27738 Posts Quote:
 Originally Posted by axn 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 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.  2018-08-03, 14:11   #73
GP2

Sep 2003

5·11·47 Posts Quote:
 Originally Posted by R. Gerbicz 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  2018-08-03, 17:48   #74
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

1,531 Posts Quote:
 Originally Posted by GP2 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).  2018-08-04, 00:51   #75
GP2

Sep 2003

5×11×47 Posts 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 gerbicz_quickprp_ver2.tar.gz (3.16 MB, 122 views)  2018-08-04, 02:16   #76
axn

Jun 2003

10100100110112 Posts Quote:
 Originally Posted by R. Gerbicz 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
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  2018-08-04, 10:32   #77
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

1,531 Posts Quote:
 Originally Posted by axn 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 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 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.  Thread Tools Show Printable Version Email this Page Similar Threads Thread Thread Starter Forum Replies Last Post a nicol Math 3 2017-11-15 20:23 ryanp Wagstaff PRP Search 26 2013-10-18 01:33 Batalov GMP-ECM 9 2012-08-24 10:26 diep Math 27 2010-01-13 20:18 eepiccolo Math 6 2006-03-28 20:53

All times are UTC. The time now is 06:01.

Fri Jan 28 06:01:30 UTC 2022 up 189 days, 30 mins, 2 users, load averages: 1.42, 1.46, 1.49