Fixed miller-rabin test

[JASory]

The MR test was discovered to be flagging a composite as a prime, consequently the witness bases were adjusted to correct for this. Note that the primality of the new base (193) is purely thematic, many composites also satisfy the conditions.

[cuviper]

Note that the primality of the new base (193) is purely thematic, many composites also satisfy the conditions.

I don't have a good sense of how to evaluate the choice of 193. Do you have citations for this?

Wikipedia says this: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases

Using the work of Feitsma and Galway enumerating all base 2 pseudoprimes in 2010, this was extended (see OEIS: A014233), with the first result later shown using different methods in Jiang and Deng:[11]

• if n < 3,825,123,056,546,413,051, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, and 23.
• if n < 18,446,744,073,709,551,616 = 2⁶⁴, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37.

So that does also say that our current test up to prime 23 is not sufficient for all u64. (Ah, and that value that limits 23 is exactly the counterexample in your new test.) But if we're going to choose new values, I'd like some formal basis for that, as I don't know how to check that myself.

[cuviper]

FWIW, here's another source that I've used for larger values in my own Project Euler code: https://primes.utm.edu/prove/prove2_3.html

[cuviper]

I merged and published #47 as a simpler update, without worrying about alternate bases for the moment.

We can leave this open if you want to explore alternatives further, but I'd like some explanation, or at least a citation. For example, while changing that 23 to 193 might solve that particular composite test case, I don't know if that introduces any other errors.

[JASory]

I merged and published #47 as a simpler update, without worrying about alternate bases for the moment.

We can leave this open if you want to explore alternatives further, but I'd like some explanation, or at least a citation. For example, while changing that 23 to 193 might solve that particular composite test case, I don't know if that introduces any other errors.

It doesn't.

A correct deterministic Miller-Rabin test requires the following.

1. A correct strong fermat test (verified by inspection)
2. All primes are mapped to coprime bases (verified by inspection again)

3. All composites are filtered.

   In the last case as base-2 is applied to all integers, we can use Feitsma and Galway's table of 2-PRPs (weak) and select bases that filter them all. While I do not have the computational resources to fully verify F&G's table, it is almost surely correct based on statistical tests I've run and the fact that the upper 2-SPRP distribution closely mirrors an extension I have computed. It is also the basis for nearly every deterministic 64-bit primality test (this is why they all have base-2, 2 is actually only marginally stronger than average (except in very special configurations)).

FYI, while a jacobi & 2-SPRP test is terrific it is effectively deterministic which is why in ENT (number-theory) I add a random base (in addition to a bunch of other fast pseudoprime checks).

[JASory]

But if we're going to choose new values, I'd like some formal basis for that, as I don't know how to check that myself.

I simply made the minimal change necessary for correctness. The practice of using sequential primes is primarily historical as the finding the first sprp to a set of sequential prime bases is much faster than searching for the minimum bases. (See Sorenson and Webster's paper). Due to F&G's table this is no longer strictly necessary (it hasn't been for over a decade but it's not that well known). One can construct a far smaller set of bases utilizing specially weighted hashing ( I use only 2 checks although computing it is non-trivial and requires other fast form checks, note that using fewer than two checks requires enumerating all pseudoprimes and very large memory requirements so 2 is effectively optimal).

[cuviper]

I've experimented with different tables for reducing the number of guaranteed witnesses, but in benchmarking some of my own use cases, it seems to make things worse on average. My hunch is that it ends up having to test more witnesses to rule out a lot of composites, but I'm not sure about that. For now, I think we should leave it alone.

[JASory]

I've experimented with different tables for reducing the number of guaranteed witnesses, but in benchmarking some of my own use cases, it seems to make things worse on average. My hunch is that it ends up having to test more witnesses to rule out a lot of composites, but I'm not sure about that.

This is absolutely not true. Some basic statistical analysis will tell you why. Base 2 pseudoprimes are far rarer than primes, meaning that anything greater than a single sprp test is mostly wasted on primes, not composites. This is almost certainly erroneous benchmarking, as there is no mathematical explanation for this claimed phenomenon (and is even completely contradictory to my own benchmarks).

[cuviper]

If you want to submit new values with benchmarks to back it up, I'll take a look.

[JASory]

If you want to submit new values with benchmarks to back it up, I'll take a look. I'm not really interested in that at this point. I've largely stopped collaborating with professional software developers, since it results in tedious arguments of disbelief like this. If you are really concerned about benchmarks then check red-primality (which uses a nearly identical approach as yours simply with a smaller set of bases), or librapidrust which was a shitty implementation of Forisek-Jancina's table for 3 base checks (written by myself for @NervousNullPtr). Both of those roughly use the same arithmetic and are directly comparable to primal-check, except with fewer bases. Fewer base checks being more efficient is a clearly true statement.

[cuviper]

I don't think it's too much to ask you to justify this, but of course you're not obligated either.

I put together a quick test myself, using a PRNG to generate 1000 odd u64 as input:

test lib_rapid     ... bench:     140,749 ns/iter (+/- 1,951)
test primal_check  ... bench:     514,546 ns/iter (+/- 6,936)
test red_primality ... bench:     569,640 ns/iter (+/- 5,150)

So red-primality came out 10% slower, despite a smaller set of bases, but it's possible there are other implementation differences holding it back. On the other hand, lib_rapid won by a long shot. (I did not verify any of their results.)

Code

```rust #![feature(test)] #[cfg(test)] extern crate test; #[cfg(test)] fn run(b: &mut test::Bencher, f: impl Fn(&u64) -> bool) { use rand::prelude::*; b.iter(|| { let mut rng = StdRng::seed_from_u64(0x1111_2222_3333_4444); std::iter::repeat_with(move || rng.gen::() | 1) .take(1000) .filter(&f) .count() }); } #[bench] fn primal_check(b: &mut test::Bencher) { run(b, |&n| primal_check::miller_rabin(n)); } #[bench] fn red_primality(b: &mut test::Bencher) { run(b, |&n| red_primality::is_u64_prime(n)); } #[bench] fn lib_rapid(b: &mut test::Bencher) { use lib_rapid::math::primes::Primality; run(b, |n| Primality::is_prime(n)); } ```

[JASory]

I don't think it's too much to ask you to justify this, but of course you're not obligated either.

I put together a quick test myself, using a PRNG to generate 1000 odd u64 as input:

test lib_rapid     ... bench:     140,749 ns/iter (+/- 1,951)
test primal_check  ... bench:     514,546 ns/iter (+/- 6,936)
test red_primality ... bench:     569,640 ns/iter (+/- 5,150)

So red-primality came out 10% slower, despite a smaller set of bases, but it's possible there are other implementation differences holding it back. On the other hand, lib_rapid won by a long shot. (I did not verify any of their results.) Code

That's because LibRapidRust only uses 3 bases, also 1000 is a tiny input. I use millions to billions when benchmarking. Additionally, using a PRNG incurs variable overhead (possibly even larger than primality testing), even if you compute it beforehand and store the random integers because moving memory also takes time. Iterating sequentially over 63-bit integers is a more practical benchmark, if you are concerned about values being skewed by any preliminary trial division then you can start from the primorial p#53 to get an even distribution of numbers coprime to the first primes; but this is a miniscule difference that is almost certainly outweighed by cpu speed variance.

FYI, you can use the recently published fermat-analysis library to justify the base selection or select your own bases. It's how I computed that value (193), all the values I use in my software, and verified most (all?) of the public primality testing libraries in rust.