C# 中的 Miller Rabin 素数测试

Question

欢迎。我正在尝试实施 MillerRabin 测试来检查给定的大数是否是素数。这是我的代码：

 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
        }

它适用于小型 Biginteger。但是如果我的程序需要计算超过 16 位的数字，程序就会冻结。例如，在成功检查数字是否为素数后，程序突然没有响应。我不明白这怎么可能。如果它检查了一个大数字，那么再次检查另一个应该没有问题。甚至调试器也没有帮助，因为 step options 消失了。如果需要，我可以分享更多功能代码。上面的函数对小数字正常工作。

编辑。更改 BigInteger.ModPow 的模函数有帮助。不幸的是，现在对于更大的数字，超过 3000 位，它永远不会返回质数，这是不可能的。还是真的 prme 数字很难找到？

Answer 1

这是我的代码，您可以在其中检查从 0 到小数的素数。MaxValue=79228162514264337593543950335

更新

我做了一些调整以使程序更快

在：
Intel(R) Atom(TM) @ 1.60GHz
2.00GB 内存
32位操作系统

结果：
1. UInt32.MaxValue = 4294967295
UInt32.MaxValue 以下的最大素数是 4294967291
经过的时间是 0.015600 秒
2. ulong.MaxValue = UInt64.MaxValue = 18446744073709551615
ulong.MaxValue 以下的最大素数是 18446744073709551533
经过的时间是 3 分 57.6059176 秒
3. decimal.MaxValue = 79228162514264337593543950335
小数点以下的最大数字。测试的最大值为 79228162514264337593543950319 但不知道 79228162514264337593543950319 是否为质数，因为我在经过 3 小时 40 分钟后中断了程序的运行（需要用高规格的笔记本电脑测试）

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 2

嗯，我的工作站（Core i5 3.2GHz，IA64 .Net 4.5）大约需要 5 秒 来测试等于 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;
    }
  }

测试和基准测试

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