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.