The Ultimate Primality Test

Figuring out whether or not a number is prime has been a long quest by mathematicians and computer scientists, especially in the last 50 years with the development of computer cryptography, in particular the RSA algorithm. The common, centuries-old method for checking whether a number N is prime is very straightforward: test all the numbers "i" between 2 and SQRT(N), and if N is divisible by "i" then N is not prime, otherwise it is prime. There is a number of optimizations to this algorithm, but its core remains the same, and so does the lower bound execution time of O(SQRT(N)). Which works well for N = 32212254719, since SQRT(N) ~= 179478. However - try this:
N = 225251798594466661409915431774713195745814267044878909733007331390393510002687. This is a special prime number, one of the Woodall prime numbers. Out of curiosity, its root square is 4.74606994E38, which becomes intractable.
What are the options? The best primality test algorithm was actually discovered in 2003 by IIT researchers (I've always been a big fan of IIT as an institution and its research team). The breakthrough paper entitled "Primes is in P" sent shockwaves thru the theoretical mathematical community since it proclaimed to have found a solution to the primality test that runs in O(Log(N)), which for the above case would mean running in a mere 78 steps. However, its implementation is really, really convoluted to a point that I haven't seen any readable implementation of the algorithm anywhere - to be very frank, I haven't found not even an unreadable one.
It turns out that for many decades there has been another method that has been vastly overseen by the community: the Miller-Rabin Primality Test. Miller-Rabin is not very complicated to understand once you get the concept of how to implement the Witness function. After that, all it is needed is to check the witness for few known primes (all primes from 2 to 500), and the algorithm is over. It runs extremely fast, as fast as the "Primes is in P" algorithm, or even faster. The only reason it has been overlooked is because it relies on an unproven hypothesis, the Generalized Riemann hypothesis. There is very little doubt among the mathematical community that the Riemann hypothesis is true, but since it hasn't been proved yet, the Miller-Rabin algorithm has not been given the respect it deserves. Which is a shame - it is easy to implement, fast as light, and although there are claims that there might be cases where the algorithm does not work (because of the Riemann hypothesis thing), no one has ever produced a prime number that the algorithm failed to find. Hence, next time that you need to check whether a number, a big one, is prime, feel free to use this implementation. Btw, take a look at how fast it runs:

225251798594466661409915431774713195745814267044878909733007331390393510002687 is prime (time: 0.0029934 seconds)

Code is here (C#, using System.Numerics for BigInteger class):

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;

namespace TheUltimatePrimalityTest
{
    class Program
    {
        static void Main(string[] args)
        {
            if (args.Length != 1)
            {
                Console.WriteLine("Usage: TheUltimatePrimalityTest.exe <number>");
                return;
            }

            long ticksBefore = DateTime.Now.Ticks;
            BigInteger number = StringToBigInteger(args[0]);
            long ticksAfter = DateTime.Now.Ticks;
            double latency = (ticksAfter - ticksBefore) / 10000000.0;

            if (IsPrimeMillerRabin(number))
            {
                Console.WriteLine("{0} is prime (time: {1} seconds)", number, latency);
            }
            else
            {
                Console.WriteLine("{0} is not prime (time: {1} seconds)", number, latency);
            }
        }

        public static BigInteger StringToBigInteger(string str)
        {
            if (str == null || str == "")
            {
                return 0;
            }
            else
            {
                BigInteger n = 0;
                foreach (char c in str)
                {
                    n = 10 * n + (int)(c - '0');
                }
                return n;
            }
        }

        public static bool IsPrimeMillerRabin(BigInteger n)
        {
            //It does not work well for smaller numbers, hence this check
            int SMALL_NUMBER = 1000;

            if (n <= SMALL_NUMBER)
            {
                return IsPrime(n);
            }

            int MAX_WITNESS = 500;
            for (long i = 2; i <= MAX_WITNESS; i++)
            {
                if (IsPrime(i) && Witness(i, n) == 1)
                {
                    return false;
                }
            }

            return true;
        }

        public static BigInteger SqRtN(BigInteger N)
        {
            /*++
             *  Using Newton Raphson method we calculate the
             *  square root (N/g + g)/2
             */
            BigInteger rootN = N;
            int count = 0;
            int bitLength = 1;
            while (rootN / 2 != 0)
            {
                rootN /= 2;
                bitLength++;
            }
            bitLength = (bitLength + 1) / 2;
            rootN = N >> bitLength;

            BigInteger lastRoot = BigInteger.Zero;
            do
            {
                if (lastRoot > rootN)
                {
                    if (count++ > 1000)                   // Work around for the bug where it gets into an infinite loop
                    {
                        return rootN;
                    }
                }
                lastRoot = rootN;
                rootN = (BigInteger.Divide(N, rootN) + rootN) >> 1;
            }
            while (!((rootN ^ lastRoot).ToString() == "0"));
            return rootN;
        }
       
        public static bool IsPrime(BigInteger n)
        {
            if (n <= 1)
            {
                return false;
            }

            if (n == 2)
            {
                return true;
            }

            if (n % 2 == 0)
            {
                return false;
            }

            for (int i = 3; i <= SqRtN(n) + 1; i += 2)
            {
                if (n % i == 0)
                {
                    return false;
                }
            }
            return true;
        }

        private static int Witness(long a, BigInteger n)
        {
            BigInteger t, u;
            BigInteger prev, curr = 0;
            BigInteger i;
            BigInteger lln = n;

            u = n / 2;
            t = 1;
            while (u % 2 == 0)
            {
                u /= 2;
                t++;
            }

            prev = BigInteger.ModPow(a, u, n);
            for (i = 1; i <= t; i++)
            {
                curr = BigInteger.ModPow(prev, 2, lln);
                if ((curr == 1) && (prev != 1) && (prev != lln - 1)) return 1;
                prev = curr;
            }
            if (curr != 1) return 1;
            return 0;
        }
    }
}