5

Possible Duplicate:
Fastest algorithm for primality test

Would appreciate a reference to sample code for fast primality testing in C#, preferably using BigInteger or other variable size type.

Community
  • 1
  • 1
Halfdan Faber
  • 223
  • 1
  • 3
  • 4

1 Answers1

17

This is a Miller Rabin test in c#:

    bool MillerRabin(ulong n)
    {
        ulong[] ar;
        if (n < 4759123141) ar = new ulong[] { 2, 7, 61 };
        else if (n < 341550071728321) ar = new ulong[] { 2, 3, 5, 7, 11, 13, 17 };
        else ar = new ulong[] { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
        ulong d = n - 1;
        int s = 0;
        while ((d & 1) == 0) { d >>= 1; s++; }
        int i, j;
        for (i = 0; i < ar.Length; i++)
        {
            ulong a   = Math.Min(n - 2, ar[i]);
            ulong now = pow(a, d, n);
            if (now == 1) continue;
            if (now == n - 1) continue;
            for (j = 1; j < s; j++)
            {
                now = mul(now, now, n);
                if (now == n - 1) break;
            }
            if (j == s) return false;
        }
        return true;
    }

    ulong mul(ulong a, ulong b, ulong mod)
    {
        int i;
        ulong now = 0;
        for (i = 63; i >= 0; i--) if (((a >> i) & 1) == 1) break;
        for (; i >= 0; i--)
        {
            now <<= 1;
            while (now > mod) now -= mod;
            if (((a >> i) & 1) == 1) now += b;
            while (now > mod) now -= mod;
        }
        return now;
    }

    ulong pow(ulong a, ulong p, ulong mod)
    {
        if (p == 0) return 1;
        if (p % 2 == 0) return pow(mul(a, a, mod), p / 2, mod);
        return mul(pow(a, p - 1, mod), a, mod);
    }
Saeed Amiri
  • 22,252
  • 5
  • 45
  • 83
  • Saeed: Thanks much. I converted the code to take the BigInteger type. Works great and extremely fast. – Halfdan Faber Nov 26 '10 at 05:15
  • This is beautiful. I was wondering how you came up with the ``ar`` values? Any references I could look up? – tweaksp Apr 30 '14 at 00:50
  • @Chris, take a look at the sequence [A014233](http://oeis.org/A014233). – Saeed Amiri May 23 '14 at 12:31
  • This is not a primality test, but in fact a compositeness test. Just because a number passes does not mean it is prime. – Dubslow Nov 30 '15 at 22:49
  • 1
    @Dubslow, For the ulong range this is correct primality test, I think you didn't notice the usage of "ar" array. This is not just miller rabin test. If there wasn't any array which excludes some important exceptions where miller rabin can' t find them then you were right. See the sequences of composite numbers that miller rabin fails and compare it by provided algorithm for the range of ulong. You can find the sequence in my other comment above. – Saeed Amiri Dec 01 '15 at 15:43
  • Okay, for a certain upper bound this algorithm and given bases are known to be a primality test. But in general it is not (even if the upper bound is beyond the non-bignum range) – Dubslow Dec 01 '15 at 19:17
  • @Dubslow, Everyone by looking at first few sentences in wiki link provided in my answer, will know that M.R in general is false positive algorithm. – Saeed Amiri Dec 02 '15 at 05:28
  • The multiplication code shown here breaks for 64-bit numbers (i.e. it works only for moduli where bit 63 is not set), and it is excruciatingly slow because it computes `now % mod` via repeated subtraction instead of using operator %. For bigger inputs it can take several minutes for a single call... Using the builtin operator % fixes that problem somewhat (1 microsecond per call for 63-bit inputs) but the problem with bit 64-bit inputs remains and hence the code can only be used for inputs up to 2^63 - 1. That's good enough for many code challenges that go only up to 10^18, though. – DarthGizka Oct 09 '16 at 15:06
  • @SaeedAmiri, When I run the code it is reporting true for some even numbers like 12,14,16. I added "if (n % 2 == 0 || n == 1) return false;" to the top and seemed to correct it. – SunsetQuest Jun 22 '17 at 15:49
  • @Sunsetquest Actually that skips 1 and 2 being reported as prime. So it should be something such as "if (n == 1 || n == 2) return true; if ((n & 1) == 0) return false;" – H. de Jonge Apr 26 '20 at 13:58
  • @DarthGizka, yes it breaks for 64 bit values, if a+b is bigger than 64bits. But you can use big integer for that case (of course it will be slower). But modulo function is extremely fast, and it is a well-known technique to compute huge multiplication modulo like this (not my invention). It doesn't repeatedly subtract, you may have a look at this explanation: https://www.geeksforgeeks.org/multiply-large-integers-under-large-modulo/ If we compute a mod z then b mod z and multiply the two, it may overflow (even for 60 bits) and we get wrong result. So, the above code resolves that issue. – Saeed Amiri Oct 23 '22 at 10:42