5

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)?

Niklas B.
  • 92,950
  • 18
  • 194
  • 224
user3717225
  • 99
  • 1
  • 5

5 Answers5

14

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 264, 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}.

Niklas B.
  • 92,950
  • 18
  • 194
  • 224
  • What's with the template? The question is tagged **c**. – Cristian Ciupitu Jun 07 '14 at 11:32
  • Reference / Reasoning? – John Dvorak Jun 07 '14 at 11:33
  • 2
    @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. – CodesInChaos Jun 07 '14 at 11:33
  • @JanDvorak Click on the first link of my answer? – Niklas B. Jun 07 '14 at 11:33
  • 1
    @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`. – Niklas B. Jun 07 '14 at 11:34
  • 1
    Also the implementations of `MultiplyMod` and `PowerMod` can be looked up in the complete source file which I linked – Niklas B. Jun 07 '14 at 12:37
  • @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 ;) – Niklas B. Mar 18 '16 at 11:51
  • @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 – Niklas B. Mar 18 '16 at 15:34
  • @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! – Brett Hale Mar 18 '16 at 15:38
  • BTW - you don't need (9) base tests, [Sinclair](http://miller-rabin.appspot.com/) 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:) – Brett Hale Mar 18 '16 at 16:21
  • @BrettHale sure, but that's not the first 7 primes ;) They are easier to type without making mistakes :P – Niklas B. Mar 18 '16 at 16:52
  • I was looking at some older [questions](http://stackoverflow.com/questions/9620181/a-way-to-find-the-nearest-prime-number-to-an-unsigned-long-integer-32-bits-wid/9620391#9620391), 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. – Brett Hale Jun 12 '16 at 11:09
  • Accuracy formula is at https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Accuracy_2 – Sam Ginrich Jul 03 '22 at 13:44
6

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.

user448810
  • 17,381
  • 4
  • 34
  • 59
  • 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? – Todd Lehman Jul 30 '14 at 05:25
  • It's 325, not 235. And you test all seven bases. – user448810 Jul 30 '14 at 11:14
  • 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*. – Todd Lehman Jul 30 '14 at 19:14
  • 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. – Brett Hale Mar 18 '16 at 15:56
  • See my [blog](https://programmingpraxis.files.wordpress.com/2012/09/primenumbers.pdf). – user448810 Jul 16 '16 at 00:54
1

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.

CodesInChaos
  • 106,488
  • 23
  • 218
  • 262
1

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!

Brett Hale
  • 21,653
  • 2
  • 61
  • 90
  • 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). – DarthGizka Mar 06 '16 at 10:17
  • @DarthGizka - 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. – Brett Hale Mar 06 '16 at 10:56
  • 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. – DarthGizka Mar 06 '16 at 21:21
  • @DarthGizka - 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:) – Brett Hale Mar 07 '16 at 01:32
-2

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.

rossum
  • 15,344
  • 1
  • 24
  • 38
  • 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 – Brett Hale Mar 18 '16 at 16:11