n_isprime() returning true for composites!!! about flint

[danaj]

It looks like SymPy found this in Jan 2010: http://code.google.com/p/sympy/issues/detail?id=1789, and their solution was add a giant list of exceptions. I have a SymPy branch that I worked on last year that replaces all that with small deterministic M-R followed by BPSW. It is a huge help for large numbers (both in correctness for known examples as well as speed). Bit-rotting at the moment :(.

Could someone else verify they get this result, to make sure this isn't me just doing incorrect usage?

#include <stdio.h>
#include <stdlib.h>
#include "flint.h"
#include "ulong_extras.h"
int main(void) {
   mp_limb_t n = UWORD(2007193456621);
   printf("%lu is %s\n", n, n_is_prime(n) ? "prime" : "composite");
   return 0;
}

Denis said he did not, but I see it on both Ubuntu and Fedora's FLINT libraries as well as 2.4.3. It's going through the 5-base test for these numbers for me (2007193456621, 2744715551581, 9542968210729, 17699592963781, 19671510288601, etc.).

Revisiting the tests wouldn't hurt. For 32-bit, the {31,73} and {2,7,61} tests should cover all cases, and perhaps would be clearer to just separate out the !FLINT64 and FLINT64 cases instead of mixing them. For 64-bit it helps to have the M-R function able to deal with bases >= n, which lets you use any of the appspot base sets -- e.g. the 1-base (n < 341531) and 2-base (n < 1050535501) sets, then measure whether the 3-base (n < 350269456337) is faster in that range than BPSW. Alternately use step-up sets (so the bases are always smaller than the numbers being tested).

I have not tested with FLINT, but the usual way to handle large bases, given base a and a number n:

  if (a >= n)  a % = n;
  if (a <= 1 || a == n-1) continue;   /* base out of bounds, consider it passing */

Any solution will need to be tested. It isn't hard to exhaustively test the 1- and 2-base versions, and not hard to test anything using base 2. The large 3-base solutions are difficult.

[wbhart]

I can verify that it claims prime on my machine. Definitely a bug.

On 9 May 2014 00:35, Dana Jacobsen [email protected] wrote:

It looks like SymPy found this in Jan 2010:
http://code.google.com/p/sympy/issues/detail?id=1789, and their solution
was add a giant list of exceptions. I have a SymPy branch that I worked on
last year that replaces all that with small deterministic M-R followed by
BPSW. It is a huge help for large numbers (both in correctness for known
examples as well as speed). Bit-rotting at the moment :(.

Could someone else verify they get this result, to make sure this isn't me
just doing incorrect usage?

#include <stdio.h>#include <stdlib.h>#include "flint.h"#include "ulong_extras.h"int main(void) {
mp_limb_t n = UWORD(2007193456621);
printf("%lu is %s\n", n, n_is_prime(n) ? "prime" : "composite");
return 0;}

Denis said he did not, but I see it on both Ubuntu and Fedora's FLINT
libraries as well as 2.4.3. It's going through the 5-base test for these
numbers for me (2007193456621, 2744715551581, 9542968210729,
17699592963781, 19671510288601, etc.).

Revisiting the tests wouldn't hurt. For 32-bit, the {31,73} and {2,7,61}
tests should cover all cases, and perhaps would be clearer to just separate
out the !FLINT64 and FLINT64 cases instead of mixing them. For 64-bit it
helps to have the M-R function able to deal with bases >= n, which lets you
use any of the appspot base sets -- e.g. the 1-base (n < 341531) and 2-base
(n < 1050535501) sets, then measure whether the 3-base (n < 350269456337)
is faster in that range than BPSW. Alternately use step-up sets (so the
bases are always smaller than the numbers being tested).

I have not tested with FLINT, but the usual way to handle large bases,
given base a and a number n:

if (a >= n) a % = n;
if (a <= 1 || a == n-1) continue; /* base out of bounds, consider it passing */

Any solution will need to be tested. It isn't hard to exhaustively test
the 1- and 2-base versions, and not hard to test anything using base 2. The
large 3-base solutions are difficult.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/69#issuecomment-42614829
.

[danaj]

I made a branch with changes:
https://github.com/danaj/flint2/tree/issue69-isprime
Let me know if I should do a pull request then discuss.

• I added large base support to is_strong_probabprime_precomp.c. Arguably the same should be done in the preinv case, but I'm not using that so it wasn't necessary.
• I added a simple divisibility check to is_prime.c. It could be removed or adjusted as desired. It speeds up the generic case.
• I make the decision of M-R vs. BPSW in is_probabprime, not in is_prime. The cost is an extra function call, the advantage is not having a magic number in is_prime.c.
• I removed the n_is_perfect_power235(n)) call. I believe calling it is a performance loss on average, but maybe there is some other reason it was there?
• I have not tested the !FLINT64 code on a non-64-bit machine.
• BPSW was very slightly faster than the 3-base M-R test for me, and given that testing its correctness is much easier than the 3-base test, I went for BPSW. I left the 3-base test code ifdef'd out -- it can be removed or left.

Let me know what you think.

I should add some of the composites to the test suite. In t-is_prime.c or t_is_probabprime.c or both?

[wbhart]

Hi Dana,

thanks for this.

On 13 May 2014 10:41, Dana Jacobsen [email protected] wrote:

I made a branch with changes:
https://github.com/danaj/flint2/tree/issue69-isprime
Let me know if I should do a pull request then discuss.

• I added large base support to is_strong_probabprime_precomp.c.
  Arguably the same should be done in the preinv case, but I'm not using that
  so it wasn't necessary.

Yes, I agree it should be added in the preinv code too. The difference is
that the precomp code uses a double precomputed inverse, so is only good
for moduli up to 53 bits. That's likely good enough for the purposes for
which we use it, but for consistency the preinv case should do the same.

In the precomp case, the a%=n could probably be replaced with a call to
n_mod2_precomp (or n_mod_precomp restricted to the cases where it is
actually applicable). I doubt the compiler ever uses a precomputed inverse
in this code, even for constant bases.

• I added a simple divisibility check to is_prime.c. It could be
  removed or adjusted as desired. It speeds up the generic case.

This looks good, given that the generic use for this code is to check
whether something is prime or not, rather than to prove something prime
(assuming it very likely is).

• I make the decision of M-R vs. BPSW in is_probabprime, not in
  is_prime. The cost is an extra function call, the advantage is not having a
  magic number in is_prime.c.

Ok.

• I removed the n_is_perfect_power235(n)) call. I believe calling it
  is a performance loss on average, but maybe there is some other reason it
  was there?

I don't recall why it was there. From memory, there was a point at which
flint was certifying perfect cubes as prime. I seem to recall it was added
at that stage. I probably should have added a regression test.

• I have not tested the !FLINT64 code on a non-64-bit machine.

No problem. It will be tested at the release. I propose we do a very quick
point release to fix this problem fairly soon. I will make sure it gets
tested on a 32 bit machine (if I can find one outside a museum).

• BPSW was very slightly faster than the 3-base M-R test for me, and
  given that testing its correctness is much easier than the 3-base test, I
  went for BPSW. I left the 3-base test code ifdef'd out -- it can be removed
  or left.

That's quite surprising. But also good news.

We tend to not leave commented out code in (there are a handful of
exceptions). I would simply remove it.

Let me know what you think.

I think it is great. Thanks very much for the contribution. Somehow I think
you were the most qualified to fix these issues.

I should add some of the composites to the test suite. In t-is_prime.c or
t_is_probabprime.c or both?

I think it is a good idea. Just a static table should be enough. You may
need to do the !FLINT64 thing if the constants are too large. Note that to
support windows, constants have to be defined using the WORD/UWORD macros
(see flint.h -- I just noticed these are not discussed in
code_conventions.txt; I will add a ticket to fix that).

Bill.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/69#issuecomment-42930576
.

[wbhart]

I forgot to add: please do issue a pull request, if you feel you have the
time.

Bill.

On 13 May 2014 12:17, Bill Hart [email protected] wrote:

Hi Dana,

thanks for this.

On 13 May 2014 10:41, Dana Jacobsen [email protected] wrote:

I made a branch with changes:
https://github.com/danaj/flint2/tree/issue69-isprime
Let me know if I should do a pull request then discuss.

• I added large base support to is_strong_probabprime_precomp.c.
  Arguably the same should be done in the preinv case, but I'm not using that
  so it wasn't necessary.

Yes, I agree it should be added in the preinv code too. The difference is
that the precomp code uses a double precomputed inverse, so is only good
for moduli up to 53 bits. That's likely good enough for the purposes for
which we use it, but for consistency the preinv case should do the same.

In the precomp case, the a%=n could probably be replaced with a call to
n_mod2_precomp (or n_mod_precomp restricted to the cases where it is
actually applicable). I doubt the compiler ever uses a precomputed inverse
in this code, even for constant bases.

• I added a simple divisibility check to is_prime.c. It could be
  removed or adjusted as desired. It speeds up the generic case.

This looks good, given that the generic use for this code is to check
whether something is prime or not, rather than to prove something prime
(assuming it very likely is).

• I make the decision of M-R vs. BPSW in is_probabprime, not in
  is_prime. The cost is an extra function call, the advantage is not having a
  magic number in is_prime.c.

Ok.

• I removed the n_is_perfect_power235(n)) call. I believe calling it
  is a performance loss on average, but maybe there is some other reason it
  was there?

I don't recall why it was there. From memory, there was a point at which
flint was certifying perfect cubes as prime. I seem to recall it was added
at that stage. I probably should have added a regression test.

• I have not tested the !FLINT64 code on a non-64-bit machine.

No problem. It will be tested at the release. I propose we do a very
quick point release to fix this problem fairly soon. I will make sure it
gets tested on a 32 bit machine (if I can find one outside a museum).

• BPSW was very slightly faster than the 3-base M-R test for me, and
  given that testing its correctness is much easier than the 3-base test, I
  went for BPSW. I left the 3-base test code ifdef'd out -- it can be removed
  or left.

That's quite surprising. But also good news.

We tend to not leave commented out code in (there are a handful of
exceptions). I would simply remove it.

Let me know what you think.

I think it is great. Thanks very much for the contribution. Somehow I
think you were the most qualified to fix these issues.

I should add some of the composites to the test suite. In t-is_prime.c or
t_is_probabprime.c or both?

I think it is a good idea. Just a static table should be enough. You may
need to do the !FLINT64 thing if the constants are too large. Note that to
support windows, constants have to be defined using the WORD/UWORD macros
(see flint.h -- I just noticed these are not discussed in
code_conventions.txt; I will add a ticket to fix that).

Bill.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/69#issuecomment-42930576
.