暴力破解、单线程素数分解

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);
}

感谢

编辑：我添加了更多的注释。将向量按引用传递的原因是为了允许在调用之间重复使用向量并避免动态分配。函数中向量不被清空的原因是为了满足将当前“num”的因子追加到向量中已有因子的奇特要求。

这个函数本身不太美观，可以进行重构，但问题是如何使算法更快。因此，请不要提出如何让函数更漂亮、可读或类似C++的建议。这是小儿科。提高此算法的速度以更快地找到（已证明的）因子才是难点。

更新：Potatoswatter 到目前为止提供了一些很好的解决方案，请确保查看他的 MMX 解决方案。

- Michael Goldshteyn

---

你的问题是什么？ - Paul R

在不影响效率的情况下，你可以将此代码更改为接受输出迭代器的模板，从而消除对vector的依赖。用户仍然可以使用vector，例如 GetFactors(n, back_inserter(v))。 - Jon Purdy

最快大约需要6秒钟才能分解最难的64位无符号整数（例如具有非常大的因子的质数和合数）。我想知道是否可以改进这个简单算法，使其更加快速。 - Michael Goldshteyn

@Michael：您可以尝试将循环展开一次，这样就不需要条件+2或+4了。相反，在循环的中间位置无条件地增加2，在末尾增加4。这可能没有任何影响，但可能会50-100％的时间预防分支错误预测。谁知道呢。 - Steve Jessop

@Michael：另外，Miller-Rabin有一种确定性变体，但我不确定它与你的代码相比有多快。此外，它依赖于广义黎曼猜想...这可能只是在你为技术论文编写此算法时才需要考虑的问题。顺便说一下，概率素性测试产生的错误对你来说不应该是主要关注的问题，假设你正确使用了这样的素性测试。请记住，通常情况下，您可以获得快速结果，并且仍然具有比PC故障引入的错误概率更低的概率。 - Brian

Answer 1

比较这种方法和（预先生成的）筛法。取模很耗费时间，因此这两种方法本质上都要做两件事：生成潜在因子并执行取模运算。任何一种程序都应该在比取模更少的周期内合理地生成新的候选因子，因此任何一种程序都是取模限制的。

所给出的方法过滤掉所有整数的一个恒定比例，即2和3的倍数，或75%。四分之一（如所给）的数字被用作取模运算符的参数。我称之为跳过滤器。

另一方面，筛法仅使用素数作为取模运算符的参数，连续素数之间的平均差异由素数定理控制为1/ln(N)。例如，e^20略小于5亿，因此超过5亿的数字有不到5%的几率是素数。如果考虑到所有小于2^32的数字，5%是一个好的经验法则。

因此，筛法在div操作上花费的时间将比您的跳过滤器少5倍。要考虑的下一个因素是筛法生成素数的速度，即从内存或磁盘中读取它们的速度。如果获取一个素数比4个div操作更快，那么筛法就更快。根据我的表格，在我的Core2上div吞吐量最多为每12个周期执行一次。这些将是困难的除法问题，因此让我们保守地预算每个素数50个周期。对于一个2.5 GHz的处理器，那是20纳秒。

在20纳秒内，一个50 MB/秒的硬盘可以读取大约1个字节。简单的解决方案是每个质数使用4个字节，因此驱动器会变慢。但是，我们可以更聪明一些。如果我们想按顺序编码所有的质数，我们只需要编码它们的差异。同样，预期的差异是1/ln(N)。而且，它们都是偶数，这可以节省一个额外的位。它们永远不会为零，这使得扩展到多字节编码是免费的。因此，使用一个字节来存储每个质数，在一个字节中可以存储高达512的差异，这使我们能够达到303371455241，根据那个维基百科文章。

因此，根据硬盘的不同，存储的质数列表应该与验证素性的速度大致相等。如果它可以存储在RAM中（它是203MB，所以后续运行可能会命中磁盘缓存），那么问题就完全消失了，因为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

- Potatoswatter

---

理论听起来很有趣，但最终还是要看它的表现。所以，当你修复了这个 bug，请在这里发布你的结果。我特别想知道它如何处理非常大的质数（例如18446744073709542713），以及如何处理仅由两个大质数因子组成的大数字（例如18446743721522234449）。 - Michael Goldshteyn

...而且这需要我25秒的时间。哇，你有一台快速的机器。 - Potatoswatter

所以结论是我们可以预料的：如果您偶尔想要分解一个大数，那么Mike的代码看起来更好，但如果您有兴趣分解许多数字并且磁盘空间不是问题，筛法会更快（对于20位数字，即两个10位质数的乘积，速度快五倍）。比较常数成本（生成质数文件）与一次因式分解的成本乘以调用次数。 - Wok

@wok: 差不多就是这样。然而，拥有表格的开销非常低。生成该表格需要75秒，不包括I/O但包括一些操作系统开销。我还没有尝试使用普通数组而不是streambuf_iterator来运行素数生成器，但如果您以这种方式运行它，速度应该会更快一些。现在我正在尝试优化，使其不会占用太多内存。平衡点已经只有3个质数需要验证，我可能会将其降至2个。 - Potatoswatter

@Michael：我不确定SO是否能够很好地处理同一人的多个答案，所以请确保您查看页面底部的我的新程序。它需要SSE2但不会在堆上放置任何东西。 - Potatoswatter

Answer 2

我的另一个答案比较长，与这个答案有很大不同，所以我来说点别的。

这个程序并不仅仅过滤掉第一、二个质数的倍数，也不是将所有相关的质数都编码成一个字节。相反，它过滤掉了适合在八位中存储的所有质数的倍数，具体来说就是2到211。因此，它只将大约10%的数字传递给除法运算符，而不是33%。

这些质数被保存在三个SSE寄存器中，它们与运行计数器的模数被保存在另外三个寄存器中。如果任何一个质数与计数器的模数等于零，则计数器不是质数。而且，如果任何模数等于1，那么(counter+2)就不可能是质数，以此类推，直到(counter+30)。偶数被忽略，像+3、+6和+5这样的偏移量也被跳过。向量处理允许同时更新16个字节大小的变量。

经过对其进行了各种微小的优化（但没有比内联指令更具体的平台），我获得了1.78倍的性能提升（在我的笔记本电脑上，24秒与13.4秒相比）。如果使用大数库（即使是非常快的库），优势会更大。请参见下面一个更易读的、进行了预优化的版本。

/*
 *  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调用改为非分支数组写入。
• 完全展开最终循环，并确保函数始终内联。
  • 消除最后一个展开迭代，因为候选项中总是有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;
}

- Potatoswatter

---

我认为除非有革命性的性能提升，否则我不会给这个特定版本打勾，但我确实为额外的努力投了一票，并希望其他人也这样做！ - Michael Goldshteyn

@Michael：使用微小优化进行了更新，现在快了25%…我认为剩余可能的收益上限大约在10%的复合增长率左右。 - Potatoswatter

@Michael：我又更新了一下。我低估了剩余的改进（可能是积极诅咒）……我又节省了3秒，这相当于16％的复合效果或比你的程序少38％的时间。也就是说，接近两倍的速度。此外，PGO让我只用了14.5秒，但我好像无法重现这个结果。 - Potatoswatter

很有可能步长为42（237）会更快，从12个数字中产生最多10个候选数，但我将把这留给读者作为练习。（提示：它只能过滤掉37以下的质数。） - Potatoswatter

@Michael：用更多的质数再次更新了程序，运行时间又减少了1.5秒。我还没有通过分析器运行它，但我怀疑通过平衡整数和向量ALU之间的OR操作，可以获得进一步的收益。 - Potatoswatter

Answer 3

费马分解法对于找到大质因数的一对来说是简单快速的，只要在过程变慢之前停止它。然而，在我对随机数的测试中，这种情况太少了，以至于看不到任何改进。

...没有使用概率方法（例如Miller-Rabin）来确定素性

使用均匀分布，你的75%输入将需要十亿次循环迭代，因此值得先花费一百万次操作在不那么确定的技术上，即使你得到一个无法确定的答案，仍需退回到试除法。

我发现Pollard's Rho算法的Brent变体非常好，尽管编码和理解更加复杂。我见过的最好的例子来自于这个论坛讨论。该方法依赖于运气，但足够常见，值得一试。

Miller-Rabin素性检验实际上是确定性的，可以检验到10^15左右，这可以节省你进行无用搜索的麻烦。

我尝试了几十种变体，并选择了以下方法来分解int64值：

1. 使用前8000个预计算的质数进行小因子试除。
2. 使用Pollard's Rho进行10次尝试，每次使用16次迭代。
3. 使用试除法直到sqrt(n)。

请注意，Pollard's Rho可以找到不一定是质数的因子，因此可以使用递归来分解它们。

- xan

Answer 4

这段代码速度较慢，我相信我知道原因。它不是非常慢，但肯定比正常速度慢了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);
  }
}

- Omnifarious

---

代码在2和4之间交替，以跳过可被3整除的数，因为所有3的因数已经被分解出来。感谢您的努力。我会在今天晚些时候尝试您的方法，并与原始方法进行时间比较，以查看是否存在速度上的渐进改进。 - Michael Goldshteyn

消除 sqrt == 好。将 workingNum == 1 和除法移入循环条件 == 不好。 - Ben Voigt

@Ben Voigt - 很不幸，我认为去除sqrt并不能真正解决问题。现代处理器可以在单个指令中快速执行此操作。 - Omnifarious

@Steve Jessop：一个优化编译器可能足够聪明，能够自行找出这一点。 - Greg Hewgill

你们过于关注平方根计算的作用。但是，如已经提到的，这仅在每个因子上进行一次，最多可能存在63个因子与数十亿个可能的模数测试相比。所以，请专注于“热”代码部分的瓶颈。 - Michael Goldshteyn

Answer 5

结合 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);
}

- Ben Voigt

---

你需要交换 n += 和 nn += 这两行代码 -- 先增加 n 会导致 nn 的值稍微有点大，这可能会使循环结束得太早。试着运行你的程序来分解 121（11^2）看看。 - Chris Dodd

@Chris：你说得对，发现得好。也可以通过改变最后一项的符号来修复它（n^2 = (n-b)^2 + 2nb - b^2），但这个顺序的好处是它不依赖于刚刚计算出的n的值，因此可以更好地进行流水线处理。 - Ben Voigt

@PigBen：是的，我刚发现问题出在括号缺失。你能再试一次吗？ - Ben Voigt

取消了踩，重新测试速度，看是否应该点赞。 - Benjamin Lindley

@Ben Voigt - 不幸的是，根据我的测试，你的速度比我的慢，并且几乎比 OP 的慢了近一倍。 我的算法大约慢10％。 在现代CPU上，sqrt 是一种快速的单指令操作，因此删除它并不能带来太大的改进。 - Omnifarious

奇怪。它与核弹“nn”和“nextShift”相比如何？在while条件中只需执行“n * n”？我猜术语“(1<<nextShift*2)”实际上不需要在运行时计算，也可以使用异或在两个值（4、16）之间切换。此外，还必须查看编译器为环绕测试编写的代码。 - Ben Voigt

Answer 6

我不确定这些会有多大的效果，但是可以尝试一下：

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）获取一个工作寄存器，但可以通过进行初始模块检查（强制将其放入%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 放置在前面可能会稍微提高性能，这一点也不奇怪。

- Kizaru

Answer 7

自然的推广是使用更多的已知质数来预计算跳过值，而不仅仅是2和3。比如2、3、5、7、11等，模式周期为2310（哦，好数字）。也许还有更多，但收益递减--运行时间的图表可以确定预计算开始产生负面影响的确切位置，但当然这取决于要分解的数字数量...

哈哈，编码细节就留给你们了。:-)

祝福&帮助。

- Alf

- Cheers and hth. - Alf

---

退化收益点可能由icache决定，其展开因子要小得多，可能不到2000。 - Ben Voigt

注意，对于已知质数2和3的模式周期为6，只有2个跳过数字，即2和4，因此预计算后保留的数据要少得多。我认为唯一能做的就是测量。这种事情很难预测！ :-) - Cheers and hth. - Alf

我也有同样的想法，但是每增加一组跳跃，效率提高的幅度显著降低（我猜测是n^-x），所以我认为边际收益点会相当快地到来。你可能可以添加5个或者7个，但我敢打赌在此之后它就不值得了。不过，正如你所指出的，唯一的方法就是进行衡量。 - Omnifarious

我尝试了很多想法，但它们只会减慢函数的速度。我认为字面值（2或4）的选择可以预编译到代码中，而不需要使用非寄存器或字面内存。这对于查找表来说并不正确，因为它会产生缓存命中 - 明显比寄存器或字面值嵌入实际指令要慢得多。 - Michael Goldshteyn