如何在 C# 中生成数百万个素数
使用 C# 生成数百万个素数的代码。
引言
这是我希望成为一系列文章的开篇之作,旨在描述我的一些关于密码学和素数的编码实践。你可以在互联网上(某个地方)找到这些内容,但我试图为新手提供一个起点。在密码学和素数测试中,你需要做的第一件事就是获得一些素数。一千个是一个不错的开始,但一百万个或更多更好。而且一旦解决了 bug,我更信任电脑和代码,而不是我自己手动输入一百万个素数的能力。这可是个大工程。另一种选择是购买这样的列表,但我比较穷。而且也有免费的列表,但它们存在格式问题。(而且,文件中可能存在错误。虽然不太可能,但也有可能。)
最初,我以为我可以从 7 开始,取它的平方根,然后检查从 2 到该平方根的所有除数。如果没有找到除数,这个数就是素数,将其放入一个 List
。如果找到除数,就跳到下一个数,然后重复。对于一个不太大的单个数字,或者(我想)一千个素数的列表来说,这还可以,但以这种方式检查大量数字会变得很慢。
深夜,我在一个密码学论坛上读到一个非常有才的想法。设定一个你想达到的最高数字的上限,然后从最小的素数开始,将其平方数标记为非素数,并以该素数的两倍递增 march 到最高数字。当这个最小的素数完成后,转到下一个未标记的素数,重复平方和标记循环。(如果这听起来很复杂,看看下面的图片。图片是我理解复杂问题的方式。)
至于你运行速度能得到什么... 这将取决于你的电脑(CPU 速度和内存大小及速度,同时运行了多少程序等)。如果代码运行太慢,就降低 iMax 的数字。
背景
避免测试和标记所有能被 2、3 和 5 整除的数字。这可以通过取每个数字的 30 取模来实现。(这会将需要检查的数字数量从 30 个减少到 8 个。) Rick Oden 的文章(这里)中有一个关于这个想法的不错图示。想象一个 30 列的表格,移除所有偶数列以及能被 3 或 5 整除的列。前 6 行是这样的:
这些数字中的每一个都有可能是素数。为了论证起见,跳过数字“1”。我只是把它放进表格,因为否则会出现一个空位,顶行包含了它下面的每个单元格的“N % 30”,看起来更美观。(一些聪明人会反对:但“1”不是素数!我说,放轻松。它是一个特殊情况。“1”是乘法的单位元。现在先不管它。)
所以从数字“7”开始,将其平方得到“49”。我们将其标记为非素数。我们还以 7(这是一个偶数,不在我们的表格中)为增量前进到“56”。跳过它,再前进 7 到达“77”。为了节省时间,我们不会尝试以 7 为增量。(奇数加上奇数总是偶数(不在我们的表格中)。奇数加上该素数的两倍将是奇数,并且可能在表格中。)我们加上 14 来到达下一个有因子 7 的数字,即“91”……就这样,以 14 的间隔前进,并将该数字标记为非素数。如果一个数字不在我们的表格中,我们仍然继续前进,直到达到最大数字。我们的表格看起来会是这样的:
表格中的下一个未标记的数字是“11”。它的平方是“121”。将其标记为非素数。加上 22 到达“143”,标记为非素数。依此类推,以 22 为步长……
如果我们前进到一个已标记的单元格,它就是合数,可以跳过。(“49”是我们的第一个这样的单元格。我们跳过它。下一个单元格是 53,将其平方并标记,然后以 106 为增量前进,就又回到正常操作了。)
要编写此代码,请分配一个大的字节数组来工作。Bytes_Allocated * 30 = N = 要 march 到的最大素数。数组中的每个字节对应 N / 30。字节中的每个位对应 N % 30 的 8 个可能值。(1、7、11、13、17、19、23、29)对应 8 个可能值。每个位对应一个值。所以,如果我们分配 1000 字节,我们可以覆盖直到 30000 的素数。代码完成后,任何未设置的位都将代表一个素数。在这种情况下,我将找到的每个素数添加到 List
中。
Using the Code
以下 C# 代码可以直接复制粘贴到 Visual Studio 中。将其集成到你自己的按钮点击事件中。之后,一切照旧——自行承担使用代码的风险,我不对你可能沉迷于追逐素数或素数测试负责。
这里没有花哨的类。只是为了速度而使用了一些全局变量、按引用传递,以及上述优化,以确定哪些才需要被考虑为素数。
List<UInt32> listPrimes = new List<UInt32>();
private void button8_Click(object sender, EventArgs e)
{
UInt32 iBump = 2;
UInt32 iPrime = 0;
UInt32 iPrime2 = 0;
UInt32 iPrimeSqr = 0;
UInt32 iMax = 20000; // 0 sec
iMax = 200000; // 0.2 sec
iMax = 2000000; // 0.8 sec
iMax = 20000000; // 8.3 sec; gens List of first 8.8M primes
//iMax = 200000000; // 103 sec; gens List of first 70M primes
UInt32 iMax8 = iMax * 8;
UInt64 iiMax8 = iMax * 8;
if ( iiMax8 > UInt32.MaxValue )
{
System.Diagnostics.Debug.WriteLine (
"Too many numbers requested for this data type!");
return;
}
// Start out all = false
byte[] aPrimesBool = new byte [ iMax ];
UInt32 iFactor = 0;
DateTime dt0 = DateTime.Now;
DateTime dt1 = DateTime.Now;
// Clear the old primes list
listPrimes.Clear();
listPrimes.Add ( 2 );
listPrimes.Add ( 3 );
listPrimes.Add ( 5 );
// Every odd number is potentially prime; set every byte (8 bits at a time)
for ( iFactor = 0; iFactor < iMax; iFactor += 1 )
{
aPrimesBool [ iFactor ] = 0xff;
}
dt1 = DateTime.Now;
System.Diagnostics.Debug.WriteLine ( "start after : " +
( dt1 - dt0).TotalSeconds + " sec");
// Init _iGR and i_GM to 7
_iGW = 0; // == 7 / 30
_iGR = 7; // == 7 % 30
for ( iPrime = 7; iPrime < iMax8; iPrime += iBump )
{
if ( iPrime % 1000000 == 1 )
{
dt1 = DateTime.Now;
System.Diagnostics.Debug.WriteLine ( "reached: " +
iPrime + " after : " + ( dt1 - dt0).TotalSeconds + " sec");
}
if ( GetCodePrime ( ref aPrimesBool ) ) // iPrime passed via _iGR and _iGW
{
listPrimes.Add ( iPrime );
iPrime2 = iPrime * 2;
iPrimeSqr = iPrime * iPrime;
// Every multiple of this prime is not prime, loop and mark them as not prime
for ( iFactor = iPrimeSqr; iFactor < iMax8; iFactor += iPrime2 )
{
ClrCodePrime ( ref iFactor, ref aPrimesBool );
}
}
// Alternate upping to the next possible prime by 2 and 4
// ( 7 + 4, 11 + 2, 13 + 4, 17 + 2, 19 + 4,
// 23 + 2, 25 + 4, 29 + 2, 31 + 4, 35 + 2, ... )
if ( iBump == 2 )
iBump = 4;
else
iBump = 2;
// Counters for GetCodePrime()
if ( _iGR < 29 )
{
_iGR += iBump; // up the remainder
}
else
{
_iGR = 1;
_iGW++;
}
}
dt1 = DateTime.Now;
System.Diagnostics.Debug.WriteLine ( "done: " + iPrime );
System.Diagnostics.Debug.WriteLine ( "count = " + listPrimes.Count );
// wipe out the byte array
aPrimesBool = null;
System.GC.WaitForPendingFinalizers();
}
UInt32 _iCR = 0; // Clear Remainder
UInt32 _iCW = 0; // Clear Whole
private void ClrCodePrime ( ref UInt32 iNum, ref byte[] bA )
{
_iCR = iNum % 30;
switch ( _iCR )
{
case 1:
_iCW = iNum / 30; // Only calculate
// the whole number if the remainder is worth it
bA [ _iCW ] &= 0xfe; // Clear bit 2^0 == 1
return;
case 7:
_iCW = iNum / 30;
bA [ _iCW ] &= 0xfd; // Clear bit 2^1 == 2
return;
case 11:
_iCW = iNum / 30;
bA [ _iCW ] &= 0xfb; // Clear bit 2^2 == 4
return;
case 13:
_iCW = iNum / 30;
bA [ _iCW ] &= 0xf7; // Clear bit 2^3 == 8
return;
case 17:
_iCW = iNum / 30;
bA [ _iCW ] &= 0xef; // Clear bit 2^4 == 16
return;
case 19:
_iCW = iNum / 30;
bA [ _iCW ] &= 0xdf; // Clear bit 2^5 == 32
return;
case 23:
_iCW = iNum / 30;
bA [ _iCW ] &= 0xbf; // Clear bit 2^6 == 64
return;
case 29:
_iCW = iNum / 30;
bA [ _iCW ] &= 0x7f; // Clear bit 2^7 == 128
return;
}
// if remainder is anything else, ignore the number because
// it is divisible by 2, 3, or 5.
return;
}
UInt32 _iGR = 0; // Get remainder
UInt32 _iGW = 0; // Get whole
private bool GetCodePrime ( ref byte[] bA )
{
//
// Returns true if ( _iGR + 30 * _iGW ) is prime
//
switch ( _iGR )
{
case 1: // 1 + 30 * _iGW
return ( bA[_iGW] & 0x01 ) == 0x01;
case 7: // 7 + 30 * _iGW
return ( bA[_iGW] & 0x02 ) == 0x02;
case 11: // 11 + 30 * _iGW
return ( bA[_iGW] & 0x04 ) == 0x04;
case 13: // 13 + 30 * _iGW
return ( bA[_iGW] & 0x08 ) == 0x08;
case 17: // 17 + 30 * _iGW
return ( bA[_iGW] & 0x10 ) == 0x10;
case 19: // 19 + 30 * _iGW
return ( bA[_iGW] & 0x20 ) == 0x20;
case 23: // 23 + 30 * _iGW
return ( bA[_iGW] & 0x40 ) == 0x40;
case 29: // 29 + 30 * _iGW
return ( bA[_iGW] & 0x80 ) == 0x80;
}
return false;
}
为了更容易理解,我避免了使用多线程。如果你找到了让这段代码运行更快的方法,请务必分享你做了哪些修改。
关注点
在 `ClrCodePrime` 例程中进行了一些巧妙的位清除操作。而且,因为它不是良好的类结构,在更改东西时要小心。非常有可能不小心弄坏某些东西。(这部分解释了为什么类比这种“恐龙风格”的代码更好。我坚持认为老式代码运行速度更快。但如果你有更快的类实现,请分享,我们可以一起决定它是否更好。)
在生成大量素数时,单独计算每个素数(并检查它是否有因子)的效率较低,而让每个素数设置它上面所有可因数化的合数到某个请求的最大值的效率更高。
历史
- 初始文章发布日期 = 2011/9/16。