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