整数分解：令人畏惧的素数列表
使用大型素数列表和试除法算法，
下载
目录
背景
本文是一系列文章的一部分。每篇文章都着重于整数分解的一个方面。
- 整数分解：试除法算法[^]：关注试除法算法的变体及其各自的效率。
- 整数分解：令人头疼的素数列表[^]：关注通过压缩处理大量素数列表的方法。
- 整数分解：优化小因子检查[^]：关注比除法更快地检查小因子的方法。
- 整数分解：逆转乘法[^]：关注优化试除法算法的另一种方法。
引言
试除法算法是通过检查所有可能的因子来分解一个整数，而所有有趣的因子都是素数。问题是如何处理大量的素数列表。
第一种方法：简单的素数列表
当使用素数列表时，第一种方法是使用一个包含素数的简单整数数组。
在 64K 以下，一个素数可以存储在 int16（2 字节）中，超过这个范围，一个素数则更适合存储在 int32（4 字节）中。
简单素数列表的数字
| 上限 | 素数数量 | 数组大小 |
| 65,536 | 6,542 | 13,084 |
| 100,000 | 9,592 | 38,368 |
| 500,000 | 41,538 | 166,152 |
| 1,000,000 | 78,498 | 313,992 |
| 5,000,000 | 348,513 | 1,394,052 |
| 10,000,000 | 664,579 | 2,658,316 |
| 100,000,000 | 5,761,455 | 23,045,820, |
| 2^32-1 4,294,967,295 | 203,280,221 | 813,120,884 |
正如你所见，当列表非常庞大时，情况会变得非常糟糕。
示例代码
// Trial Division: Square Root + Prime list + Wheel
// Check all numbers until Square Root, start with a list of primes
// and Skip useless factors with a wheel
long long TD_SRPW(long long Cand) {
long long SPrimes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191,
193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347,
349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421,
431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
599, 601, 0 };
long long Wheel[] = { 6, 4, 2, 4, 2, 4, 6, 2, 0 };
long long Div;
Count = 0;
long long Top = sqrt(Cand);
// Check small primes
for (long long Ind = 0; SPrimes[Ind] != 0; Ind++) {
Div = SPrimes[Ind]; // for debug purpose
if (SPrimes[Ind] > Top) {
return Cand;
}
Count++;
if (Cand % SPrimes[Ind] == 0) {
return SPrimes[Ind];
}
}
// Start the Wheel
Div = 601;
while (Div <= Top) {
for (long long Ind = 0; Wheel[Ind] != 0; Ind++) {
if (Div > Top) {
return Cand;
}
Count++;
if (Cand % Div == 0) {
return Div;
}
Div += Wheel[Ind];
}
}
return Cand;
}
TrialDivLP.cpp 在 TrialDivLP.zip 中.
第二种方法：通过位字段数组编码压缩素数列表
与其使用素数列表，不如采用编码一个数字是否为素数的方法，这是布尔信息，可以存储在一个位中。
采用这种方法，可以应用经典试除法算法变体中的相同优化，并且节省巨大。
所有数字的位字段
非常简单的想法，每个数字 1 位，或者一个字节编码 8 个数字。
数字如下
| 上限 | 素数数量 | 数组大小 |
| 65,536 | 6,542 | 8,192 |
| 100,000 | 9,592 | 12,500 |
| 500,000 | 41,538 | 62,500 |
| 1,000,000 | 78,498 | 125,000 |
| 5,000,000 | 348,513 | 625,000 |
| 10,000,000 | 664,579 | 1,250,000 |
| 100,000,000 | 5,761,455 | 12,500,000, |
| 2^32-1 4,294,967,295 | 203,280,221 | 536,870,912 |
所有奇数的位字段
稍微优化一下，每个奇数 1 位，或者一个字节编码 16 个数字（8 个奇数）。
数字如下所示
| 上限 | 素数数量 | 数组大小 |
| 65,536 | 6,542 | 4,096 |
| 100,000 | 9,592 | 6,250 |
| 500,000 | 41,538 | 31,250 |
| 1,000,000 | 78,498 | 62,500 |
| 5,000,000 | 348,513 | 312,500 |
| 10,000,000 | 664,579 | ,625,000 |
| 100,000,000 | 5,761,455 | 6,250,000, |
| 2^32-1 4,294,967,295 | 203,280,221 | 268,435,456 |
使用“30轮”筛选的所有数字的位字段
然后更加优化，使用“轮”（wheel），每 30 个数字只有 8 个可能的素数，或者一个字节编码 30 个数字（8 个可能的素数）。
数字是
| 上限 | 素数数量 | 数组大小 |
| 65,536 | 6,542 | 2,188 |
| 100,000 | 9,592 | 3,333 |
| 500,000 | 41,538 | 16,667 |
| 1,000,000 | 78,498 | 33,333 |
| 5,000,000 | 348,513 | 166,667 |
| 10,000,000 | 664,579 | 333,334 |
| 100,000,000 | 5,761,455 | 3,333,334, |
| 2^32-1 4,294,967,295 | 203,280,221 | 143,165,577 |
利用素数列表的代码（存储在头文件中）。
// Trial Division: Square Root + Long Prime list + Wheel
// Check all numbers until Square Root, start with a list of primes
// and Skip useless factors with a wheel
long long TD_SRLPW(long long Cand) {
long long WPrimes[] = { 2, 3, 5, 0 };
long long Wheel[] = { 6, 4, 2, 4, 2, 4, 6, 2, 0 };
// unsigned char EPrimes[] = { 0xfe, 0xdf, 0xef, 0x7e, 0xb6, 0xdb, 0x3d, ...
// Encoded list moved to header file
long long Div;
Count = 0;
long long Top = sqrt(Cand);
// Check small primes
for (long long Ind = 0; WPrimes[Ind] != 0; Ind++) {
Div = WPrimes[Ind]; // for debug purpose
if (WPrimes[Ind] > Top) {
return Cand;
}
Count++;
if (Cand % WPrimes[Ind] == 0) {
return WPrimes[Ind];
}
}
// start the wheel with encoded primes
Div = 1;
int EInd = 0;
while (Div <= Top) {
unsigned char Flag = EPrimes[EInd];
unsigned char Mask = 1;
// if (Flag == 0) {
if (EInd >= EPrimesLen) {
break;
}
for (long long Ind = 0; Wheel[Ind] != 0; Ind++) {
if (Div > Top) {
break;
}
if (Flag & Mask) {
Count++;
if (Cand % Div == 0) {
return Div;
}
}
Div += Wheel[Ind];
Mask = Mask << 1;
}
EInd++;
}
// Start the Wheel after encoded primes
while (Div <= Top) {
for (long long Ind = 0; Wheel[Ind] != 0; Ind++) {
if (Div > Top) {
return Cand;
}
Count++;
if (Cand % Div == 0) {
return Div;
}
Div += Wheel[Ind];
}
}
return Cand;
}
TrialDivLP.cpp 在 TrialDivLP.zip 中.
编码素数列表
30轮（Wheel 30）有 8 个槽，一个字节有 8 位，所以两者是同步的，这简化了它的使用。每个槽分配一个位。一个 char 的值是这些槽中是素数的槽位的总和。
列表的编码方式如下
0x01 0x02 0x04 0x08 0x10 0x20 0x40 0x80 code
7 11 13 17 19 23 29 0xfe
31 37 41 43 47 53 59 0xdf
61 67 71 73 79 83 89 0xef
97 101 103 107 109 113 0x7e
127 131 137 139 149 0xb6
151 157 163 167 173 179 0xdb
181 191 193 197 199 0x3d
211 223 227 229 233 239 0xf9
241 251 257 263 269 0xd5
271 277 281 283 293 0x4f
307 311 313 317 0x1e
331 337 347 349 353 359 0xf3
367 373 379 383 389 0xea
397 401 409 419 0xa6
421 431 433 439 443 449 0xed
457 461 463 467 479 0x9e
487 491 499 503 509 0xe6
521 523 0x0c
541 547 557 563 569 0xd3
571 577 587 593 599 0xd3
显示列表编码方式的代码如下
void TD_EncodeDtl(int Max) {
long long Wheel[] = { 6, 4, 2, 4, 2, 4, 6, 2, 0 };
long long Cand = 1;
cout << "0x01\t0x02\t0x04\t0x08\t0x10\t0x20\t0x40\t0x80\tcode" << endl;
while (Cand < Max) {
int Code = 0;
int Index = 1;
int Ind = 0;
do {
if ((TD_SRW(Cand) == Cand) && (Cand != 1)) {
Code |= Index;
cout << dec << Cand;
}
cout << "\t";
Cand += Wheel[Ind];
Index = Index << 1;
Ind++;
} while (Wheel[Ind] != 0);
cout << "0x" << setfill('0') << setw(2) << hex << Code << dec << endl;
}
}
TrialDivLP.cpp 在 TrialDivLP.zip 中.
生成数组的代码如下
// encode Prime Numbers as C array
void TD_EncodeArray(int Max) {
long long Wheel[] = { 6, 4, 2, 4, 2, 4, 6, 2, 0 };
long long Cand = 1;
cout << "Encodage liste nombres premiers compressée" << endl;
cout << "{";
while (Cand < Max) {
int Code = 0;
int Index = 1;
int Ind = 0;
do {
if ((TD_SRW(Cand) == Cand) && (Cand != 1)) {
Code |= Index;
}
Cand += Wheel[Ind];
Index = Index << 1;
Ind++;
} while (Wheel[Ind] != 0);
cout << "0x";
cout << setfill('0') << setw(2) << hex << Code << ",";
}
cout << "0}" << dec << endl;
}
TrialDivLP.cpp 在 TrialDivLP.zip 中.
基准测试
我们计算除法/取模次数来衡量变体的效率。编码的素数列表高达 1,000,000，存储在一个 33kB 的数组中。
Number TD_SR TD_SREF TD_SRW TD_SRLPW Delta
1,009 30 16 11 11 0
2,003 43 22 14 14 0
3,001 53 27 17 16 1
4,001 62 32 19 18 1
5,003 69 35 20 19 1
6,007 76 39 23 21 2
7,001 82 42 25 23 2
8,009 88 45 26 24 2
9,001 93 47 27 24 3
10,007 99 50 28 25 3
20,011 140 71 40 34 6
30,011 172 87 49 40 9
40,009 199 100 56 46 10
49,999 222 112 62 48 14
1,000,003 999 500 268 168 100
4,000,021 1,800 901 483 279 204
9,000,011 2,999 1,500 802 430 372
16,000,057 3,999 2,000 1,068 550 518
25,000,009 4,999 2,500 1,336 669 667
36,000,007 5,999 3,000 1,602 783 819
49,000,027 6,999 3,500 1,868 900 968
64,000,031 7,999 4,000 2,136 1,007 1,129
81,000,001 8,999 4,500 2,402 1,117 1,285
100,000,007 9,999 5,000 2,668 1,229 1,439
121,000,003 10,999 5,500 2,936 1,335 1,601
144,000,001 11,999 6,000 3,202 1,438 1,764
169,000,033 12,999 6,500 3,468 1,547 1,921
196,000,001 13,999 7,000 3,736 1,652 2,084
225,000,011 14,999 7,500 4,002 1,754 2,248
256,000,001 15,999 8,000 4,268 1,862 2,406
289,000,001 16,999 8,500 4,536 1,960 2,576
10,000,000,019 99,999 50,000 26,668 9,592 17,076
1,000,000,000,039 999,999 500,000 266,668 78,498 188,170
Delta 用于突出庞大素数列表带来的节省。
long long Test[] = { 101, 1009, 2003, 3001, 4001, 5003, 6007, 7001, 8009,
9001, 10007, 20011, 30011, 40009, 49999, 1000003, 4000021, 9000011,
16000057, 25000009, 36000007, 49000027, 64000031, 81000001,
100000007, 121000003, 144000001, 169000033, 196000001, 225000011,
256000001, 289000001, 10000000019, 1000000000039, 0 };
cout << "Comparaison Variantes" << endl;
cout << "Number\tTD_SR\tTD_SREF\tTD_SRW\tTD_SRLPW\tDelta" << endl;
// test after first 0
for (Ind = 0; Test[Ind] != 0; Ind++) {
cout << Test[Ind] << "\t";
TD_SR(Test[Ind]);
cout << Count << "\t";
TD_SREF(Test[Ind]);
cout << Count << "\t";
TD_SRW(Test[Ind]);
int C1 = Count;
cout << Count << "\t";
TD_SRLPW(Test[Ind]);
cout << Count << "\t";
cout << C1 - Count << endl;
}
cout << endl;
TrialDivLP.cpp 在 TrialDivLP.zip 中.
128位整数
GCC 通过使用 2 个 64 位寄存器来提供 128 位变量的支持，这对于计算来说没问题，但不支持作为源代码中的常量或用于 I/O 操作。需要一些编程来绕过这些限制。
//#define int64
// GCC have provisions to handle 128 bits integers, but not for constants in source code
// and not IO
#ifdef int64
// 64 bits integer
typedef long long MyBigInt;
void output(MyBigInt Num, int Size) {
cout << setfill(' ') << setw(Size);
}
#else
// 128 bits integer
typedef __int128 MyBigInt;
void output(MyBigInt Num, int Size) {
if (Num <= (long long) 0x7fffffffffffffff) {
if (Size > 0) {
cout << setfill(' ') << setw(Size);
}
cout << (long long) Num;
} else {
output(Num / 1000, Size - 4);
long long Tmp = Num % 1000;
cout << "," << setfill('0') << setw(3) << Tmp;
}
}
#endif
TrialDivLP.cpp 在 TrialDivLP.zip 中.
Number Fartor Timing
10,000,000,000,000,061 1 365 ms
100,000,000,000,000,003 1 1,173 ms
1,000,000,000,000,000,003 1 3,866 ms
10,000,000,000,000,000,051 1 11,988 ms
100,000,000,000,000,000,039 1 82,326 ms
1,000,000,000,000,000,000,117 1 272,376 ms
每次数字乘以 10，运行时间大约乘以 3（sqrt(10)）。注意当数字跨越 64 位边界时运行时间的差距。
#ifndef int64
long long Test2[][3] = { { 10, 16, 61 }, { 10, 17, 3 }, { 10, 18, 3 }, { 10,
19, 51 }, { 10, 20, 39 }, { 10, 21, 117 }, { 0, 0, 0 } };
// for each triplet A B C, the number is
// Cand= A ^ B + C
cout.imbue(locale(cout.getloc(), new gr3));
for (Ind = 0; Test2[Ind][0] > 0; Ind++) {
for (int Offset = 0; Offset < 1000; Offset += 2) {
MyBigInt Cand = 1;
for (int Scan = 0; Scan < Test2[Ind][1]; Scan++) {
Cand *= Test2[Ind][0];
}
Cand += Test2[Ind][2] + Offset;
output(Cand, 30);
cout << "\t";
auto start2 = high_resolution_clock::now();
long long Tmp = TD_SRLPW(Cand);
auto stop2 = high_resolution_clock::now();
if (Tmp == 1) {
output(Tmp, 0);
} else {
cout << Tmp;
cout << "\t";
output((Cand / Tmp) * Tmp, 0);
}
cout << "\t";
auto duration = duration_cast<milliseconds>(stop2 - start2);
cout << " " << duration.count() << " ms";
cout << endl;
if (Tmp == 1) {
break;
}
}
}
cout.imbue(locale(cout.getloc(), new gr0));
#endif
TrialDivLP.cpp 在 TrialDivLP.zip 中.
使用压缩的素数列表加速 IsPrime 判断
我们可以利用庞大的素数列表并将其用作查找表。
使用普通列表，我们需要进行二分查找。
使用编码列表，几乎可以进行直接读取。
// IsPrime taking advantage of long list of primes table
// returns 1 if prime, returns 0 otherwise
int IsPrimeLP(MyBigInt Cand) {
const unsigned long long Wheel[] = { 6, 4, 2, 4, 2, 4, 6, 2, 0 };
const unsigned long long WOffset[] = { 1, 7, 11, 13, 17, 19, 23, 29, 0 };
if (Cand <= EPrimesMax) {
Count = 1;
// Cand is within list of primes
long long x = Cand / WSize; // position in wheel list
long long y = Cand % WSize; // offset in wheel
for (int Index = 0; WOffset[Index] != 0; Index++) {
if (WOffset[Index] == y) {
return ((EPrimes[x] >> Index) & 0x01);
}
}
if (Cand == 2)
return 1;
else if (Cand == 3)
return 1;
else if (Cand == 5)
return 1;
else
return 0; // multiple of 2, 3 or 5
}
// Search a Factor
return (TD_SRLPW(Cand) == 1);
}
TrialDivLP.cpp 在 TrialDivLP.zip 中.
关注点
随着素数列表的增大，试除法算法会得到改进，而压缩技术在这方面帮助很大。
这项工作是原创的。如果您知道类似的工作，请在下面的论坛中提供链接。
链接
历史
- 2020年12月18日：初版
- 2020年12月25日：第二版 - 修正和改进
- 2021年2月6日：添加 IsPrime 代码 + 修正
- 2022年8月7日：添加目录
- 2022年8月20日：源代码中的拼写错误