Re: Questions about primality certification in Pari

[William Hart <hart_wb@yahoo.com>]

--- On Mon, 8/31/09, Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr> wrote:

> From: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
> Subject: Re: Questions about primality certification in Pari
> To: "William Hart" <hart_wb@yahoo.com>
> Cc: pari-dev@list.cr.yp.to
> Date: Monday, August 31, 2009, 3:05 PM
> * 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.

In fact I had a student work on some code to try and extend this. So far it'll still take 1000 years even on a cluster. There is a project under way somewhere to extend this though. I think it is headed by Galway.

There are certainly figures available to take you to 10^16 currently, and I believe unpublished figures to 10^17, though I haven't been able to verify the latter.

> 
> > 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).

I'm told there is an x^(1/4) version too, but it's heresay, I don't have a reference for you. I've really only just begun to look at primality proving seriously. Another interesting algorithm is the pseudosquares algorithm, see Theorem 2.7 of "Some
   results on pseudosquares" by Lukes, Patterson and Williams,
   Math. Comp. vol 65, No. 213. pp 361-372. See
   http://www.ams.org/mcom/1996-65-213/S0025-5718-96-00678-3/S0025-5718-96-00678-3.pdf

In practice I found this slower than PL, but much simpler than APRCL for small integers. I have not done timings against APRCL as I don't have an implementation in my own code yet. Of course APRCL is for general integers, not just small ones.

> 
> > 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.
> 

That makes sense. And yes, I see that now. Clever trick!!

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

Got it.

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

Yes, PL doesn't need it. 

> This code looks fine to me.

Great, thanks for clarifying. That's an excellent trick. I suppose in practice the BSW pseudoprimes are always prime anyway, so it never ends up saving time though.

> 
> > 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 !

No problem.

Bill.