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

using System;
using System.Numerics;
using System.Threading;
using System.Diagnostics;

class Program
{
    static void Main()
    {
        TestPrime(BigInteger.Pow(2, 17) - 1, 10); // Small Mersenne prime (2^17 - 1)
        TestPrime(new BigInteger(1009), 10); // Small prime number
        TestPrime(new BigInteger(1010), 10); // Small non-prime number
        TestPrime(BigInteger.Pow(2, 1279) - 1, 10); // Large Mersenne prime (2^1279 - 1)
        TestPrime(BigInteger.Pow(2, 3217) - 1, 10); // Very large Mersenne prime (2^3217 - 1)
        TestPrime(new BigInteger(long.MaxValue) * 2, 10); // Large non-prime number
    }

    static void TestPrime(BigInteger value, int certainty)
    {
        Stopwatch sw = new Stopwatch();
        sw.Start();
        bool isPrime = value.IsProbablyPrime(certainty);
        sw.Stop();
        Console.WriteLine($"Testing: {value}");
        Console.WriteLine(isPrime ? "Probably prime" : "Not prime");
        Console.WriteLine($"Elapsed time: {sw.ElapsedMilliseconds} ms");
        Console.WriteLine("-----------------------------------------");
    }
}




public static class BigIntegerExtensions
{
    private static Random random = new Random();

    public static bool IsProbablyPrime(this BigInteger source, int certainty)
    {
        if (source == 2 || source == 3)
            return true;
        if (source < 2 || source % 2 == 0)
            return false;

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

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

        for (int i = 0; i < certainty; i++)
        {
            BigInteger a = RandomBigInteger(2, source - 2);
            BigInteger x = BigInteger.ModPow(a, d, source);
            if (x == 1 || x == source - 1)
                continue;

            for (int r = 1; r < s; r++)
            {
                x = BigInteger.ModPow(x, 2, source);
                if (x == 1)
                    return false;
                if (x == source - 1)
                    break;
            }

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

        return true;
    }

    private static BigInteger RandomBigInteger(BigInteger minValue, BigInteger maxValue)
    {
        if (minValue > maxValue)
            throw new ArgumentException("minValue must be less than or equal to maxValue");

        BigInteger range = maxValue - minValue + 1;
        int length = range.ToByteArray().Length;
        byte[] buffer = new byte[length];

        BigInteger result;
        do
        {
            random.NextBytes(buffer);
            buffer[buffer.Length - 1] &= 0x7F; // Ensure non-negative
            result = new BigInteger(buffer);
        } while (result < minValue || result >= maxValue);

        return result;
    }
}

Testing: 131071
Probably prime
Elapsed time: 1 ms
-----------------------------------------
Testing: 1009
Probably prime
Elapsed time: 0 ms
-----------------------------------------
Testing: 1010
Not prime
Elapsed time: 0 ms
-----------------------------------------
Testing: 10407932194664399081925240327364085538615262247266704805319112350403608059673360298012239441732324184842421613954281007791383566248323464908139906605677320762924129509389220345773183349661583550472959420547689811211693677147548478866962501384438260291732348885311160828538416585028255604666224831890918801847068222203140521026698435488732958028878050869736186900714720710555703168729087
Probably prime
Elapsed time: 183 ms
-----------------------------------------
Testing: 259117086013202627776246767922441530941818887553125427303974923161874019266586362086201209516800483406550695241733194177441689509238807017410377709597512042313066624082916353517952311186154862265604547691127595848775610568757931191017711408826252153849035830401185072116424747461823031471398340229288074545677907941037288235820705892351068433882986888616658650280927692080339605869308790500409503709875902119018371991620994002568935113136548829739112656797303241986517250116412703509705427773477972349821676443446668383119322540099648994051790241624056519054483690809616061625743042361721863339415852426431208737266591962061753535748892894599629195183082621860853400937932839420261866586142503251450773096274235376822938649407127700846077124211823080804139298087057504713825264571448379371125032081826126566649084251699453951887789613650248405739378594599444335231188280123660406262468609212150349937584782292237144339628858485938215738821232393687046160677362909315071
Probably prime
Elapsed time: 1894 ms
-----------------------------------------
Testing: 18446744073709551614
Not prime
Elapsed time: 0 ms
-----------------------------------------

Try it online!

Answer 2

here is my code where you can check the primality from 0 to decimal.MaxValue=79228162514264337593543950335

UPDATE

i made some adjustments to make the program faster

in a:
Intel(R) Atom(TM) @ 1.60GHz
2.00GB RAM
32-bit Operating System

results:
1. UInt32.MaxValue = 4294967295
largest prime number below UInt32.MaxValue is 4294967291
elapsed time is 0.015600 second
2. ulong.MaxValue = UInt64.MaxValue = 18446744073709551615
largest prime number below ulong.MaxValue is 18446744073709551533
elapsed time is 3 minutes and 57.6059176 seconds
3. decimal.MaxValue = 79228162514264337593543950335
largest number below decimal.MaxValue tested is 79228162514264337593543950319 but not known if 79228162514264337593543950319 is a prime number because i interrupted the running of the program after elapsed time is 3 hours and 40 minutes (need to be tested with high specifications laptops)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace PrimalityTest
{
    class Program
    {
        static void Main(string[] args)
        {
            Console.Write("Enter a number: ");
            decimal decimal_number = Convert.ToDecimal(Console.ReadLine());
            DateTime date = DateTime.Now;
            ulong ulong_a;
            ulong ulong_b;
            if (decimal_number <= ulong.MaxValue)
            {
                ulong ulong_number = Convert.ToUInt64(decimal_number);
                if (ulong_number < 2)
                {
                    Console.WriteLine(ulong_number + " is not a prime number");
                }
                else if (ulong_number == 2 || ulong_number == 3)
                {
                    Console.WriteLine(ulong_number + " is a prime number");
                }
                else if (ulong_number % 2 == 0)
                {
                    Console.WriteLine(ulong_number + " is not a prime number and is divisible by 2");
                }
                else
                {
                    ulong_a = Convert.ToUInt64(Math.Ceiling(Math.Sqrt(ulong_number)));
                    for (ulong_b = 3; ulong_b <= ulong_a; ulong_b += 2)
                    {
                        if (ulong_number % ulong_b == 0)
                        {
                            Console.WriteLine(ulong_number + " is not a prime number and is divisible by " + ulong_b);
                            goto terminate_ulong_primality_test;
                        }
                    }
                    Console.WriteLine(ulong_number + " is a prime number");
                }
                terminate_ulong_primality_test:
                {
                }
            }
            else
            {
                if (decimal_number % 2 == 0)
                {
                    Console.WriteLine(decimal_number + " is not a prime number and is divisible by 2");
                }
                else
                {
                    ulong_a = Convert.ToUInt64(Math.Ceiling(Math.Sqrt(ulong.MaxValue) * Math.Sqrt(Convert.ToDouble(decimal_number / ulong.MaxValue))));
                    for (ulong_b = 3; ulong_b <= ulong_a; ulong_b += 2)
                    {
                        if (decimal_number % ulong_b == 0)
                        {
                            Console.WriteLine(decimal_number + " is not a prime number and is divisible by " + ulong_b);
                            goto terminate_decimal_primality_test;
                        }
                    }
                    Console.WriteLine(decimal_number + " is a prime number");
                }
                terminate_decimal_primality_test:
                {
                }
            }
            Console.WriteLine("elapsed time: " + (DateTime.Now - date));
            Console.ReadKey();
        }
    }
}

Answer 3

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