整数分解：试除算法
试除算法的小回顾
目录
背景
本文是系列文章的一部分。每篇文章都着重介绍整数分解的某个方面。
- 整数分解：试除法算法[^] : 侧重于试除法算法的变体及其各自的效率。
- 整数分解：令人头疼的素数列表[^] : 侧重于通过压缩处理大型素数列表的方法。
- 整数分解：优化小因子检查[^] : 侧重于一种比除法更快地检查小因子的方法。
- 整数分解：反转乘法[^] : 侧重于优化试除法算法的另一种方法。
xPrompt 下载和 prg 文件
xPrompt.zip 包含 xHarbour Scripting Engine for Windows 的二进制文件。这是运行 prg 代码所需的应用程序。
用法
- 将 prg 文件复制到 xPrompt 的同一目录下
- 启动 xPrompt
- 输入 "
do FileName.prg" 来运行代码
xHarbour 是免费软件，可在 xHarbour.org[^] 下载
xHarbour 是 xBase 系列语言的一部分：xBase - Wikipedia[^]
流程图
流程图与本文相关：使用树形图支持的自顶向下分析方法[^]
如果图片太小，请在新标签页中打开并缩放，这是一个 svg 文件。
引言
在开始玩弄整数分解时，尝试所有可能的因子是第一个想法，该算法称为 **试除法**。
该算法有两个目的——找到一个素数因子，或者通过不找到素数因子来判断一个整数是否为素数。
由于该算法是关于查找因子的，所以最坏的情况是待分解的整数是素数。
试除法：经典变体
很明显，最有效的方法是只检查素数，但处理素数列表本身就是一个问题。经典的变体是为了避免这个问题而存在的。随着方法的不断完善，它们会去除更多无用的非素数，从而更有效。
暴力破解
一些新手会测试小于 n 的所有单个数字。这是多余的。只要待分解的整数不是素数，低效性就会隐藏起来。
n 的最大成本是 n-2 次除法。复杂度为 O(n)。
// Trial Division: Brute Force 1
// Check all numbers until Cand - 1
long long TD_BF1(long long Cand) {
Count = 0;
long long Top = Cand - 1;
for (long long Div = 2; Div <= Top; Div++) {
Count++;
if (Cand % Div == 0) {
return Div;
}
}
return Cand;
}
// Trial Division Brute Force 1
// May 2019
function TD_BF1(Prod)
Local D, Top
Top= Prod-1
for D= 2 to Top
if Prod % D = 0
return D
endif
next
return Prod
在 TrialDiv.cpp 和 TrialDiv.prg 中，位于 TrialDiv.zip 内.
TD_BF1.graphml 和 TD_BF1.svg 在 TrialDiv.zip 内.
测试直到 n/2 的所有数字会更好，但同样是多余的。就像以前的方法一样，只要待分解的整数不是素数，低效性就会隐藏起来。
n 的最大成本是 n / 2 次除法。复杂度为 O(n)。
// Trial Division: Brute Force 2
// Check all numbers until Cand / 2
long long TD_BF2(long long Cand) {
Count = 0;
long long Top = Cand / 2;
for (long long Div = 2; Div <= Top; Div++) {
Count++;
if (Cand % Div == 0) {
return Div;
}
}
return Cand;
}
// Trial Division Brute Force 2
// May 2019
function TD_BF2(Prod)
Local D, Top
Top= Prod/2
for D= 2 to Top
if Prod % D = 0
return D
endif
next
return Prod
在 TrialDiv.cpp 和 TrialDiv.prg 中，位于 TrialDiv.zip 内.
TD_BF2.graphml 和 TD_BF2.svg 在 TrialDiv.zip 内.
基本版本：平方根
稍微分析一下就知道，平方根之后就没有需要检查的因子了。
n 的最大成本是 √n 次除法。复杂度为 O(√n)。
// Trial Division: Square Root
// Check all numbers until Square Root
long long TD_SR(long long Cand) {
Count = 0;
long long Top = sqrt(Cand);
for (long long Div = 2; Div <= Top; Div++) {
Count++;
if (Cand % Div == 0) {
return Div;
}
}
return Cand;
}
// Trial Division Square Root
// May 2019
function TD_SR(Prod)
Local D, Top
Top= int(sqrt(Prod))
for D= 2 to Top
if Prod % D = 0
return D
endif
next
return Prod
在 TrialDiv.cpp 和 TrialDiv.prg 中，位于 TrialDiv.zip 内.
TD_SR.graphml 和 TD_SR.svg 在 TrialDiv.zip 内.
去除偶数因子
进一步分析可知，除了 2 之外，不会有偶数因子。
n 的最大成本是 (√n) * 1 / 2 => (√n) * 0.50 次除法。复杂度为 O(√n)。
// Trial Division: Square Root + Even Factors
// Check all numbers until Square Root and Skip Even Factors
long long TD_SREF(long long Cand) {
Count = 0;
// check 2 the only Even Prime
Count++;
if (Cand % 2 == 0) {
return 2;
}
long long Top = sqrt(Cand);
for (long long Div = 3; Div <= Top; Div += 2) {
Count++;
if (Cand % Div == 0) {
return Div;
}
}
return Cand;
}
// Trial Division Square Root + Even Factors
// May 2019
function TD_SREF(Prod)
Local D, Top
if Prod % 2 = 0
return 2
endif
Top= int(sqrt(Prod))
for D= 3 to Top step 2
if Prod % D = 0
return D
endif
next
return Prod
在 TrialDiv.cpp 和 TrialDiv.prg 中，位于 TrialDiv.zip 内.
TD_SREF.graphml 和 TD_SREF.svg 在 TrialDiv.zip 内.
轮子
轮子是先前优化的扩展。轮子由小的素数（例如 2 和 3）构建，轮子的大小是 2 * 3 = 6。当在一张 6 列的表格中写下一系列数字并去掉 2 和 3 时，可以看出素数只出现在第一列和第五列。这就是轮子，我们只需要检查第 1 列和第 5 列的数字。
| 1 | 2 | 3 | 4 | 5 | 6 |
| 7 | 8 | 9 | 10 | 11 | 12 |
| 13 | 14 | 15 | 16 | 17 | 18 |
| 19 | 20 | 21 | 22 | 23 | 24 |
| 25 | 26 | 27 | 28 | 29 | 30 |
| 31 | 32 | 33 | 34 | 35 | 36 |
| 37 | 38 | 39 | 40 | 41 | 42 |
| 43 | 44 | 45 | 46 | 47 | 48 |
| 49 | 50 | 51 | 52 | 53 | 54 |
| 55 | 56 | 57 | 58 | 59 | 60 |
n 的最大成本是 (√n) * (1 / 2) * (2 / 3) => (√n) * (1 / 3) => (√n) * 0.33 次除法。复杂度为 O(√n)。
// Trial Division Square Root + Wheel
// May 2019
<pre lang="c++">
// Trial Division: Square Root + Wheel
// Check all numbers until Square Root and Skip useless factors with a wheel
long long TD_SRW(long long Cand) {
long long SPrimes[] = { 2, 3, 0 };
long long Wheel[] = { 4, 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 = 1;
while (Div <= Top) {
for (long long Ind = 0; Wheel[Ind] != 0; Ind++) {
Div += Wheel[Ind];
if (Div > Top) {
break;
}
Count++;
if (Cand % Div == 0) {
return Div;
}
}
}
return Cand;
}
function TD_SRW(Prod)
local D, Top, SPrimes, Wheel, W
// Check small primes
SPrimes= {2, 3}
Wheel= {4, 2}
for each D in SPrimes
if Prod % D = 0
return D
endif
next
// Start the wheel
D= 1
Top= int(sqrt(Prod))
while D <= Top
for each W in wheel
D += W
if Prod % D = 0
return D
endif
next
enddo
return Prod
在 TrialDiv.cpp 和 TrialDiv.prg 中，位于 TrialDiv.zip 内.
TD_SRW.graphml 和 TD_SRW.svg 在 TrialDiv.zip 内.
轮子可以包含更多素数。对于 2、3 和 5，轮子的大小是 30。
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 |
| 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 |
| 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 |
n 的最大成本是 (√n) * (1 / 2) * (2 / 3) * (4 / 5) => (√n) * (4 / 15) => (√n) * 0.27 次除法。复杂度为 O(√n)。
long long SPrimes[] = { 2, 3, 5, 0 };
long long Wheel[] = { 6, 4, 2, 4, 2, 4, 6, 2, 0 };
// only those 2 lines need change
SPrimes= {2, 3, 5}
Wheel= {6, 4, 2, 4, 2, 4, 6, 2}
或者更大。对于 2、3、5 和 7，轮子的大小是 210。
n 的最大成本是 (√n) * (1 / 2) * (2 / 3) * (4 / 5) * (6 / 7) => (√n) * (24 / 105) => (√n) * 0.23 次除法。复杂度为 O(√n)。
素数列表
基本上，它是一个轮子变体，其中小素数列表超出了轮子所需的数量，当列表结束时，我们会回到轮子。与简单轮子的优点是，它避免了测试轮子未过滤的非素数因子。只要我们在素数列表中，工作量就是最佳的。
只需注意在素数列表末尾正确设置轮子。
n 的最大成本略好于轮子，素数列表只是节省了轮子中的非素数除数。因此，素数列表越长，效果越好。复杂度为 O(√n)。
// 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, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661,
673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757,
761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853,
857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021,
1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093,
1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259,
1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433,
1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831,
1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913,
1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087,
2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269,
2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347,
2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417,
2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,
2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767,
2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953,
2957, 2963, 2969, 2971, 2999, 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 = 3001;
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;
}
// Trial Division Square Root + Wheel + list of primes
// with a large list of primes
// May 2019
function TD_SRW2(Prod)
local D, Top, SPrimes, Wheel, W
// Check small primes
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}
Wheel= {6, 4, 2, 4, 2, 4, 6, 2}
for each D in SPrimes
if Prod % D = 0
return D
endif
next
// Start the wheel
D= 301
Top= int(sqrt(Prod))
while D <= Top
for each W in wheel
D += W
if Prod % D = 0
return D
endif
next
enddo
return Prod
在 TrialDiv.cpp 和 TrialDiv.prg 中，位于 TrialDiv.zip 内.
基准测试
基准测试是为了突出变体的效率。
由于除法/模运算是耗时的，因此计算它们足以进行基准测试。
基准测试 1
第一个基准测试是为了强调暴力破解变体的低效性。
long long Test[] = { 11, 31, 53, 71, 97, 0, 101, 1009, 2003, 3001, 4001,
5003, 6007, 7001, 8009, 9001, 10007, 20011, 30011, 40009, 49999,
1000003, 4000021, 9000011, 0 };
cout << "Comparaison Variantes avec Brute Force" << endl;
cout << "Number\tTD_BF1\tTD_BF2\tTD_SR\tTD_SREF\tTD_SRW\tTD_SRPW" << endl;
for (Ind = 0; Test[Ind] != 0; Ind++) {
cout << Test[Ind] << "\t";
TD_BF1(Test[Ind]);
cout << Count << "\t";
TD_BF2(Test[Ind]);
cout << Count << "\t";
TD_SR(Test[Ind]);
cout << Count << "\t";
TD_SREF(Test[Ind]);
cout << Count << "\t";
TD_SRW(Test[Ind]);
cout << Count << "\t";
TD_SRPW(Test[Ind]);
cout << Count << endl;
}
cout << endl;
Test = { 11, 31, 53, 71, 97, 101, 1009, 2003, 3001, 4001, 5003, 6007, 7001, 8009, 9001, 10007, 20011, 30011, 40009, 49999, 1000003, 4000021, 9000011 }
? "Comparaison Variantes avec Brute Force"
? "Number","TD_BF1","TD_BF2","TD_SR","TD_SREF","TD_SRW","TD_SRPW"
for Ind = 1 to 5
? Test[Ind]
TD_BF1(Test[Ind])
?? Count
TD_BF2(Test[Ind])
?? Count
TD_SR(Test[Ind])
?? Count
TD_SREF(Test[Ind])
?? Count
TD_SRW(Test[Ind])
?? Count
TD_SRPW(Test[Ind])
?? Count
next
?
Number TD_BF1 TD_BF2 TD_SR TD_SREF TD_SRW TD_SRPW
11 9 4 2 2 2 2
31 29 14 4 3 3 3
53 51 25 6 4 4 4
71 69 34 7 4 4 4
97 95 47 8 5 4 4
基准测试 2
这个扩展基准测试用于暴力破解之外的变体。
cout << "Comparaison Variantes sans Brute Force" << endl;
cout << "Number\tTD_SR\tTD_SREF\tTD_SRW\tTD_SRPW\tDelta" << endl;
// test after first 0
for (Ind++; 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_SRPW(Test[Ind]);
cout << Count << "\t";
cout << C1 - Count << endl;
}
cout << endl;
? "Number","TD_SR","TD_SREF","TD_SRW","TD_SRPW"
for Ind = 6 to len(Test)
? Test[Ind]
TD_SR(Test[Ind])
?? Count
TD_SREF(Test[Ind])
?? Count
TD_SRW(Test[Ind])
?? Count
TD_SRPW(Test[Ind])
?? Count
next
?
Comparison Variants sans Brute Force
Number TD_SR TD_SREF TD_SRW TD_SRPW Delta
101 9 5 4 4 0
1009 30 16 11 11 0
2003 43 22 14 14 0
3001 53 27 17 16 1
4001 62 32 19 18 1
5003 69 35 20 19 1
6007 76 39 23 21 2
7001 82 42 25 23 2
8009 88 45 26 24 2
9001 93 47 27 24 3
10007 99 50 28 25 3
20011 140 71 40 34 6
30011 172 87 49 40 9
40009 199 100 56 46 10
49999 222 112 62 48 14
1000003 999 500 268 168 100
4000021 1800 901 483 279 204
9000011 2999 1500 802 430 372
Delta 显示了素数列表如何改进性能，并且随着素数列表的增长而变得更好。
检查变体正确性
Rgis 代码通过比较变体结果来检查代码的正确性。
cout << "Vérification Variantes" << endl;
for (long long Cand = 3; Cand < 10000000; Cand += 10) {
long long d1 = TD_SREF(Cand);
long long d2 = TD_SRPW(Cand);
if (d1 != d2) {
cout << Cand << "\t" << d1 << "\t" << d2 << endl;
}
}
关注点
试除法是暴力破解，因此可以看出暴力破解也有好坏之分。
链接
- 整数因式分解：整数因式分解 - 维基百科[^]
- 试除法：试除法 - Wikipedia[^]
- 暴力搜索：暴力搜索 - Wikipedia[^]
历史
- 2019年4月1日：第一个版本
- 2020年11月27日：第二个版本
- 2020年12月20日：第三个版本：一些清理，将下载移至顶部
- 2020年12月25日：第四个版本：一些清理和更正
- 2021年1月31日：更正
- 2022年1月30日：添加目录和少量更新
- 2023年9月16日：添加流程图和更正。
- 2020年12月15日：拼写错误