Karim Belabas on Mon, 31 Aug 2009 08:52:43 +0200

 Re: Questions about primality certification in Pari

```* William Hart [2009-08-31 05:29]:
> I have two questions, the second one regarding what may be a bug:

First a remark: your questions concern pari-stable, but you should
rather study the current SVN code which is a bit different (and more
efficient).

> Q1) In arith2.c we have:
>
> int
> BSW_isprime(GEN x)
> {
>   pari_sp av = avma;
>   long l, res;
>   GEN F, p, e;
>   if (BSW_isprime_small(x)) return 1;
>   F = auxdecomp(subis(x,1), 0);
>   l = lg(gel(F,1))-1; p = gcoeff(F,l,1); e = gcoeff(F,l,2); F=gel(F,
> 1);
>   if (cmpii(powgi(p, shifti(e,1)), x)<0)
>     res = isprimeSelfridge(mkvec2(x,vecslice(F,1,l-1))); /* half-
> smooth */
>   else if (BSW_psp(p))
>     res = isprimeSelfridge(mkvec2(x,F)); /* smooth */
>   else
>     res = isprimeAPRCL(x);
>   avma = av; return res;
> }
>
> long
> isprime(GEN x)
> {
>   return BSW_psp(x) && BSW_isprime(x);
> }
>
> Now the BSW_psp test I understand. It is the Bailley-Selfridge-
> Wagstaff "pseudo"primality test. When x is small enough (< 10^13 in
> pari) it is a guaranteed primality test.

Actually, < 10^15 is still fine. Should be doable (and useful) to extend this
further, to 2^64 say.

> Now if x is too large for this to guarantee you have a prime, pari then
> (as can be seen) does some factoring of x-1. The factors it finds, it
> sticks in the vector F. There might be some cofactor p^e. The first
> thing it does is check whether (p^e)^2 < x. This is because the
> Pocklington-Lehmer test can then be applied by looking for what Pari
> calls witnesses for each of the factors of x-1.

N.B. Current code needs partial factorization up to x^(1/3), instead of x^(1/2).

> It is the next step I do not understand. Supposing p^e >= sqrt(x) they
> it checks if p is a BSW_psp. Again it basically does a Pocklington-
> Lehmer test if this is the case (they call it a Selfridge test).
>
> Now my question is, why can one get away with a psp test rather than a full primality test here?

isprimeSelfridge receives an integer N and a list of pseudoprime
divisors of N-1. It performs the BSW test proving the primality of N
under the assumption that the pseudoprimes are true primes, then
calls the primality prover recursively to prove the assumption.

Notice that in the second case ( p^(2e) >= x ) we do stick p in the
vector of pseudoprimes. Its primality will then be certified.

In the first one (p^(2e) < x), the vecslice() removes p from the list.

This code looks fine to me.

> Q2) Also in arith2.c:
>
> int
> BSW_isprime_small(GEN x)
> {
>   long l = lgefint(x);
>   if (l < 4) return 1;
>   if (l == 4)
>   {
>     pari_sp av = avma;
>     long t = cmpii(x, u2toi(0x918UL, 0x4e72a000UL)); /* 10^13 */
>     avma = av;
>     if (t < 0) return 1;
>   }
>   return 0;
> }
>
> This code is used to determine if the BSW test is sufficient to guarantee primality, by testing whether the number is < 10^13.
>
> But on a 64 bit machine, doesn't this always return 1 when x is a single limb integer? It looks like it was written for a 32 bit machine. Shouldn't it be checking if l == 3 on a 64 bit machine then checking if x < 10^13.
> My apologies if that is not right.

You are right: this one is a bug.

Thanks for reporting this !

K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite Bordeaux 1          Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`

```