可以考虑使用以下函数,该函数可以(相对快速)将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); }
谢谢
编辑:我添加了更多的评论。通过引用传递向量的原因是允许向量在调用之间重用,并避免动态分配。向量未在函数中清空的原因是允许将当前“数字”因数附加到向量中已存在的因数的奇怪要求。
函数本身并不漂亮,可以重构,但是问题在于如何使算法更快。因此,请不要提出任何有关如何使函数更美观,可读性更好或C ++的建议。那是孩子的戏。改进此算法,使其能够更快地找到(证明)因素是困难的部分。
更新:到目前为止, Potawswattswatter 拥有一些出色的解决方案,请确保也从底部检查他的MMX解决方案。
将这种方法与(预生成的)筛网进行比较。模数昂贵,因此这两种方法本质上都做两件事:生成潜在因子并执行模数运算。任一程序都应在比模数更少的周期内合理地生成新的候选因子,因此任一程序都是模数界。
给定的方法过滤掉所有整数的恒定比例,即2和3的倍数,即75%。四分之一(给定的数字)用作取模运算符的参数。我将其称为跳过过滤器。
另一方面,筛子仅使用质数作为模运算符的参数,并且连续质数之间的平均差由质数定理控制为1 / ln(N)。例如, e ^ 20不到5亿,因此超过5亿的数字有5%的机会成为质数。如果考虑不超过2 ^ 32的所有数字,则5%是一个很好的经验法则。
因此,筛网的div操作时间比跳过过滤器少5倍。下一个要考虑的因素是筛子产生质子的速度,即从内存或磁盘读取质子的速度。如果获取一个质数的速度快于4 divs,则筛分速度会更快。根据我的表div,Core2上的吞吐量最多为每12个周期1个。这些将是困难的除法问题,因此让我们保守地估计每个素数50个周期。对于2.5 GHz处理器,则为20纳秒。
div
在20 ns中,一个50 MB /秒的硬盘驱动器可以读取大约一个字节。简单的解决方案是每个素数使用4个字节,因此驱动器将变慢。但是,我们可以变得更加聪明。如果我们要按顺序对所有素数进行编码,则可以对它们的差异进行编码。同样,期望的差是1 / ln(N)。而且,它们都是偶数,这节省了额外的费用。而且它们永远不会为零,这使得扩展到多字节编码成为免费的。因此,每个素数使用一个字节,一个字节中最多可以存储512个差异,根据Wikipedia文章,这最多可以使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