Miller Rabin Primality test accuracy

Question

I know the Miller–Rabin primality test is probabilistic. However I want to use it for a programming task that leaves no room for error.

Can we assume that it is correct with very high probability if the input numbers are 64-bit integers (i.e. long long in C)?

Invulnerable

It seems to me that long long is defined as having a minimum of 64 bits. – Albright

Answer 1

Miller–Rabin is indeed probabilistic, but you can trade accuracy for computation time arbitrarily. If the number you test is prime, it will always give the correct answer. The problematic case is when a number is composite, but is reported to be prime. We can bound the probability of this error using the formula on Wikipedia: If you select k different bases randomly and test them, the error probability is less than 4^-k. So even with k = 9, you only get a 3 in a million chance of being wrong. And with k = 40 or so it becomes ridiculously unlikely.

That said, there is a deterministic version of Miller–Rabin, relying on the correctness of the generalized Riemann hypothesis. For the range u up to 2⁶⁴, it is enough to check a = 2, 3, 5, 7, 11, 13, 17, 19, 23. I have a C++ implementation online which was field-tested in lots of programming contests. Here's an instantiation of the template for unsigned 64-bit ints:

bool isprime(uint64_t n) { //determines if n is a prime number
    const int pn = 9, p[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
    for (int i = 0; i < pn; ++i)
        if (n % p[i] == 0) return n == p[i];
    if (n < p[pn - 1]) return 0;
    uint64_t s = 0, t = n - 1;
    while (~t & 1)
        t >>= 1, ++s;
    for (int i = 0; i < pn; ++i) {
        uint64_t pt = PowerMod(p[i], t, n);
        if (pt == 1) continue;
        bool ok = 0;
        for (int j = 0; j < s && !ok; ++j) {
            if (pt == n - 1) ok = 1;
            pt = MultiplyMod(pt, pt, n);
        }
        if (!ok) return 0;
    }
    return 1;
}

PowerMod and MultiplyMod are just primitives to multiply and exponentiate under a given modulus, using square-and-{multiply,add}.

Endsley

What's with the template? The question is tagged c. – Hilliary

Reference / Reasoning? – Sindhi

@CristianCiupitu Just treat is as pseudo code. Niklas can hardly post a biginteger library here. It's obvious what PowerMod and MultiplyMod are, but it's annoying to implement them and doesn't help with the understanding of the algorithm. – Linoel

@JanDvorak Click on the first link of my answer? – Endsley

@CristianCiupitu It's C++ code copied from a working implemenation, but you should be able to manually instantiate the template by replacing every occurence of T by uint64_t or unsigned long long. – Endsley

Also the implementations of MultiplyMod and PowerMod can be looked up in the complete source file which I linked – Endsley

@Brett it's optimized for amount of characters to type during a competition, because you are not allowed at ICPC to access the Internet. So yeah, maybe we have different design goals ;) – Endsley

@Brett MulMod is just standard square/multiply with addition instead of multiplication. I think last time I checked __int128 was pretty slow, but maybe that has improved since – Endsley

@NiklasB. - I wouldn't worry about it. 64-bit processors can do __int128 ok with gcc / clang, but it still doesn't qualify as a portable solution. Unless you hack something like gmp's longlong.h or you're a really big Knuth fan! – Herv

BTW - you don't need (9) base tests, Sinclair has a (7) base test. It's been verified without appeals to the GRH, using GPUs, and other techniques that mad / obsessed people do with their time:) – Herv

@BrettHale sure, but that's not the first 7 primes ;) They are easier to type without making mistakes :P – Endsley

I was looking at some older questions, and it occurs to me that your 'first 9 primes' have a product that fits in a 32-bit unsigned int. That's ~ 70% of odd candidates eliminated using a (potentially much cheaper) gcd test. It's a close call, given the effectiveness of the 2-SPRP test, but still an interesting result. With a bit of work, it's substantially more effective than 9 modulo tests... I think. – Herv

Accuracy formula is at en.wikipedia.org/wiki/… – Mallorca

Answer 2

For n < 2^64, it is possible to perform strong-pseudoprime tests to the seven bases 2, 325, 9375, 28178, 450775, 9780504, and 1795265022 and completely determine the primality of n; see http://miller-rabin.appspot.com/.

A faster primality test performs a strong-pseudoprime test to base 2 followed by a Lucas pseudoprime test. It takes about 3 times as long as a single strong-pseudoprime test, so is more than twice as fast as the 7-base Miller-Rabin test. The code is more complex, but not dauntingly so.

I can post code if you're interested; let me know in the comments.

Schubert

For me, this set of bases works for all n < 2^32 except that it fails for n=4033 and n=4681, because these two numbers are pseudoprimes to bases 2 and 235. As I understand it, with Miller–Rabin, you only test bases less than n, so in the case of n=4033 and n=4681, it means testing only bases 2 and 235. How is it possible that this escaped the attention of Jim Sinclair (cited as discoverer of this 7-base set)? Or am I missing something here? – Pyretic

It's 325, not 235. And you test all seven bases. – Schubert

Ah yes, 235, not 325, thank you. Slippery fingers above. My code does say 325 (I used cut & paste to insert the set). Someone else pointed out to me that even though you are supposed to test all 7 bases, you are actually supposed to ignore bases that are multiples of n (e.g., when testing for n = 5, you are supposed to ignore bases 235, 9375, and 450775, which I wasn't doing). I just modified my program to ignore a when a ≡ 0 (mod n), and it appears to be doing the right thing now. Actually, I return true (is probable prime) if a ≡ 0 (mod n) rather than ignoring a. – Pyretic

An overwhelming proportion of composite candidates will fail the 2-SPRP test in any case. e.g., only 63 of all 24-bit odd composites pass the 2-SPRP test. I think you underplay the complexity of a Lucas pseudoprime test, and Jacobi sums. If we were talking about cryptographic integer sizes I might agree, but I think the idea that this is 'faster' on average for this range is misleading. – Herv

See my blog. – Schubert

Answer 3

In each iteration of Miller-Rabin you need to choose a random number. If you are unlucky this random number doesn't reveal certain composites. A small example of this is that 2^341 mod 341 = 2, passing the test

But the test guarantees that it only lets a composite pass with probability <1/4. So if you run the test 64 times with different random values, the probability drops below 2^(-128) which is enough in practice.

You should take a look at the Baillie–PSW primality test. While it may have false positives, there are no known examples for this and according to wikipedia has been verified that no composite number below 2^64 passes the test. So it should fit your requirements.

Linoel

Can you quote the algorithm here? Also, a link to wikipedia would be nice. Then I'll consider this a good answer. – Sindhi

@JanDvorak Markdown didn't accept the dash in the wikipedia link. – Linoel

Thanks @Linoel ,atleast you answered !! – Invulnerable

Answer 4

There are efficient deterministic variants of the MR test for 64-bit values - which do not rely on the GRH - having been exhaustively tested by exploiting GPUs and other known results.

I've listed the pertinent sections of a C program I wrote that tests the primality of any 64-bit value: (n > 1), using Jaeschke's and Sinclair's bases for the deterministic MR variant. It makes use of gcc and clang's __int128 extended type for exponentiation. If not available, an explicit routine is required. Maybe others will find this useful...

#include <inttypes.h>

/******************************************************************************/

static int sprp (uint64_t n, uint64_t a)
{
    uint64_t m = n - 1, r, y;
    unsigned int s = 1, j;

    /* assert(n > 2 && (n & 0x1) != 0); */

    while ((m & (UINT64_C(1) << s)) == 0) s++;
    r = m >> s; /* r, s s.t. 2^s * r = n - 1, r in odd. */

    if ((a %= n) == 0) /* else (0 < a < n) */
        return (1);

    {
        unsigned __int128 u = 1, w = a;

        while (r != 0)
        {
            if ((r & 0x1) != 0)
                u = (u * w) % n; /* (mul-rdx) */

            if ((r >>= 1) != 0)
                w = (w * w) % n; /* (sqr-rdx) */
        }

        if ((y = (uint64_t) u) == 1)
            return (1);
    }

    for (j = 1; j < s && y != m; j++)
    {
        unsigned __int128 u = y;
        u = (u * u) % n; /* (sqr-rdx) */

        if ((y = (uint64_t) u) <= 1) /* (n) is composite: */
            return (0);
    }

    return (y == m);
}

/******************************************************************************/

static int is_prime (uint64_t n)
{
    const uint32_t sprp32_base[] = /* (Jaeschke) */ {
        2, 7, 61, 0};

    const uint32_t sprp64_base[] = /* (Sinclair) */ {
        2, 325, 9375, 28178, 450775, 9780504, 1795265022, 0};

    const uint32_t *sprp_base;

    /* assert(n > 1); */

    if ((n & 0x1) == 0) /* even: */
        return (n == 2);

    sprp_base = (n <= UINT32_MAX) ? sprp32_base : sprp64_base;

    for (; *sprp_base != 0; sprp_base++)
        if (!sprp(n, *sprp_base)) return (0);

    return (1); /* prime. */
}

/******************************************************************************/

Note that the MR (sprp) test is slightly modified to pass values on an iteration where the base is a multiple of the candidate, as mentioned in the 'remarks' section of the website

---

Update: while this has fewer base tests than Niklas' answer, it's important to note that the bases: {3, 5, 7, 11, 13, 17, 19, 23, 29} provide a cheap test that allows us to eliminate candidates exceeding: 29 * 29 = 841 - simply using the GCD.

For (n > 29 * 29), we can clearly eliminate any even value as prime. The product of the small primes: (3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29} = 3234846615, fits nicely in a 32-bit unsigned value. A gcd(n, 3234846615) is a lot cheaper than a MR test! If the result is not (1), then (n) > 841 has a small factor.

Merten's (?) theorem suggests that this simple gcd(u64, u64) test eliminates ~ 68% of all odd candidates (as composites). If you're using M-R to search for primes (randomly or incrementally), rather than just a 'one-off' test, this is certainly worth while!

Herv

Two out of four answers is not 'most'; and the other two answers mentioned the 'deterministic MR' years ago. This means that your opening remarks are thoroughly misleading for pundits just browsing this topic. However, your answer does offer something new because you provided an excellent C rendition where the others offered only links or C++, and hence the +1 (mostly because I'm a sucker for clean, efficient professional code). – Efthim

@Efthim - We'll have to agree to disagree. The idea of using a sledgehammer like the BPSW or Lucas tests is misleading - as well as being outside the scope of the OPs request. The implicit complications of Jacobi sum implementations - not a trivial exercise in itself - and the fact that SPRP-2 eliminates the overwhelming majority of composites in any case; not to mention factual errors, like needing 9 bases, or its reliance on the Riemann hypothesis, which is only relevant as a generalised bound - not as an exhaustive test due to Sinclair and others. – Herv

You are right - after taking a closer look at Niklas B's post I do have to concede the majority (Niklas was thinking in the right direction but unlike user448810 he didn't quite get there). And thanks again for posting your code - a short, tight piece like that can give a lot more useful information and insights than reams of prose. – Efthim

@Efthim - For what it's worth, my opening remarks are a bit strong, and my defence is probably a little pedantic. It's always easy to misinterpret tone with writing. I should have been a bit more diplomatic:) – Herv

Answer 5

Your computer is not perfect; it has a finite probability of failing in such a way as to produce an incorrect result to a calculation. Providing the probability of the M-R test giving a false result is greatly less than the probability of some other computer failure, then you are fine. There is no reason to run the M-R test for less than 64 iterations (a 1 in 2^128 chance of error). Most examples will fail in the first few iterations, so only the actual primes will be thoroughly tested. Use 128 iterations for a 1 in 2^256 chance of error.

Tribunal

Unfortunately, this answer is not informed about the reality of deterministic M-R tests in this range. The probalistic logic is roughly correct, but not really relevant to a determininstc test in a 2^64 range – Herv