William Hart on Mon, 31 Aug 2009 05:29:09 +0200

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Questions about primality certification in Pari

I have two questions, the second one regarding what may be a bug:

Q1) In arith2.c we have:

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, 
  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 */ 
    res = isprimeAPRCL(x); 
  avma = av; return res; 

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. 

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. 

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?

Q2) Also in arith2.c:

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. I tried to trace the source code through, but got into a complete muddle. At least when I change that to what I think it should be, it still works, but goes much slower.