Primality testing

[Bill Hart]

Are there any primality testing experts on the list. I'm working on
some code for this and I've noticed Pari does some clever things. But
one thing I do not understand.

Here is the relevant piece of pari code:

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
then check if p is a BSW_psp. Again they basically do a Pocklington-
Lehmer test if this is the case (they call it a Selfridge test). ****

Obviously if x-1 cannot be factored in an appropriate way, they fall
back to a hard core primality test, in this case APRCL. I have
something else implemented in its place in FLINT which should be
faster and is certainly easier to implement.

Anyhow, I'd like to do something similar to what they've done at the
step I've labelled **** above. Does anyone know why this works, or why
it is permissible.

Obviously there is no reason why we can't have x = 2pq+1 for some
large primes p and q. Then pq is a factor x-1 which BSW_psp might
declare prime. It is not clear to me that the Pocklington-Lehmer test
works in this case.

Bill.

[Bill Hart]

Another question. The code which determines if a BSW pseudoprimality
test is sufficient or not to absolutely guarantee primality is:

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;
}

But on a 64 bit machine, I'm pretty sure this always returns 1. It
looks like it was written for a 32 bit machine. That should be
checking if l == 3 on a 64 bit machine then checking if x < 10^13.

The first limb of a gen stores the type, the second the number of
limbs and some other bits of info and the third limb is the first limb
of actual data for an int.

I'm pretty sure that is a bug. So pari is not correctly certifying
single limb integers prime. No wonder it is so damned fast at
certifying primality.

Bill.

[William Stein]

On Sun, Aug 30, 2009 at 7:19 PM, Bill Hart <goodwi...@googlemail.com> wrote:

Another question. The code which determines if a BSW pseudoprimality
test is sufficient or not to absolutely guarantee primality is:

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;
}

But on a 64 bit machine, I'm pretty sure this always returns 1. It
looks like it was written for a 32 bit machine. That should be
checking if l == 3 on a 64 bit machine then checking if x < 10^13.

The first limb of a gen stores the type, the second the number of
limbs and some other bits of info and the third limb is the first limb
of actual data for an int.

I'm pretty sure that is a bug. So pari is not correctly certifying
single limb integers prime. No wonder it is so damned fast at
certifying primality.

Argh!  That's really disappointing (if true).

Any change you could also send these emails to the pari development list, where they also naturally belong?
 

--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

[Bill Hart]

I want to verify first before reporting bugs. I'll download their
latest stable source, make the changes (if they haven't already) and
check that my "fix" works. Pari is pretty complex, and it would be
easy to be wrong.

Bill.

On 31 Aug, 03:40, William Stein <wst...@gmail.com> wrote:

[Bill Hart]

I checked that what I thought should work, does in fact work, and that
the bugs are in the latest stable release code.

It's ok, my "fix" only slows primality testing down by a factor of 92
though. I'm glad however, as I had just spent days trying to figure
out why FLINT was 4 times slower than Pari. It's now 2 times slower,
but I was beginning to despair. I think with a lot of extra code I
could make it 30% faster again, but when I looked at what Pari was
actually doing, I kept thinking, "wait a minute, Pari doesn't even
have *that* trick, how am I not faster than them."

Oh well, I could still turn out to be wrong. It wouldn't be the first
time. So probably no need to become too concerned until I verify it
with the Pari developers. To be honest, I don't know much about the
internal data format or the lgefint macro in Pari, so I could easily
be mistaken.

Bill.

[John Cremona]

2009/8/31 Bill Hart <goodwi...@googlemail.com>:

>
> I checked that what I thought should work, does in fact work, and that
> the bugs are in the latest stable release code.
>
> It's ok, my "fix" only slows primality testing down by a factor of 92

Is that a typo for 2?

I would talk to Karim now, not wait until you have sorted everything out.

John

[Bill Hart]

I checked with Karim. There are three important points:

* Firstly that lgefint == 4 is a bug. It should be 3 on a 64 bit
system (in both places) and should check that n < 10^13. Fixing this
is a substantial slowdown (no, the 92 was not a typo for 2, but see
below).

* Pari svn is more efficient at primality testing than pari stable. Do
we use pari svn yet in Sage?

* Karim believes the first issue above is not a bug. He claims that
the Selfridge primality test passes pseudoprime factors of n-1 and
then only checks later if they are actual primes. I haven't worked out
where it does this, but maybe the svn code is different or maybe I've
missed something. I'll look into it more carefully.

Bill.

On 31 Aug, 12:13, John Cremona <john.crem...@gmail.com> wrote:
> 2009/8/31 Bill Hart <goodwillh...@googlemail.com>:

[Bill Hart]

Just a quick followup, Karim is of course correct about the first
thing I mentioned above not being a bug. The pseudoprimes are indeed
later verified to be prime, inside the pocklington lehmer code.

So only the second issue above is actually a bug. According to my
timings fixing it makes pari-stable 71 times slower at certifying
primality of random probable primes that fit in a limb with random
numbers of bits from 2 to 64. Of course pari-svn is slightly faster
and Karim may well find some easy speedups.

Note that the nextprime function in pari does not certify primality by
default, by the way. Like MPIR it merely returns the next probable
prime (hence the change of name in MPIR).

Bill.