c++ - 蛮力，单线程素数分解

Question

需要考虑的是以下函数，该函数可用于(相对快速地)将 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);
}

谢谢

编辑:我添加了更多评论。通过引用传入 vector 的原因是允许在调用之间重用 vector 并避免动态分配。 vector 没有在函数中清空的原因是允许将当前“num's”因子附加到 vector 中已经存在的因子上的奇怪要求。

函数本身并不漂亮，可以重构，但问题是如何让算法更快。所以，请不要提出关于如何使函数更漂亮、更易读或 C++ish 的建议。那是儿戏。改进这个算法，使其能够更快地找到(证明)因子是困难的部分。

更新:Potatoswatter 到目前为止有一些出色的解决方案，请务必在底部附近查看他的 MMX 解决方案。

Answer 1

将这种方法与(预先生成的)筛子进行比较。取模很昂贵，所以这两种方法本质上都做两件事:生成潜在因子，并执行取模运算。任何一个程序都应该在比模数更少的周期内合理地生成一个新的候选因子，因此任何一个程序都是模数绑定(bind)的。

给定的方法过滤掉所有整数的恒定比例，即 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 快，则筛子更快。根据我的表格，Core2 上的 div 吞吐量最多为每 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