私の他の答えはかなり長く、これとはかなり異なるので、ここに別のものがあります.
最初の 2 つの素数の倍数を単に除外したり、関連するすべての素数をそれぞれ 1 バイトにエンコードしたりするのではなく、このプログラムは、8 ビット (具体的には 2 から 211) に収まるすべての素数の倍数を除外します。したがって、33% を渡す代わりに、数値の場合、これは約 10% を除算演算子に渡します。
素数は 3 つの SSE レジスタに保持され、ランニング カウンターを使用したモジュラスは別の 3 つのレジスタに保持されます。カウンターの素数のモジュラスがゼロに等しい場合、カウンターは素数にはなりません。また、いずれかの法が 1 に等しい場合、(counter+2) は (counter+30) まで素数にすることはできません。偶数は無視され、+3、+6、+5 などのオフセットはスキップされます。ベクトル処理では、16 バイト サイズの変数を一度に更新できます。
マイクロ最適化を完全に実行した後 (ただし、インライン ディレクティブ以外にプラットフォーム固有のものは何もありません)、パフォーマンスが 1.78 倍向上しました (ラップトップで 13.4 秒に対して 24 秒)。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 の間の素数は、1 つのストライドで複数の倍数を持つ可能性があるため、ベクター内のさまざまな位相に複数のコピーを配置します。
- 因数分解し
remove_factorます。
- 条件付きで予測不可能な
remove_factor呼び出しを非分岐配列書き込みに変更します。
- 呼び出しで最終ループを完全に展開
remove_factorし、関数が常にインライン化されていることを確認します。
- 候補の中には常に 7 の倍数が存在するため、展開された最後の繰り返しを削除します。
- 十分に小さい残りのすべての素数を含む別のベクトルを追加します。
- カウンターにバイアスを追加してスペースを増やし、さらに別のベクトルを追加します。これで、最大 16 ビットに達することなくフィルタリングできる可能性のある素数はあと 6 つしかなく、レジスタも不足しています。ループには、3 つの素数ベクトル、3 つのモジュラス ベクトル、検索する 8 つの定数、およびそれぞれ 1 つの定数が必要です。インクリメントして範囲チェックを行います。すると16になります。
- このアプリケーションではゲインは最小限 (ただし正) ですが、この手法の本来の目的は、他の回答のふるいの素数をフィルター処理することでした。どうぞお楽しみに…</li>
読み取り可能なバージョン:
/*
* 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;
}