Karim Belabas on Mon, 31 Aug 2009 08:52:43 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Questions about primality certification in Pari |
* 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] `