MillerRabin primality test in C#

问题:

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?

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

回答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();
        }
    }
}