欢迎。我正在尝试实施 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
}
对于小大整数来说它工作得很好。但如果我的程序需要计算包含超过 16 位的数字，程序就会冻结。例如，在成功检查数字是否为素数后，程序突然没有响应。我不明白这怎么可能。如果它检查了一个大数字，那么再次检查另一个数字应该没有问题。即使调试器也没有帮助，因为
step options编辑。更改 BigInteger.ModPow 的模函数有帮助。不幸的是，现在对于更大的数字，超过 3000 位，它永远不会返回质数，这是相当不可能的。或者说真正重要的数字很难找到？
嗯，在我的工作站（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);
这是我的代码，您可以在其中检查从0到小数的素数。MaxValue=79228162514264337593543950335
更新
我做了一些调整以使程序更快
在：
英特尔(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.小数.MaxValue = 79228162514264337593543950335
小数以下的最大数字。测试的MaxValue是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();
}
}
}
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
-----------------------------------------