gpuOwL: an OpenCL program for Mersenne primality testing

[Prime95]

So gpuOwl is printing non-standard residues for a Fermat test. That's not particularly good.

Here's my proposal. A Fermat probable prime test for N is 3^(N-1) = 1 mod N. Iteration 0 is 3. Using a standard left-to-right exponentiation algorithm, we get iteration 1 = 3^(top 2 bits of N), iteration 2 = 3^(top 3 bits of N), and so forth.

I believe this is what LLR, PFGW, and prime95 all do for interim and final residues.

You can still do the no-mul-by-3's implementation. That is, N = 2^P-1, and you will calculate 3^(N+1) mod N. To print a standard iteration 1 residue, output 3^(top 3 bits of N+1)/3, etc. That is, you'll do divide-by-3 to output interim residues and divide-by-9 to get the final residue. This is what I'm doing in the next prime95 release -- testing it now.

[axn]

Originally Posted by Prime95

So gpuOwl is printing non-standard residues for a Fermat test. That's not particularly good.

Here's my proposal. A Fermat probable prime test for N is 3^(N-1) = 1 mod N. Iteration 0 is 3. Using a standard left-to-right exponentiation algorithm, we get iteration 1 = 3^(top 2 bits of N), iteration 2 = 3^(top 3 bits of N), and so forth.

I believe this is what LLR, PFGW, and prime95 all do for interim and final residues.

You can still do the no-mul-by-3's implementation. That is, N = 2^P-1, and you will calculate 3^(N+1) mod N. To print a standard iteration 1 residue, output 3^(top 3 bits of N+1)/3, etc. That is, you'll do divide-by-3 to output interim residues and divide-by-9 to get the final residue. This is what I'm doing in the next prime95 release -- testing it now.

What all types of test are you planning to cover under the "Gerbicz error check" in P95? Only LL or the general k*b^n+/-1 form? Because, your approach will not work for bases other than 2.

[preda]

Originally Posted by Prime95

So gpuOwl is printing non-standard residues for a Fermat test. That's not particularly good.

Here's my proposal. A Fermat probable prime test for N is 3^(N-1) = 1 mod N. Iteration 0 is 3. Using a standard left-to-right exponentiation algorithm, we get iteration 1 = 3^(top 2 bits of N), iteration 2 = 3^(top 3 bits of N), and so forth.

I believe this is what LLR, PFGW, and prime95 all do for interim and final residues.

You can still do the no-mul-by-3's implementation. That is, N = 2^P-1, and you will calculate 3^(N+1) mod N. To print a standard iteration 1 residue, output 3^(top 3 bits of N+1)/3, etc. That is, you'll do divide-by-3 to output interim residues and divide-by-9 to get the final residue. This is what I'm doing in the next prime95 release -- testing it now.

Apparently what I compute is a Strong PRP https://primes.utm.edu/prove/prove2_3.html
which is stronger than a Fermat PRP. Given no additional cost and some benefit, then why not do the stronger test?

I do not calculate 3^(N+1) mod N, but 3^((N+1)/2) mod N. I think that's exactly the difference that makes it "strong", by halving the probability of false positive.

[Prime95]

Originally Posted by axn

What all types of test are you planning to cover under the "Gerbicz error check" in P95? Only LL or the general k*b^n+/-1 form? Because, your approach will not work for bases other than 2.

I'll support all (k*b^n+c)/known_factors.

You are correct, the no-mul-by-3's code is only applies to Mersenne numbers with no known factors. In all other cases, I'll do a standard Fermat PRP test with Gerbicz recovery.

[axn]

Originally Posted by Prime95

I'll support all (k*b^n+c)/known_factors.

You are correct, the no-mul-by-3's code is only applies to Mersenne numbers with no known factors. In all other cases, I'll do a standard Fermat PRP test with Gerbicz recovery.

Gerbicz method can only be applied to a sequence of powerings (a^(b^N)), so standard PRP is not compatible with it. Which is why gpuOwl is implementing 3^(N+1) so that the computation sequence is an uninterrupted sequence of squarings only. Now, as luck would have it, we can flip back and forth between standard PRP and powerings-only mode with a simple division/multiplication of 3 (in the case of base 2 only). But that will not work in case of other bases, nor will it work in the case where there is a factor.

[preda]

Originally Posted by Prime95

So gpuOwl is printing non-standard residues for a Fermat test. That's not particularly good.

Here's my proposal. A Fermat probable prime test for N is 3^(N-1) = 1 mod N. Iteration 0 is 3. Using a standard left-to-right exponentiation algorithm, we get iteration 1 = 3^(top 2 bits of N), iteration 2 = 3^(top 3 bits of N), and so forth.

I believe this is what LLR, PFGW, and prime95 all do for interim and final residues.

You can still do the no-mul-by-3's implementation. That is, N = 2^P-1, and you will calculate 3^(N+1) mod N. To print a standard iteration 1 residue, output 3^(top 3 bits of N+1)/3, etc. That is, you'll do divide-by-3 to output interim residues and divide-by-9 to get the final residue. This is what I'm doing in the next prime95 release -- testing it now.

When doing the mul-by-3 implementation, then for Mersenne numbers you should probably test 3^((N-1)/2) == -1, which is stronger than 3^(N-1) == 1.

[R. Gerbicz]

Originally Posted by axn

Because, your approach will not work for bases other than 2.

No and yes!
No, because my method works for any b, but yes because when we want to use this type
of ladder check we need to use b-th powermod calculation, this gives
that I will use (at least!) n*ceil(log2(b)) squaremod/mulmod, while
at a standard PRP test you are using (much) fewer, only n*log2(b) squaremods.

the core of the algorithm in b!=2:

u[i]=a^(b^i) mod N
d[t]=prod(i=0,t,u[L*i]) mod N
d[t]=d[t-1]*u[L*t] mod N
d[t]=a*d[t-1]^(b^L) mod N  (to do the error check)

With this a^(k*b^n)==u[n]^k mod N. (we delayed the k-th powermod).

New idea: you can save some work here: use say e=1000 and set B=b^e, L=1, but use error
checking after each 500 iteration. With such a huge B we reduce the overhead of the
b-th powermod to B-th powermod with (optimal or close to that) additon chain (https://en.wikipedia.org/wiki/Additi...exponentiation)
And you need to do floor(n/e) iterations, after that do an additional (n mod e)
b-th powermod to get a^(b^n) mod N, and final k-th powermod.

But I'm not sure that even with this it is worth to do this error checking in the b!=2 test.

Ofcourse when b=2 we are already in an optimal case.
Multiple the final residue by a^(c-1) mod N if you want the Fermat residue if your number is N=k*b^n+c

Originally Posted by axn

But that will not work in case of other bases, nor will it work in the case where there is a factor.

No. When you know a nontrivial d>1 factor of N=k*b^n+c then to test the primality of N/d, just
calculate a^N mod N what you would do in the d=1 case.
If (N/d) is Fermat prp number using base=a, then
a^(N/d)==a mod (N/d) this gives that
a^N==a^d mod (N/d),
so in the test get a^N mod N, then reduce it mod (N/d).
if it is not a^d mod (N/d), then N/d is not a Fermat prp number, the other direction is not true, because
a^N==a^d mod (N/d) doesn't give that
a^(N/d)==a mod (N/d), but it isn't a much weaker test.
Even in the thousand digits of range you wouldn't see a false hit in your life.
If this type of test passed then you can still do a true Fermat test for (N/d).

The above is a known trick from the standard Suyama test, what we use for example for
Fermat numbers, in the test, if we know the residue we get a fast prp test for N/d but that is also
weaker than a (single) Fermat test.

Originally Posted by preda

When doing the mul-by-3 implementation, then for Mersenne numbers you should probably test 3^((N-1)/2) == -1, which is stronger than 3^(N-1) == 1.

For Mersenne prime(!) these are equivalent tests, for others I don't know even a single false hit, so at least in practice these are the same.

[axn]

Originally Posted by R. Gerbicz

No and yes!.

No, I don't mean that we can't implement your test. I mean that we can't do what George wants to do while also implementing this test. What George wants to do is, IIRC, to output residues (final & intermmediate) that corresponds to the computation 3^(N-1) for arbitrary N=(k*b^n+c)/f. Your test requires that we evaluate (3^k)^(b^n), and then equate it to 3^(f-c). These two sequences produce residues that are incompatible -- except for the case of f=1, b=2, where we can use George's trick to transform between the two.

[R. Gerbicz]

Originally Posted by axn

No, I don't mean that we can't implement your test. I mean that we can't do what George wants to do while also implementing this test. What George wants to do is, IIRC, to output residues (final & intermmediate) that corresponds to the computation 3^(N-1) for arbitrary N=(k*b^n+c)/f. Your test requires that we evaluate (3^k)^(b^n), and then equate it to 3^(f-c). These two sequences produce residues that are incompatible -- except for the case of f=1, b=2, where we can use George's trick to transform between the two.

Huge loss, but do you know that for f!=1 we have no (very) fast true primality test, like for example we still don't have a fast test for (10^n-1)/9 .

Btw my test could be even faster than George's test, if he really computes the Fermat residue for f>1, why? Because for b=2 I'm not doing mult by 2. *

Originally Posted by axn

Your test requires that we evaluate (3^k)^(b^n), and then equate it to 3^(f-c).

I'm prefering to write "a" instead of 3, because for Proth numbers we have a true primality test, and in those cases it is possible that 3 is a quadratic residue, so not usable.

* ps. the slowdown is much smaller if you do only a 64 bits mult once in 64 iterations. Maybe there could be other methods to reduce more this cost.

[Prime95]

Originally Posted by R. Gerbicz

the core of the algorithm in b!=2:

u[i]=a^(b^i) mod N

Yes, I overlooked that this step requires the costly full mulmod rather than the easy mul-by-small-const. So, I concur that for b!=2 prime95 will not support Gerbicz error checks.

[Prime95]

Here is my PRP residue proposal:

There are (at least) 5 PRP residue types:
1; Fermat PRP. N is PRP if a^(N-1) = 1 mod N
2: SPRP variant, N is PRP if a^((N-1)/2) = +/-1 mod N
3: Fermat PRP variant. N is PRP if a^(N+1) = a^2 mod N
4: SPRP variant. N is PRP if a^((N+1)/2) = +/-a mod N
5: Fermat PRP cofactor variant, N/d is PRP if a^(N-1) = a^d mod N/d

I encourage programs to return type 1 residues as that has been the standard for prime95, PFGW, LLR for many years.

I think Prime95 v29.4 will have options to support all of the above residue types for doing double-checks. PrimeNet will store residue type along with the residue -- need another JSON argument. When issuing PRP double-check assignments, the server will send the PRP base and desired residue type to the client. For first-time PRP assignments, the server will not specify a PRP base or residue type.

Interim residues are application dependent. However, it is strongly encouraged to report Fermat PRP residues using a standard binary powering algorithm. That is, iteration 0 is a, iteration 1 is a^(top 2 bits of N-1), etc.

The server will soon have about 50 type-4 residues from gpuOwl, which prime95 should be able to verify someday. While I'd prefer gpuOwl converts to type-1 residues, it is not a requirement.

The server also has several thousand type-1 PRP cofactor residues. I'm thinking of making type-5 the default for v29.4 as it has some advantages. 1) It allows Gerbicz error checking, 2) The residue remains the same as more and more factors are found for a Mersenne number - no double-check needed. 3) If a user decides to archive his PRP save files he can retest with a new known factor quite quickly.