Fast and robust error checking on Proth/Pepin tests

[R. Gerbicz]

In the recent days an interesting new error checking for LL test appeared, using Jacobi test.
I'll give another type of checking, that detects much more errors (we see later what this means), but works on Proth numbers.

The conditions: we allow at most O(log(N)) more memory to solve this,
and at most eps overhead in running time (say eps~1/1000 or even smaller).
We are seeking a high performance error checking algorithm.

quick background::
N=k*2^n+1 is Proth number if k<2^n, the base=a proves the primality of N,
if a^((N-1/2)==-1 mod N. In special case k=1 for Fermat number, and in
    this case we can choose a=3 (Pepin).
(if we choose a=quadratic nonresidue modulo N, then a^((N-1)/2)!=-1 mod N proves that
N is composite).

the setup::

Let 
L=2000 constant (could even depend on n), and L2=L^2.
u(t)=(a^k)^(2^t) mod N  [0]
d(t)=u(0)*u(L)*u(2*L)*...*u(t*L) mod N  [1]
with this 
d(t+1)=d(t)*u((t+1)*L) mod N  [2]
but 
d(t+1)=u(0)*d(t)^(2^L) mod N  [3]
is also true.

in the standard Proth test (usually for small k), first we compute u(0), then with n-1 iterative squaring we get u(n-1), with this u(n-1)==-1 mod N iff N is prime (for a=quadratic nonresidue).

We store only the last term of the d sequence to use identity [2].
(when we compute the next term, then two terms are available,
after a possible check we can delete the last but one term).
and store u(0)=a^k mod N,
the last d[z],u[z], where z is divisible by L2.

At each L-th term of the d sequence we check the identity of [3], if this does not hold, then
we roll back, notice that it is also possible a computation error of the d
sequence, so if we would roll back too much (say 100 times to the same term)
then we just restart completely the computation.

At the last few squarings in u, we also force an error checking computation of [3]
(in that i when i is divisible by L and i+L>=n, this means only one extra checking of [3].)

This leaves all potential erros in the (at most) last L squarings in u,
or very unlikely errors earlier in u or d.

The overhead is n/L mulmods in [2] and
n/L2*L=n/L squaremod in [3] and n/L2 mulmods in [3].
so over the n-1 mulmods of the Proth test there is approx. n/1000 mulmods, if we count in mulmods everything. So the overhead is 0.1% in time.

And we see why we haven't checked all terms of d, we could do that, but in that case the overhead in the error checking would be n/L*(L+1)>n mulmods, and that is a lot, slightly more time what we spend on the Proth test squarings.

using PARI-GP:

myrand(r,N) returns s randomly from [0,N) for that s!=0 and s!=r (the r,s are in Z_N).

myrand(r,N)={local(tmp);while(1,tmp=random(N);if(tmp!=0&&lift(tmp+r)!=0,return(tmp+r)))}

we test N=k*2^n+1 Proth number, we use L at error checking, making errors in the i-th squaring
with 50% chance if errpos[i]!=0 (note that if we return to the same i multiple times,
then we choose the making error independently from the previous choices already done)
we use base=a, this is optional, if you don't give, we find one suitable,
(here a=quadratic nonresidue)
if printmsg!=0, then we print out some additional info,
the return value is (1+a^((N-1)/2) mod N), note that for prime the 
return value is zero. [not forget the 1+ in the formula]
If you would not give a Proth number or errpos's length is too small, then the
return value is (-1).

fncheck(N,n,L,errpos,a=0,printmsg=1)={
k=(N-1)/(2^n);
if(k>=2^n||k<1||type(k)!="t_INT",if(printmsg,print("Not a Proth number"));return(-1));
if(length(errpos)<n-1,if(printmsg,print("The errpos array's length should be at least n-1"));return(-1));
if(a==0,a=2;while(kronecker(a,N)==1,a+=1));
numerr=0;
L2=L^2;
u0=Mod(a,N)^k;
prev_d=u0;
saved_u=u0;
saved_d=u0;
saved_i=0;
i=0;res=u0;while(i<n-1,i+=1;
res=res^2;
if(errpos[i]&&random(2),res=myrand(res,N));
if(i%L==0,d=prev_d*res;set_d=0;
if(i%L2==0||(i%L==0&&i+L>=n),
if(d!=u0*prev_d^(2^L),
numerr+=1;if(printmsg,print("Found error at iteration=",i,", roll back to iteration=",saved_i));
i=saved_i;res=saved_u;prev_d=saved_d;set_d=1,
saved_i=i;saved_u=res;saved_d=d));
if(!set_d,prev_d=d)));
res+=1;
if(printmsg,if(lift(res)==0,print("N is prime."),print("N is composite."));
print("Number of errors (corrected and) detected=",numerr));
return(lift(res))}

Proof that this detects the error exactly if there is at most one error
in each L*L squarings block:
suppose that the first L*z squarings was good, and after that in one of the next L squarings there was
an error, we got r0!=u0^(e^(z+1)) mod N.
(for shortening here we write e=2^L)

prev_d=u0*u0^e*u0^(e^2)*u0^(e^3)*...*u0^(e^z)*r0*r0^e*...*r0^(e^w)
d =u0*u0^e*u0^(e^2)*u0^(e^3)*...*u0^(e^z)*r0*r0^e*...*r0^(e^w)*r0^(e^(w+1))

d/prev_d^e/u0=r0/u0^(e^(z+1))!=1 mod N, so
d!=u0*prev_d^e mod N, so it will be detected.

What about multiple errors in the same block, say
we've chosen N=13*2^82+1 Proth number and we make errors at i=11,13,43 (at each position
with 50% chance), used L=3, this means that i=11,13 are in the same L2=9 length block,
so with 25% chance we have double errors.
Make the test for 50 times:

N=13*2^82+1;n=82;
errpos=vector(n-1,i,0);errpos[11]=1;errpos[13]=1;errpos[43]=1;
cnt=0;for(h=1,50,cnt+=(fncheck(N,n,3,errpos,0,1)==0);print());
print("Found prime ",cnt," times.");

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=7

Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=1

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=9

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=6

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=9

Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=3

Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=1

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=10

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=5

Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=1

Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=5

Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=4

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=10

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=7

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=4

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=4

Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=3

N is prime.
Number of errors (corrected and) detected=0

Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=1

Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=4

Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=9

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=14

N is prime.
Number of errors (corrected and) detected=0

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=14

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=2

N is prime.
Number of errors (corrected and) detected=0

Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=3

N is prime.
Number of errors (corrected and) detected=0

Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=1

N is prime.
Number of errors (corrected and) detected=0

N is prime.
Number of errors (corrected and) detected=0

Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=1

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=10

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=5

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=4

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=4

N is prime.
Number of errors (corrected and) detected=0

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=9

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=4

Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=1

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=5

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
N is prime.
Number of errors (corrected and) detected=2

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=3

Found error at iteration=18, roll back to iteration=9
Found error at iteration=18, roll back to iteration=9
Found error at iteration=45, roll back to iteration=36
N is prime.
Number of errors (corrected and) detected=3

? Found prime 50 times.
?

Impressive, detected ALL errors, but this is still NOT the theory, though
we are not that far from the truth, lets repeat it with a much smaller number
and with more times to see the true error rate:
test N=5*2^13+1 Proth number; 1 million times (takes less than a minute), using L=2.

N=5*2^13+1;n=13;
errpos=vector(n-1,i,0);errpos[5]=1;errpos[7]=1;
sum(h=1,1000000,fncheck(N,n,2,errpos,0,0)==0)

%12 = 999947

so there were only 53 false results out of a million Proth test, that is an error rate 53/10^6 ~ 2/N, quite good.

making (at most) two errors, but in different blocks:

N=5*2^13+1;n=13;
errpos=vector(n-1,i,0);errpos[7]=1;errpos[10]=1;
sum(h=1,1000000,fncheck(N,n,2,errpos,0,0)==0)

%14 = 1000000

This is exactly the theory, all errors detected, because they were in different blocks.

We can easily generalize this error checking to the computation of a^(k*b^t) mod N (for small k). So we can use not only b=2.

[GP2]

More discussion of this idea and its possible applicability to Mersenne numbers is in the thread from the GPU Computing subforum.

[preda]

Originally Posted by R. Gerbicz

I'll give another type of checking, that detects much more errors (we see later what this means), but works on Proth numbers.

Robert, do you have some idea about an error checking for stage 1 of P-1?

e.g. if the P-1 is done with left-to-right binary exponentiation.

[kriesel]

Originally Posted by preda

Robert, do you have some idea about an error checking for stage 1 of P-1?

e.g. if the P-1 is done with left-to-right binary exponentiation.

He referred to a P-1 check in http://mersenneforum.org/showpost.ph...&postcount=110
See also http://mersenneforum.org/showpost.ph...&postcount=245 and http://mersenneforum.org/showpost.ph...&postcount=252

[R. Gerbicz]

Originally Posted by kriesel

He referred to a P-1 check in
See also http://mersenneforum.org/showpost.ph...&postcount=245

That is good.
Trivial extension and speedup: with Euler-Fermat theorem
a^e==a^(e mod eulerphi(d)) mod N, if gcd(a,N)=1, so if d|N is also true then
a^e==a^(e mod eulerphi(d)) mod d
should be also true, if it doesn't hold then we made an error.
Ofcourse here calculation of eulerphi(d) is easy, if you know the prime factorization of d.

Another small check: if e is even then (res / N) Jacobi symbol should be 1, where res=a^e mod N.
And if you make any multiplication error then you can discover it at any iteration with 50% chance.

Not counting these I don't know other checks. With a so general e, there is no such ladder trick, unfortunately.

ps. you could also post these to the longer error check's thread.

[kriesel]

There was also the generic error checking and detection thread http://www.mersenneforum.org/showthread.php?t=23467
and a P-1 specific Jacobi check thread http://www.mersenneforum.org/showthread.php?t=23470

[kriesel]

Originally Posted by preda

Robert, do you have some idea about an error checking for stage 1 of P-1?

e.g. if the P-1 is done with left-to-right binary exponentiation.

Not sure what you mean by left to right binary exponentiation.
Are you contemplating implementing P-1 in OpenCL?
If so, that would fill a prominent vacancy in the available software versus algorithm and computing platform grid in part one of the attachment at http://www.mersenneforum.org/showpos...91&postcount=2

[science_man_88]

Originally Posted by kriesel

Not sure what you mean by left to right binary exponentiation.
Are you contemplating implementing P-1 in OpenCL?
If so, that would fill a prominent vacancy in the available software versus algorithm and computing platform grid in part one of the attachment at http://www.mersenneforum.org/showpos...91&postcount=2

http://primes.utm.edu/glossary/xpage...entiation.html

[preda]

Thank you Robert and @kriesel! I'll think about these ideas. In my setup, I don't have a divisor d|N. Also the "up to 50%" of Jacobi seems weak IMO.

[preda]

Originally Posted by kriesel

Not sure what you mean by left to right binary exponentiation.
Are you contemplating implementing P-1 in OpenCL?

AKA top-to-bottom, starting with the most significant bit and moving towards the least significant.

Yes, contemplating. P-1 is magically attractive.

[kriesel]

Originally Posted by preda

... the "up to 50%" of Jacobi seems weak IMO.

It is what it is, including better than nothing. Two things can make it better than it first seems.

1) It appears, from a test I wrote and ran, that for the early stage 1 powering computations, before the mod Mp kicks in, induced errors don't matter. I think that means they are equivalent to starting with a different base, which is allowed.
2) For tasks like computing the exponent (or partial computation) to which to raise the base, the Jacobi of the product and the product of the Jacobis of the terms can be computed and compared twice, with different denominators, going from ~50% to ~75% detection. I tried denominators M31 and M61 and confirmed that. That's extensible to 3 or more, with diminishing returns and increasing overhead.

Whether the checking overhead is worthwhile depends on the error rate. Making error checking optional and the default might be good. One of the challenges with P-1 is it's hard to tell from the results produced whether it's all working correctly. Starting by error checking a new installation is advisable.