蛮力，单线程素数分解
Posted
技术标签:
【中文标题】蛮力，单线程素数分解【英文标题】：Brute-force, single-threaded prime factorization 【发布时间】：2010-10-12 20:57:31 【问题描述】：需要考虑的是以下函数，该函数可用于（相对快速地）将 64 位无符号整数分解为其质因数。请注意，因式分解不是概率的（即，它是精确的）。在现代硬件上，该算法已经足够快，可以在几秒钟内找到一个数是素数或几乎没有非常大的因数。
问题： 可以对所提出的算法进行任何改进，同时保持单线程，以便它可以更快地分解（任意）非常大的无符号 64 位整数，最好不使用概率方法（例如 Miller-Rabin）来确定素数?
// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef unsigned __int64 ulong;
typedef std::vector<ulong> ULongVector;
// Caller needs to pass in an empty factors vector
void GetFactors(ULongVector &factors, ulong num)
// Num has to be at least 2 to contain "prime" factors
if (num<2)
return;
ulong workingNum=num;
ulong nextOffset=2; // Will be used to skip multiples of 3, later
// Factor out factors of 2
while (workingNum%2==0)
factors.push_back(2);
workingNum/=2;
// Factor out factors of 3
while (workingNum%3==0)
factors.push_back(3);
workingNum/=3;
// If all of the factors were 2s and 3s, done...
if (workingNum==1)
return;
// sqrtNum is the (inclusive) upper bound of our search for factors
ulong sqrtNum=(ulong) sqrt(double(workingNum+0.5));
// Factor out potential factors that are greate than or equal to 5
// The variable n represents the next potential factor to be tested
for (ulong n=5;n<=sqrtNum;)
// Is n a factor of the current working number?
if (workingNum%n==0)
// n is a factor, so add it to the list of factors
factors.push_back(n);
// Divide current working number by n, to get remaining number to factor
workingNum/=n;
// Check if we've found all factors
if (workingNum==1)
return;
// Recalculate the new upper bound for remaining factors
sqrtNum=(ulong) sqrt(double(workingNum+0.5));
// Recheck if n is a factor of the new working number,
// in case workingNum contains multiple factors of n
continue;
// n is not or is no longer a factor, try the next odd number
// that is not a multiple of 3
n+=nextOffset;
// Adjust nextOffset to be an offset from n to the next non-multiple of 3
nextOffset=(nextOffset==2UL ? 4UL : 2UL);
// Current workingNum is prime, add it as a factor
factors.push_back(workingNum);
谢谢
编辑：我添加了更多的 cmets。通过引用传入向量的原因是允许在调用之间重用向量并避免动态分配。向量没有在函数中清空的原因是允许将当前“num's”因子附加到向量中已经存在的因子上的奇怪要求。
函数本身并不漂亮，可以重构，但问题是如何让算法更快。所以，请不要提出关于如何使函数更漂亮、更易读或 C++ish 的建议。那是儿戏。改进这个算法，使其能够更快地找到（证明）因子是困难的部分。
更新：Potatoswatter 到目前为止有一些出色的解决方案，请务必在底部附近查看他的 MMX 解决方案。
【问题讨论】：
你的问题是什么？ 问题是如何加速这个算法，希望通过算法改进/添加。 在不损失效率的情况下，您可以将其设为接受输出迭代器的模板，从而消除对vector 的依赖。然后，用户仍然可以使用 vector 和 GetFactors(n, back_inserter(v))。
对于最难的 64 位无符号整数（例如，具有非常大因数的素数和复合数），Fast 大约需要 6 秒。我想知道这个简单的算法是否可以改进得更快。
@Michael：您可能会尝试展开循环一次，这样您就没有条件 +2 或 +4。相反，您在循环中间无条件地增加 2，最后增加 4。可能没有任何区别，可能会在 50-100% 的时间内防止分支错误预测。谁知道呢。
【参考方案1】：
将这种方法与（预先生成的）筛子进行比较。取模很昂贵，所以这两种方法本质上都做两件事：生成潜在因子，并执行取模运算。任何一个程序都应该在比模数更少的周期内合理地生成一个新的候选因子，因此任何一个程序都是模数限制的。
给定的方法过滤掉所有整数的恒定比例，即 2 和 3 的倍数，即 75%。四分之一（如给定的）数字用作模运算符的参数。我称之为跳过过滤器。
另一方面，筛子只使用素数作为模运算符的参数，difference between successive primes 的平均值由prime number theorem 控制为 1/ln(N)。例如，e^20 略低于 5 亿，因此超过 5 亿的数字有不到 5% 的机会成为素数。如果考虑到 2^32 以内的所有数字，则 5% 是一个很好的经验法则。
因此，作为您的跳过过滤器，筛子在 div 操作上花费的时间将减少 5 倍。下一个要考虑的因素是筛子产生素数的速度，即从内存或磁盘中读取它们。如果获取一个素数比 4 divs 快，则筛子更快。根据我的表格div，Core2 上的吞吐量最多为每 12 个周期一个。这些将是硬除法问题，所以让我们保守地预算每个素数 50 个周期。对于 2.5 GHz 处理器，这是 20 纳秒。
在 20 ns 内，一个 50 MB/秒的硬盘可以读取大约一个字节。简单的解决方案是每个素数使用 4 个字节，因此驱动器会更慢。但是，我们可以更聪明。如果我们想按顺序编码所有素数，我们可以只编码它们的差异。同样，预期的差异是 1/ln(N)。而且，它们都是偶数，这节省了额外的一点。而且它们永远不会为零，这使得对多字节编码的扩展是免费的。因此，每个素数使用一个字节，最多可以将 512 的差异存储在一个字节中，根据that Wikipedia article，我们最多可以存储 303371455241。
因此，根据硬盘驱动器的不同，存储的素数列表在验证素数时的速度应该大致相同。如果它可以存储在 RAM 中（它是 203 MB，因此后续运行可能会遇到磁盘缓存），那么问题就完全消失了，因为 FSB 速度通常与处理器速度的差异小于 FSB 宽度（以字节为单位） — 即，FSB 每个周期可以传输多个素数。然后改进的因素是减少除法操作，即五倍。下面的实验结果证明了这一点。
当然，还有多线程。可以将素数或跳过过滤的候选范围分配给不同的线程，从而使这两种方法都令人尴尬地并行。没有不涉及增加并行除法器电路数量的优化，除非您以某种方式消除模数。
这是一个这样的程序。它是模板化的，因此您可以添加 bignums。
/*
* multibyte_sieve.cpp
* Generate a table of primes, and use it to factorize numbers.
*
* Created by David Krauss on 10/12/10.
*
*/
#include <cmath>
#include <bitset>
#include <limits>
#include <memory>
#include <fstream>
#include <sstream>
#include <iostream>
#include <iterator>
#include <stdint.h>
using namespace std;
char const primes_filename[] = "primes";
enum encoding_base = (1<< numeric_limits< unsigned char >::digits) - 2 ;
template< typename It >
unsigned decode_gap( It &stream )
unsigned gap = static_cast< unsigned char >( * stream ++ );
if ( gap ) return 2 * gap; // only this path is tested
gap = ( decode_gap( stream )/2-1 ) * encoding_base; // deep recursion
return gap + decode_gap( stream ); // shallow recursion
template< typename It >
void encode_gap( It &stream, uint32_t gap )
unsigned len = 0, bytes[4];
gap /= 2;
do
bytes[ len ++ ] = gap % encoding_base;
gap /= encoding_base;
while ( gap );
while ( -- len ) // loop not tested
* stream ++ = 0;
* stream ++ = bytes[ len + 1 ];
* stream ++ = bytes[ 0 ];
template< size_t lim >
void generate_primes()
auto_ptr< bitset< lim / 2 > > sieve_p( new bitset< lim / 2 > );
bitset< lim / 2 > &sieve = * sieve_p;
ofstream out_f( primes_filename, ios::out | ios::binary );
ostreambuf_iterator< char > out( out_f );
size_t count = 0;
size_t last = sqrtl( lim ) / 2 + 1, prev = 0, x = 1;
for ( ; x != last; ++ x )
if ( sieve[ x ] ) continue;
size_t n = x * 2 + 1; // translate index to number
for ( size_t m = x + n; m < lim/2; m += n ) sieve[ m ] = true;
encode_gap( out, ( x - prev ) * 2 );
prev = x;
for ( ; x != lim / 2; ++ x )
if ( sieve[ x ] ) continue;
encode_gap( out, ( x - prev ) * 2 );
prev = x;
cout << prev * 2 + 1 << endl;
template< typename I >
void factorize( I n )
ifstream in_f( primes_filename, ios::in | ios::binary );
if ( ! in_f )
cerr << "Could not open primes file.\n"
"Please generate it with 'g' command.\n";
return;
while ( n % 2 == 0 )
n /= 2;
cout << "2 ";
unsigned long factor = 1;
for ( istreambuf_iterator< char > in( in_f ), in_end; in != in_end; )
factor += decode_gap( in );
while ( n % factor == 0 )
n /= factor;
cout << factor << " ";
if ( n == 1 ) goto finish;
cout << n;
finish:
cout << endl;
int main( int argc, char *argv[] )
if ( argc != 2 ) goto print_help;
unsigned long n;
if ( argv[1][0] == 'g' )
generate_primes< (1ul<< 32) >();
else if ( ( istringstream( argv[1] ) >> n ).rdstate() == ios::eofbit )
factorize( n );
else goto print_help;
return 0;
print_help:
cerr << "Usage:\n\t" << argv[0] << " <number> -- factorize number.\n"
"\t" << argv[0] << " g -- generate primes file in current directory.\n";
在 2.2 GHz MacBook Pro 上的性能：
dkrauss$ time ./multibyte_sieve g
4294967291
real 2m8.845s
user 1m15.177s
sys 0m2.446s
dkrauss$ time ./multibyte_sieve 18446743721522234449
4294967231 4294967279
real 0m5.405s
user 0m4.773s
sys 0m0.458s
dkrauss$ time ./mike 18446743721522234449
4294967231 4294967279
real 0m25.147s
user 0m24.170s
sys 0m0.096s
【讨论】：
嗯，这个理论听起来确实很有趣，但最终，让我们看看它的表现如何。因此，当您修复错误时，请在此处发布您的结果。我特别想知道它如何处理非常大的素数（例如 18446744073709542713），以及它如何处理仅包含两个大素数因子的大数（例如 18446743721522234449）。 ……我需要 25 秒。哇，你的机器速度很快。 所以结论是我们可以期待的：如果你想偶尔分解一个大数字，mike 的代码看起来更好，但如果你对分解许多数字感兴趣并且磁盘空间不是一个问题，筛子要快得多（20位数字快五倍，这是两个10位素数的乘积）。比较恒定成本（生成素数文件）与一次因式分解的成本乘以调用次数。 @wok：差不多。不过，拥有这张桌子的开销非常低。生成该表需要 75 秒，不包括 I/O，但包括一些操作系统开销。我没有在常规数组上尝试过素数生成器而不是streambuf_iterator，但如果你这样运行它应该会快一点。现在我正在尝试优化，这样它就不会过多地消耗内存。盈亏平衡点已经只是需要验证的 3 个素数，我可能会降低到 2 个。
我想知道对筛子的测试是否可以有效地并行化。【参考方案2】：
我的另一个答案很长，与这个答案完全不同，所以这里有别的东西。
这个程序不只是过滤掉前两个素数的倍数，或者将所有相关的素数编码到一个字节中，而是过滤掉所有适合 8 位的素数的倍数，特别是 2 到 211。所以不是通过33% 的数字，这会将大约 10% 传递给除法运算符。
素数保存在三个 SSE 寄存器中，它们与运行计数器的模数保存在另外三个中。如果任何素数与计数器的模数为零，则计数器不能是素数。此外，如果任何模数等于 1，则 (counter+2) 不能是素数，依此类推，直到 (counter+30)。偶数会被忽略，并且像 +3、+6 和 +5 这样的偏移会被跳过。矢量处理允许一次更新 16 个字节大小的变量。
在对厨房水槽进行了全面的微优化（但没有比内联指令更具体的平台）之后，我的性能提升了 1.78 倍（我的笔记本电脑为 24 秒 vs 13.4 秒）。如果使用 bignum 库（甚至是非常快的库），则优势更大。请参阅下面的更易读的预优化版本。
/*
* factorize_sse.cpp
* Filter out multiples of the first 47 primes while factorizing a number.
*
* Created by David Krauss on 10/14/10.
*
*/
#include <cmath>
#include <sstream>
#include <iostream>
#include <xmmintrin.h>
using namespace std;
inline void remove_factor( unsigned long &n, unsigned long factor ) __attribute__((always_inline));
void remove_factor( unsigned long &n, unsigned long factor )
while ( n % factor == 0 )
n /= factor;
cout << factor << " ";
int main( int argc, char *argv[] )
unsigned long n;
if ( argc != 2
|| ( istringstream( argv[1] ) >> n >> ws ).rdstate() != ios::eofbit )
cerr << "Usage: " << argv[0] << " <number>\n";
return 1;
int primes[] = 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 ;
for ( int *p = primes; p < primes + sizeof primes/sizeof *primes; ++ p )
remove_factor( n, * p );
//int histo[8] = , total = 0;
enum bias = 15 - 128 ;
__m128i const prime1 = _mm_set_epi8( 21, 21, 21, 22, 22, 26, 26, 17, 19, 23, 29, 31, 37, 41, 43, 47 ),
prime2 = _mm_set_epi8( 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127 ),
prime3 = _mm_set_epi8( 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211 ),
vbias = _mm_set1_epi8( bias ),
v3 = _mm_set1_epi8( 3+bias ), v5 = _mm_set1_epi8( 5+bias ), v6 = _mm_set1_epi8( 6+bias ), v8 = _mm_set1_epi8( 8+bias ),
v9 = _mm_set1_epi8( 9+bias ), v11 = _mm_set1_epi8( 11+bias ), v14 = _mm_set1_epi8( 14+bias ), v15 = _mm_set1_epi8( 15+bias );
__m128i mod1 = _mm_add_epi8( _mm_set_epi8( 3, 10, 17, 5, 16, 6, 19, 8, 9, 11, 14, 15, 18, 20, 21, 23 ), vbias ),
mod2 = _mm_add_epi8( _mm_set_epi8( 26, 29, 30, 33, 35, 36, 39, 41, 44, 48, 50, 51, 53, 54, 56, 63 ), vbias ),
mod3 = _mm_add_epi8( _mm_set_epi8( 65, 68, 69, 74, 75, 78, 81, 83, 86, 89, 90, 95, 96, 98, 99, 105 ), vbias );
for ( unsigned long factor = 1, limit = sqrtl( n ); factor <= limit + 30; factor += 30 )
if ( n == 1 ) goto done;
// up to 2^32, distribution of number candidates produced (0 up to 7) is
// 0.010841 0.0785208 0.222928 0.31905 0.246109 0.101023 0.0200728 0.00145546
unsigned candidates[8], *cand_pen = candidates;
* cand_pen = 6;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v3 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v3 ), _mm_cmpeq_epi8( mod3, v3 ) ) ) );
* cand_pen = 10;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v5 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v5 ), _mm_cmpeq_epi8( mod3, v5 ) ) ) );
* cand_pen = 12;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v6 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v6 ), _mm_cmpeq_epi8( mod3, v6 ) ) ) );
* cand_pen = 16;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v8 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v8 ), _mm_cmpeq_epi8( mod3, v8 ) ) ) );
* cand_pen = 18;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v9 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v9 ), _mm_cmpeq_epi8( mod3, v9 ) ) ) );
* cand_pen = 22;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v11 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v11 ), _mm_cmpeq_epi8( mod3, v11 ) ) ) );
* cand_pen = 28;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v14 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v14 ), _mm_cmpeq_epi8( mod3, v14 ) ) ) );
* cand_pen = 30;
cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v15 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v15 ), _mm_cmpeq_epi8( mod3, v15 ) ) ) );
/*++ total;
++ histo[ cand_pen - candidates ];
cout << "( ";
while ( cand_pen != candidates ) cout << factor + * -- cand_pen << " ";
cout << ")" << endl; */
mod1 = _mm_sub_epi8( mod1, _mm_set1_epi8( 15 ) ); // update residuals
__m128i mask1 = _mm_cmplt_epi8( mod1, _mm_set1_epi8( 1+bias ) );
mask1 = _mm_and_si128( mask1, prime1 ); // residual goes to zero or negative?
mod1 = _mm_add_epi8( mask1, mod1 ); // combine reset into zero or negative
mod2 = _mm_sub_epi8( mod2, _mm_set1_epi8( 15 ) );
__m128i mask2 = _mm_cmplt_epi8( mod2, _mm_set1_epi8( 1+bias ) );
mask2 = _mm_and_si128( mask2, prime2 );
mod2 = _mm_add_epi8( mask2, mod2 );
mod3 = _mm_sub_epi8( mod3, _mm_set1_epi8( 15 ) );
__m128i mask3 = _mm_cmplt_epi8( mod3, _mm_set1_epi8( 1+bias ) );
mask3 = _mm_and_si128( mask3, prime3 );
mod3 = _mm_add_epi8( mask3, mod3 );
if ( cand_pen - candidates == 0 ) continue;
remove_factor( n, factor + candidates[ 0 ] );
if ( cand_pen - candidates == 1 ) continue;
remove_factor( n, factor + candidates[ 1 ] );
if ( cand_pen - candidates == 2 ) continue;
remove_factor( n, factor + candidates[ 2 ] );
if ( cand_pen - candidates == 3 ) continue;
remove_factor( n, factor + candidates[ 3 ] );
if ( cand_pen - candidates == 4 ) continue;
remove_factor( n, factor + candidates[ 4 ] );
if ( cand_pen - candidates == 5 ) continue;
remove_factor( n, factor + candidates[ 5 ] );
if ( cand_pen - candidates == 6 ) continue;
remove_factor( n, factor + candidates[ 6 ] );
cout << n;
done:
/*cout << endl;
for ( int hx = 0; hx < 8; ++ hx ) cout << (float) histo[hx] / total << " ";*/
cout << endl;
.
dkrauss$ /usr/local/bin/g++ main.cpp -o factorize_sse -O3 --profile-use
dkrauss$ time ./factorize_sse 18446743721522234449
4294967231 4294967279
real 0m13.437s
user 0m13.393s
sys 0m0.011s
以下是上述内容的初稿。包括优化
使循环计数器无条件合并（避免分支）。 通过将循环展开 15 倍，将步幅增加到 30 来获得 ILP。 灵感来自您的优化。 30 似乎是一个最佳选择，因为它可以免费删除 2、3 和 5 的倍数。 5 到 15 之间的素数可能有多个倍数，因此请在向量中以不同的相位放置多个副本。 排除remove_factor。
将有条件的、不可预测的remove_factor 调用更改为非分支数组写入。
使用remove_factor 调用完全展开最终循环，并确保函数始终是内联的。
消除最终展开的迭代，因为候选中始终存在 7 的倍数。
添加另一个向量，其中包含所有剩余的足够小的素数。
通过向计数器添加偏差来腾出更多空间，并添加另一个向量。现在只有 6 个素数可以在不增加 16 位的情况下被过滤掉，而且我的寄存器也用完了：循环需要 3 个素数向量、3 个模数向量、8 个要搜索的常数，以及每个常数递增并进行范围检查。这使得 16。
在此应用中增益最小（但为正），但此技术的最初目的是在另一个答案中过滤筛子的素数。敬请期待……
可读版本：
/*
* factorize_sse.cpp
* Filter out multiples of the first 17 primes while factorizing a number.
*
* Created by David Krauss on 10/14/10.
*
*/
#include <cmath>
#include <sstream>
#include <iostream>
#include <xmmintrin.h>
using namespace std;
int main( int argc, char *argv[] )
unsigned long n;
if ( argc != 2
|| ( istringstream( argv[1] ) >> n >> ws ).rdstate() != ios::eofbit )
cerr << "Usage: " << argv[0] << " <number>\n";
return 1;
int primes[] = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 ;
for ( int *p = primes; p < primes + sizeof primes/sizeof *primes; ++ p )
while ( n % * p == 0 )
n /= * p;
cout << * p << " ";
if ( n != 1 )
__m128i mod = _mm_set_epi8( 1, 2, 3, 5, 6, 8, 9, 11, 14, 15, 18, 20, 21, 23, 26, 29 );
__m128i const prime = _mm_set_epi8( 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 ),
one = _mm_set1_epi8( 1 );
for ( unsigned long factor = 1, limit = sqrtl( n ); factor < limit; )
factor += 2;
__m128i mask = _mm_cmpeq_epi8( mod, one ); // residual going to zero?
mod = _mm_sub_epi8( mod, one ); // update other residuals
if ( _mm_movemask_epi8( mask ) )
mask = _mm_and_si128( mask, prime ); // reset cycle if going to zero
mod = _mm_or_si128( mask, mod ); // combine reset into zeroed position
else while ( n % factor == 0 )
n /= factor;
cout << factor << " ";
if ( n == 1 ) goto done;
cout << n;
done:
cout << endl;
【讨论】：
我不认为我会给这个特定版本打勾，表明革命性的性能改进，但我确实给了它额外的努力加一票，并希望其他人也这样做！ @Michael：更新了微优化。现在快了 25%……我将剩余可能收益的上限设定在 10% 左右。 @Michael：再次更新。我低估了剩余的改进（我猜是正面的恶作剧）......我又减少了 3 秒，这比你的程序减少了 16% 的复合时间或 38% 的时间。也就是说，速度接近两倍。 PGO 也让我一次降至 14.5 秒，但我似乎无法重现它。 步幅为 42 (2*3*7)，从 12 个候选者中产生多达 10 个候选者，这很有可能会稍微快一些，但我将把它作为练习留给读者。 （提示：它只能过滤掉最多 37 个素数。） @Michael：再次更新了更多素数，将运行时间再缩短 1.5 秒。我没有通过分析器运行它，但我怀疑通过平衡整数和向量 ALU 之间的 OR 运算可能会获得进一步的收益。【参考方案3】：Fermat's factorization method 可以简单快捷地找到成对的大素因数，只要在它走得太远和变慢之前停止它。然而，在我对随机数的测试中，这种情况很少见，看不到任何改进。
...不使用概率方法（例如 Miller-Rabin）来确定素数
在均匀分布的情况下，您 75% 的输入将需要十亿次循环迭代，因此即使您得到一个不确定的答案并不得不返回到试验划分，也值得首先在不太确定的技术上花费一百万次操作。
我发现 Pollard 的 Rho 方法的 Brent 变体非常好，但编码和理解起来更复杂。我见过的最好的例子来自forum discussion。该方法依赖于运气，但通常很有帮助，值得。
Miller-Rabin 素数检验实际上是高达 10^15 左右的确定性检验，可以省去搜索无果的麻烦。
我尝试了几十种变体，并选择了以下方法来分解 int64 值：
-
小因素的审判分工。 （我使用前 8000 个预先计算的素数。）
使用 Pollard 的 Rho 进行 10 次尝试，每次使用 16 次迭代
试除以 sqrt(n)。
请注意，Pollard 的 Rho 会找到不一定是素数的因子，因此可以使用递归来分解这些因子。
【讨论】：
【参考方案4】：融合了 Omnifarious 的一些想法，以及其他改进：
// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef unsigned __int64 ulong;
typedef std::vector<ulong> ULongVector;
// Caller needs to pass in an empty factors vector
void GetFactors(ULongVector &factors, ulong num)
if (num<2)
return;
ulong workingNum=num;
// Factor out factors of 2
while (workingNum%2==0)
factors.push_back(2);
workingNum/=2;
// Factor out factors of 3
while (workingNum%3==0)
factors.push_back(3);
workingNum/=3;
if (workingNum==1)
return;
// Factor out factors >=5
ulong nextOffset=2;
char nextShift = 1;
ulong n = 5;
ulong nn = 25;
do
// Is workingNum divisible by n?
if (workingNum%n==0)
// n is a factor!
// so is the number multiplied by n to get workingNum
// Insert n into the list of factors
factors.push_back(n);
// Divide working number by n
workingNum/=n;
// Test for done...
if (workingNum==1)
return;
// Try n again
else
nn += (n << (nextShift+1)) + (1<<(nextShift*2)); // (n+b)^2 = n^2 + 2*n*b + b*2
n += nextOffset;
nextOffset ^= 6;
nextShift ^= 3;
// invariant: nn == n*n
if (n & 0x100000000LL) break; // careful of integer wraparound in n^2
while (nn <= workingNum);
// workingNum is prime, add it to the list of factors
factors.push_back(workingNum);
【讨论】：
您需要交换n += 和 nn += 行——首先增加 n 会导致 nn 值有点过大，因此您可能会过早结束循环。试试你的程序因式分解 121 (11^2) 看看。
@Chris：你说得对，很好。也可以通过更改最后一项的符号 (n^2 = (n-b)^2 + 2*n*b - b^2) 来解决它，但是这个顺序的好处是它会更好地流水线，因为它不是t 取决于刚刚计算的 n 值。
@PigBen：是的，我刚刚发现了这个问题，缺少括号。你会再试一次吗？
已删除反对票。重新测试速度，看看它是否值得投票。
@Ben Voigt - 不幸的是，根据我的测试，你的速度比我的慢，而且花费的时间几乎是 OP 的两倍。我的需要大约 10% 的时间。 sqrt 是现代 CPU 上的快速单指令操作确实使得删除它并没有那么大的改进。【参考方案5】：
自然的概括是使用比 2 和 3 更多的已知素数来预先计算跳页。比如 2、3、5、7、11，模式周期为 2310（嗯，不错的数字）。也许更多，但它的收益递减 - 运行时间图可以准确确定预计算开始产生负面影响的确切位置，但当然它取决于要考虑的数字数量......
哈，我把编码细节留给你们。 :-)
干杯，
--阿尔夫
【讨论】：
收益递减点可能由 icache 决定，这将是一个比 2000 小很多的展开因子。 @Ben：请注意，模式周期为 6（对于已知的素数 2 和 3）只有 2 个跳过数，即 2 和 4，因此预计算后保留的数据要少得多比期间。我认为唯一要做的就是衡量。很难预测这样的事情！ :-) 我也有同样的想法，但是对于您添加的每组跳过，效率提升明显更少（我怀疑是 n^-x）。所以我认为收益递减点会公平地到来。加 5 甚至加 7 可能会侥幸逃脱，但我敢打赌，那之后就不值得了。但是，正如您所指出的，唯一的方法是测量。 我尝试了很多想法，但它们只会减慢功能。我认为可以将文字（2 或 4）的选择预编译到代码中，而无需使用任何非寄存器或文字内存。查找表并非如此，它会导致缓存命中 - 比嵌入实际指令的寄存器或文字要慢得多。【参考方案6】：这段代码比较慢，我很确定我明白为什么。它并不是慢得令人难以置信，但在 10-20% 的范围内肯定会慢一些。除法不应该在每个循环中执行一次，但唯一的方法是调用sqrt 或类似的方法。
// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef std::vector<ulong> ULongVector;
void GetFactors(ULongVector &factors, ulong num)
if (num<2)
return;
ulong workingNum=num;
ulong nextOffset=2;
while (workingNum%2==0)
factors.push_back(2);
workingNum/=2;
while (workingNum%3==0)
factors.push_back(3);
workingNum/=3;
ulong n = 5;
while ((workingNum != 1) && ((workingNum / n) >= n))
// Is workingNum divisible by n?
if (workingNum%n==0)
// n is a factor!
// so is the number multiplied by n to get workingNum
// Insert n into the list of factors
factors.push_back(n);
// Divide working number by n
workingNum/=n;
else
n+=nextOffset;
nextOffset=(nextOffset==2UL ? 4UL : 2UL);
if (workingNum != 1)
// workingNum is prime, add it to the list of factors
factors.push_back(workingNum);
【讨论】：
代码在 2 和 4 之间交替以跳过可被 3 整除的数字，因为 3 的所有因数都已被分解。感谢你的付出。我将在今天晚些时候尝试您的方法，并将其与原始方法进行比较，以查看通过代码重构是否会在速度方面有所提高... 消除sqrt ==好。移动 workingNum == 1 并划分到循环条件 == bad。
@Ben Voigt - 不幸的是，我不认为消除 sqrt 真的有帮助。在现代处理器上，这可以在一条指令中相当快地完成。
@Steve Jessop：优化编译器可能足够聪明，可以自己解决这个问题。
你们太沉迷于平方根计算的作用了。但是，正如已经提到的，每个因素只进行一次，最多可以有 63 个因素，而不是数十亿个潜在的 mod 测试。所以，请专注于代码“热”部分的瓶颈。【参考方案7】：
我不确定这些效果如何，但不是
while (workingNum%2==0)
你可以这样做
while (workingNum & 1 == 0)
我不确定 gcc 或 msvc（或您正在使用的任何编译器）是否足够聪明以更改 workingNum%2 表达式，但很可能它正在执行除法并查看 edx 中的模数...
根据您的编译器，我的下一个建议可能完全没有必要，但您可以尝试将 workingNum /= 3 放在方法调用之前。 G++ 可能足够聪明，可以看到不必要的除法，只需使用 eax 中的商（您也可以在更大的循环中执行此操作）。或者，更彻底（但很痛苦）的方法是内联汇编以下代码。
while (workingNum%3==0)
factors.push_back(3);
workingNum/=3;
编译器可能将模运算转换为除法，然后查看 edx 中的模数。问题是，您正在再次执行除法（我怀疑编译器是否看到您只是在循环条件下隐式执行了除法）。所以，你可以内联组装它。这提出了两个问题：
1) push_back(3) 的方法调用。这可能会弄乱寄存器，因此完全没有必要。
2) 获取 workingNum 的寄存器，但这可以通过进行初始模块检查（强制它进入 %eax）来确定，或者在当前时刻，它将/应该在 eax 中。
您可以将循环编写为（假设 workingNum 在 eax 中，这是 32 位 AT&T 语法，只是因为我不知道 64 位汇编或 Intel 语法）
asm( "
movl $3, %ebx
WorkNumMod3Loop:
movl %eax, %ecx # just to be safe, backup workingNUm
movl $0, %edx # zero out edx
idivl $3 # divide by 3. quotient in eax, remainder in edx
cmpl $0, %edx # compare if it's 0
jne AfterNumMod3Loop # if 0 is the remainder, jump out
# no need to perform division because new workingNum is already in eax
#factors.push_back(3) call
je WorkNumMod3Loop
AfterNumMod3Loop:
movl %ecx, %eax"
);
您应该查看这些循环的程序集输出。您的编译器可能已经在进行这些优化，但我对此表示怀疑。如果在某些情况下将 workingNum /= n 放在方法调用之前可以稍微提高性能，我不会感到惊讶。
【讨论】：
以上是关于蛮力，单线程素数分解的主要内容，如果未能解决你的问题，请参考以下文章