MillerRabin primality test in C#

Question

Welcome. I am trying to implement MillerRabin test for checking if large given number is a prime. Here is my code:

 public static bool MillerRabinTest(BigInteger number)
        {

            BigInteger d;
            var n = number - 1;
            var s = FindK(n, out d);

            BigInteger a = 2;
            BigInteger y = Calc(a, d, number);  //a^d mod number
            if (y != BigInteger.One && y != n)
            {
                for (var r = 1; r <= s - 1; r++)
                {
                    y = Calc(y, 2, number);
                    if (y == 1)
                        return false;  
                }

                if (y != n)
                    return false;
            }
            return true; //it is probably prime
        }

It is working fine for small Bigintegers. But if my programs needs to evalute numbers containing of more than 16 bits, program freezes. For instance after succesful checking if number is a prime, program suddenly is not responsive. I dont understand how is that possible. If it checked one big number, it should have no problem for checking another one again. Even debugger is not being helpful ,becasue step options disappear. I can share more code of functions if needed. Above function is working correctly for small numbers.

EDIT. Changing my modulo function for BigInteger.ModPow helped. Unfortunately now for bigger numbers, more than 3000 bits it is never returning prime number which is rather impossible. Or really prme numbers are hard to find out?

Answer 1

Well, it takes about 5 seconds at my workstation (Core i5 3.2GHz, IA64 .Net 4.5) to test for being prime for numbers equals to 2**3000:

  public static class PrimeExtensions {
    // Random generator (thread safe)
    private static ThreadLocal<Random> s_Gen = new ThreadLocal<Random>(
      () => {
        return new Random();
      }
    );

    // Random generator (thread safe)
    private static Random Gen {
      get {
        return s_Gen.Value;
      }
    }

    public static Boolean IsProbablyPrime(this BigInteger value, int witnesses = 10) {
      if (value <= 1)
        return false;

      if (witnesses <= 0)
        witnesses = 10;

      BigInteger d = value - 1;
      int s = 0;

      while (d % 2 == 0) {
        d /= 2;
        s += 1;
      }

      Byte[] bytes = new Byte[value.ToByteArray().LongLength];
      BigInteger a;

      for (int i = 0; i < witnesses; i++) {
        do {
          Gen.NextBytes(bytes);

          a = new BigInteger(bytes);
        }
        while (a < 2 || a >= value - 2);

        BigInteger x = BigInteger.ModPow(a, d, value);
        if (x == 1 || x == value - 1)
          continue;

        for (int r = 1; r < s; r++) {
          x = BigInteger.ModPow(x, 2, value);

          if (x == 1)
            return false;
          if (x == value - 1)
            break;
        }

        if (x != value - 1)
          return false;
      }

      return true;
    }
  }

Test and benchmark

  BigInteger value = BigInteger.Pow(2, 3217) - 1; // Mersenne prime number (2.5e968)

  Stopwatch sw = new Stopwatch();

  sw.Start();

  Boolean isPrime = value.IsProbablyPrime(10);

  sw.Stop();

  Console.Write(isPrime ? "probably prime" : "not prime");
  Console.WriteLine();
  Console.Write(sw.ElapsedMilliseconds);