Re: Questions about primality certification in Pari

[Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>]

* 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]
`