ライフハック:1秒あたり1ギガバイトのdoubleを解析する方法





C ++コードの文字列からdouble値を読み取る方法は?



std::stringstream in(mystring);
while(in >> x) {
   sum += x;
}
      
      





GCC8.3コンパイラを搭載したIntelSkylakeでは、このコードは50 MB / sを解析します。ハードドライブは数GB /秒の速度でシーケンシャル読み取りを簡単に提供できるため、ディスクからの読み取り速度ではなく、解析速度が制限されることは間違いありません。どうすればスピードアップできますか?



それ自体を示唆する最初のことは、C ++のストリームによって提供される便利さを放棄し、strtod(3)を直接呼び出すことです。



do {
    number = strtod(s, &end);
    if(end == s) break;
    sum += number;
    s = end; 
} while (s < theend);
      
      





速度は最大90MB /秒になります。プロファイリングは、ストリームから読み取る場合、読み取り番号ごとに約1600の命令が実行され、strtodを使用する場合、番号ごとに約1100の命令が実行されることを示しています。CおよびC ++標準ライブラリは、普遍性と移植性の要件によって正当化できます。ただし、解析をdoubleとx64のみに制限すると、はるかに効率的なコードを記述できます。数値ごとに280命令で十分です。



整数の解析



サブタスクから始めましょう。8桁の10進数が与えられた場合、それをintに変換する必要があります。学校では、私たち全員がこれをサイクルで行うように教えられました。



int sum = 0;
for (int i = 0; i <= 7; i++)
{
	sum = (sum * 10) + (num[i] - '0');
}
return sum;
      
      





コンパイルされたGCC9.3.1 -O3、私にとってこのコードは3 GB / sを処理します。それをスピードアップする明白な方法は、ループを逆にすることです。しかし、コンパイラはそれ自体を行います。



8桁の数値の10進表記は、int64_t変数に収まることに注意してください。たとえば、文字列「87654321」は、int64_t値0x3132333435363738と同じ方法で格納されます(最初のバイトには最下位8ビットが含まれ、最後のバイトには-最も重要なもの)。これは、8回のメモリアクセスの代わりに、シフトを使用して各桁の値を割り当てることができることを意味します。



int64_t sum = *(int64_t*)num;
return (sum      & 15) * 10000000 +
    ((sum >> 8)  & 15) * 1000000 +
    ((sum >> 16) & 15) * 100000 +
    ((sum >> 24) & 15) * 10000 +
    ((sum >> 32) & 15) * 1000 +
    ((sum >> 40) & 15) * 100 +
    ((sum >> 48) & 15) * 10 +
    ((sum >> 56) & 15);
      
      





まだ加速はありません。それどころか、速度は1.7 GB / sに低下します!さらに進んでみましょう:そしてマスク0x000F000F000F000Fを使用すると、0x0002000400060008-偶数の位置に10進数が表示されます。それらのそれぞれを以下と組み合わせてみましょう:



sum = ((sum & 0x000F000F000F000F) * 10) + 
      ((sum & 0x0F000F000F000F00) >> 8);
      
      





その後、0x3132333435363738は0x00(12)00(34)00(56)00(78)に変わります-偶数の位置のバイトはゼロになり、奇数の位置には2桁の10進数のバイナリ表現が含まれます。



return (sum        & 255) * 1000000 +
      ((sum >> 16) & 255) * 10000 +
      ((sum >> 32) & 255) * 100 +
      ((sum >> 48) & 255);
      
      



-8桁の数値の変換を完了します。



同じトリックを2桁の数字のペアで繰り返すことができます。

sum = ((sum & 0x0000007F0000007F) * 100) +
      ((sum >> 16) & 0x0000007F0000007F);
      
      



-0x00(12)00(34)00(56)00(78)は0x0000(1234)0000(5678)になります。

結果として得られる4桁のペア:

return ((sum & 0x3FFF) * 10000) + ((sum >> 32) & 0x3FFF);
      
      



-0x00(12)00(34)00(56)00(78)は0x00000000(12345678)になります。結果のint64_tの小さい方の半分が望ましい結果です。コードが著しく単純化されているにもかかわらず(7回ではなく3回の乗算)、解析速度(2.6 GB / s)は単純な実装よりも劣っています。



ただし、次のアクションでマスク0x007F007F007F007Fが適用されることに気付いた場合でも、数値のペアの組み合わせは簡略化できます。つまり、null許容バイトにゴミが残る可能性があります。



sum = ((sum & 0x0?0F0?0F0?0F0?0F) * 10) + ((sum & 0x0F??0F??0F??0F??) >> 8) =
   = (((sum & 0x0F0F0F0F0F0F0F0F) * 2560)+ (sum & 0x0F0F0F0F0F0F0F0F)) >> 8 =
    = ((sum & 0x0F0F0F0F0F0F0F0F) * 2561) >> 8;
      
      





簡単に言えば、2つではなく1つのマスクであり、追加はありません。残りの2つの式は、同じ方法で簡略化されます。



sum = ((sum & 0x00FF00FF00FF00FF) * 6553601) >> 16;
return((sum & 0x0000FFFF0000FFFF) * 42949672960001) >> 32;
      
      





このトリックはレジスター内のSIMD(SWAR)と呼ばれます。これを使用すると、解析速度は3.6 GB / sに増加します。



同様のSWARトリックを使用して、8文字の文字列が8桁の10進数であるかどうかを確認できます。



return ((val & 0xF0F0F0F0F0F0F0F0) |
      (((val + 0x0606060606060606) & 0xF0F0F0F0F0F0F0F0) >> 4))
            == 0x3333333333333333;
      
      





ダブルデバイス



IEEE標準では、これらの数値内の仮数に52ビット、指数に11ビットが割り当てられています。これは、番号が次のように保存されることを意味します ±1.m2e どこ 0m<252<1016 そして 1022e+1023 ..。特に、これは、doubleの10進表記では、最上位の17桁のみが重要であることを意味します。最下位ビットは、double値に影響を与えることはできません。有効数字17桁でさえ、実際のアプリケーションに必要とされるよりもはるかに多くなります。非正規化数は、作業を少し複雑にします







21074 21022 それに対応して仮数の有効ビット数が少なくなります)および丸め要件(10進数は最も近い2進数で表す必要があり、数値が最も近い2進数のちょうど中間にある場合は、偶数の仮数が優先されます)。これにより、コンピューターAが数値Xを有効数字17桁の10進文字列に変換した場合、この文字列を読み取るコンピューターBは、AとBが同じであるかどうかに関係なく、ビットごとに同じ数値Xを受け取ることが保証されます。プロセッサモデル、オペレーティングシステム、およびプログラミング言語。 (丸め誤差がコンピューターごとに異なる場合、コードをデバッグするのがどれほど楽しいか想像してみてください。)丸め要件のために、最近言及された誤解が生じます。 Habréの場合:小数0.1は、それに最も近い2進分数として表されます。 7205759403792794256 正確に0.1000000000000000055511151231257827021181583404541015625に等しい; 0.2-形式で 7205759403792794255 正確に0.200000000000000011102230246251565404236316680908203125に等しい; しかし、それらの合計は0.3に最も近い2進分数ではありません:下からの近似 5404319552844595254 = 0.299999999999999988897769753748434595763683319091796875の方が正確であることがわかります。したがって、IEEE標準で 、0.3以外の結果を生成するために0.1 +0.2を追加する必要があります。



ダブル解析



整数アルゴリズムの単純な一般化は、10.0による反復的な乗算と除算で構成されます。整数の解析とは異なり、ここでは、上位桁の丸め誤差が下位桁の丸め誤差を「隠す」ように、桁を下位から上位に処理する必要があります。同時に、仮数の構文解析は整数の構文解析に簡単に縮小できます。正規化を変更して、仮数の2進点が仮数の先頭ではなく末尾になるようにするだけで十分です。結果の53ビット整数は、必要な10の累乗で乗算または除算されます。最後に、指数から52を引きます。ポイントを52ビット左に移動します-IEEE標準に準拠している必要があります。さらに、注意すべき3つの重要な事実があります。



  1. 5で乗算または除算する方法を学ぶだけで十分であり、2で乗算または除算することは、指数の増分または減分にすぎません。
  2. uint64_t 5 0xcccccccccccccccd 66 , , 0xcccccccccccccccd26615=15266<268 64 ( );
  3. 10324<21074 21024<10309 ; , 309 , 324 0xcccccccccccccccd, . ( 53 ; 128- , 53- 53- .) 633 double ( , ⅕, 7205759403792794 – 0xcccccccccccccccd, 53 ), double ; 53 , . , , 64 64 , , 128- . .


標準の丸めの複雑さは、結果が最も近い2進小数のちょうど中間にあることを確認するために、結果の54の有効ビットだけでなく、最下位のゼロビットは「すべてが正常であること」を意味します。 1つは「真ん中を打つ」ことを意味します)が、最下位ビットの後にゼロ以外の継続があったかどうかを監視する必要があります。特に、これは、数値の10進表記で最上位17桁だけを考慮するだけでは不十分であることを意味します。これらは、±1の精度でのみバイナリ仮数を定義し、目的の丸め方向を選択するには、考慮する必要があります。下の桁。たとえば、10000000000000003と10000000000000005は同じdouble値(10000000000000004に等しい)であり、10000000000000005.00(多くのゼロ)001は異なる値(10000000000000006に等しい)です。



このパーサーを発明したケベック通信大学(TÉLUQ)のDaniel Lemire教授は、githubで公開しました このライブラリは2つの関数を提供します。1 from_chars



つは文字列を解析してdoubleにし、もう1つはfloatにします。彼の実験は、99.8%のケースで19の有効な10進数(64ビット)を考慮するだけで十分であることを示しました:64ビット仮数の2つの連続した値に対して同じ最終的なdouble値が得られた場合、これは正しい値です。0.2%の場合にのみ、これら2つの値が一致しない場合、常に正しい丸めを実装する、より複雑でより普遍的なコードが実行されます。






All Articles