世界最速の atof 実装はどこで見つけられますか?
https://stackoverflow.com/questions/98586
-
01-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
US-en ロケール、ASCII、非科学表記用に最適化された IA32 上の非常に高速な atof() 実装を探しています。Windows のマルチスレッド CRT は、isdigital() を呼び出すたびにロケールの変更をチェックするため、ここで悲惨な状態に陥ります。現在の最高のパフォーマンスは、perl + tcl の atof 実装の最高の結果から得られており、msvcrt.dll の atof を一桁上回っています。もっと改善したいのですが、アイデアがありません。BCD 関連の x86 命令は有望に思えましたが、perl/tcl C コードを上回るパフォーマンスを実現することはできませんでした。SO'ers は、最高のものへのリンクを見つけることができますか?非 x86 アセンブリベースのソリューションも歓迎します。
最初の回答に基づいた説明:
このアプリケーションでは、最大 2 ulp の誤差は問題ありません。
変換される数値は、小さなバッチでネットワーク経由で ASCII メッセージとして到着するため、アプリケーションはそれらを可能な限り短い遅延で変換する必要があります。
解決
精度の要件は何ですか?本当に「正しい」（指定された 10 進数に最も近い浮動小数点値を常に取得する）必要がある場合は、標準ライブラリのバージョンに勝つのはおそらく難しいでしょう（ロケールのサポートを削除することを除いて、すでに行っています）。これには、任意精度の演算を行う必要があります。1、2 ulp のエラー (および非正規数の場合はそれ以上) を許容できる場合は、cruzer が提案した種類のアプローチが機能し、より高速になる可能性がありますが、0.5ulp 未満の出力は絶対に生成されません。整数部分と小数部分を別々に計算し、最後に分数を計算すると、精度が向上します。12345.6789 の場合、6*.1 + 7*.01 + 8*.001 + 9*0.0001) ではなく、12345 + 6789 / 10000.0 として計算します。0.1 は無理数の 2 進分数であり、0.1^ を計算すると誤差が急速に蓄積するためです。 n.これにより、浮動小数点ではなく整数を使用してほとんどの計算を実行できるようになります。
BCD 命令は (IIRC) 286 以降ハードウェアに実装されておらず、現在では単純にマイクロコード化されています。特に高性能であるとは考えられません。
他のヒント
私がコーディングを終えたこの実装は、デスクトップに組み込まれている「atof」の 2 倍の速度で実行されます。1024*1024*39 の数値入力を 2 秒で変換しますが、私のシステムの標準 GNU 'atof' では 4 秒かかります。(セットアップ時間やメモリの取得などすべてを含みます)。
アップデート：申し訳ありませんが、2 倍の速度の申請を取り消さなければなりません。変換するものがすでに文字列内にある場合は高速ですが、ハードコードされた文字列リテラルを渡す場合は、atof とほぼ同じになります。ただし、ragel ファイルとステート マシンを少し調整すると、特定の目的でより高速なコードを生成できる可能性があるため、ここではそのままにしておきます。
https://github.com/matiu2/yajp
興味深いファイルは次のとおりです。
https://github.com/matiu2/yajp/blob/master/tests/test_number.cpp
https://github.com/matiu2/yajp/blob/master/number.hpp
また、変換を行うステート マシンにも興味があるかもしれません。
役に立つかもしれないものを実装しました。atof と比較すると約 5 倍高速になり、atof と併用すると __forceinline 約10倍速くなります。もう 1 つの優れた点は、crt 実装とまったく同じ演算を行うように作られていることです。もちろん、いくつかの短所もあります。
- 単精度浮動小数点のみをサポートします。
- #INF などの特別な値はスキャンしません。
__forceinline bool float_scan(const wchar_t* wcs, float* val)
{
int hdr=0;
while (wcs[hdr]==L' ')
hdr++;
int cur=hdr;
bool negative=false;
bool has_sign=false;
if (wcs[cur]==L'+' || wcs[cur]==L'-')
{
if (wcs[cur]==L'-')
negative=true;
has_sign=true;
cur++;
}
else
has_sign=false;
int quot_digs=0;
int frac_digs=0;
bool full=false;
wchar_t period=0;
int binexp=0;
int decexp=0;
unsigned long value=0;
while (wcs[cur]>=L'0' && wcs[cur]<=L'9')
{
if (!full)
{
if (value>=0x19999999 && wcs[cur]-L'0'>5 || value>0x19999999)
{
full=true;
decexp++;
}
else
value=value*10+wcs[cur]-L'0';
}
else
decexp++;
quot_digs++;
cur++;
}
if (wcs[cur]==L'.' || wcs[cur]==L',')
{
period=wcs[cur];
cur++;
while (wcs[cur]>=L'0' && wcs[cur]<=L'9')
{
if (!full)
{
if (value>=0x19999999 && wcs[cur]-L'0'>5 || value>0x19999999)
full=true;
else
{
decexp--;
value=value*10+wcs[cur]-L'0';
}
}
frac_digs++;
cur++;
}
}
if (!quot_digs && !frac_digs)
return false;
wchar_t exp_char=0;
int decexp2=0; // explicit exponent
bool exp_negative=false;
bool has_expsign=false;
int exp_digs=0;
// even if value is 0, we still need to eat exponent chars
if (wcs[cur]==L'e' || wcs[cur]==L'E')
{
exp_char=wcs[cur];
cur++;
if (wcs[cur]==L'+' || wcs[cur]==L'-')
{
has_expsign=true;
if (wcs[cur]=='-')
exp_negative=true;
cur++;
}
while (wcs[cur]>=L'0' && wcs[cur]<=L'9')
{
if (decexp2>=0x19999999)
return false;
decexp2=10*decexp2+wcs[cur]-L'0';
exp_digs++;
cur++;
}
if (exp_negative)
decexp-=decexp2;
else
decexp+=decexp2;
}
// end of wcs scan, cur contains value's tail
if (value)
{
while (value<=0x19999999)
{
decexp--;
value=value*10;
}
if (decexp)
{
// ensure 1bit space for mul by something lower than 2.0
if (value&0x80000000)
{
value>>=1;
binexp++;
}
if (decexp>308 || decexp<-307)
return false;
// convert exp from 10 to 2 (using FPU)
int E;
double v=pow(10.0,decexp);
double m=frexp(v,&E);
m=2.0*m;
E--;
value=(unsigned long)floor(value*m);
binexp+=E;
}
binexp+=23; // rebase exponent to 23bits of mantisa
// so the value is: +/- VALUE * pow(2,BINEXP);
// (normalize manthisa to 24bits, update exponent)
while (value&0xFE000000)
{
value>>=1;
binexp++;
}
if (value&0x01000000)
{
if (value&1)
value++;
value>>=1;
binexp++;
if (value&0x01000000)
{
value>>=1;
binexp++;
}
}
while (!(value&0x00800000))
{
value<<=1;
binexp--;
}
if (binexp<-127)
{
// underflow
value=0;
binexp=-127;
}
else
if (binexp>128)
return false;
//exclude "implicit 1"
value&=0x007FFFFF;
// encode exponent
unsigned long exponent=(binexp+127)<<23;
value |= exponent;
}
// encode sign
unsigned long sign=negative<<31;
value |= sign;
if (val)
{
*(unsigned long*)val=value;
}
return true;
}
各状態が N 番目の入力桁または指数桁を処理するステート マシンに相当するものを (手動で) 構築したいようです。このステート マシンはツリーのような形になります (ループはありません!)。目標は、可能な限り整数演算を行うことと、(当然のことですが) 状態内の状態変数 (「先頭のマイナス」、「位置 3 の小数点」) を暗黙的に記憶し、そのような値の代入、保存、およびその後のフェッチ/テストを回避することです。 。入力文字のみに対する単純な古い "if" ステートメントを使用してステート マシンを実装します (そのため、ツリーはネストされた if のセットになります)。バッファ文字へのインライン アクセス。関数呼び出しをしたくない getchar あなたを遅くするために。
先頭のゼロは単純に抑制できます。途方もなく長い先行ゼロシーケンスを処理するには、ここでループが必要になる場合があります。最初の非ゼロ桁は、アキュムレータをゼロにしたり、10 を乗算したりせずに収集できます。最初の 4 ～ 9 の非ゼロ桁 (16 ビットまたは 32 ビット整数の場合) は、整数と定数値 10 の積を使用して収集できます (ほとんどのコンパイラでは、数回のシフトと加算に変換されます)。[オーバートップ:ゼロ以外の数字が見つかるまで、ゼロの数字については何の作業も必要ありません。その後、N 個の連続するゼロに対する 10^N の乗算が必要になります。これらすべてをステート マシンに接続できます]。最初の 4 ～ 9 に続く数字は、マシンのワード サイズに応じて 32 ビットまたは 64 ビットの乗算を使用して収集される場合があります。精度は気にしないので、32 ビットまたは 64 ビット相当を収集した後は単純に数字を無視できます。アプリケーションがこれらの数値を実際にどのように処理するかに基づいて、ゼロ以外の桁の固定数が確保できたら実際に停止できると思います。数字列内に小数点が見つかった場合、ステート マシン ツリーに分岐が生じるだけです。このブランチはポイントの暗黙的な位置を知っているため、後で 10 のべき乗で適切にスケールする方法を知っています。このコードのサイズが気に入らない場合は、努力すれば、いくつかのステート マシンのサブツリーを組み合わせることができるかもしれません。
[オーバートップ:整数部分と小数部分を別々の (小さな) 整数として保持します。これには、整数部分と小数部分を結合するために最後に追加の浮動小数点演算が必要になりますが、おそらくその価値はありません。
[オーバートップ:数字ペアの 2 文字を 16 ビット値に収集し、16 ビット値を検索します。これにより、メモリ アクセスと引き換えにレジスタの乗算が回避されますが、おそらく最新のマシンでは成功しません。
「E」に遭遇すると、上記のように指数を整数として収集します。事前計算された乗数 (指数に「-」記号が存在する場合は逆数) のテーブルで、事前に計算されスケーリングされた 10 のべき乗を正確に検索し、収集された仮数を乗算します。(浮動小数点除算は決して行わないでください)。各指数収集ルーチンはツリーの異なるブランチ (リーフ) にあるため、10 の累乗インデックスをオフセットすることによって、小数点の見かけ上の位置または実際の位置を調整する必要があります。
[オーバートップ:～のコストを避けることができます ptr++ 数値の文字がバッファ内に直線的に格納され、バッファの境界を越えないことがわかっている場合。木の枝に沿った k 番目の状態では、次のように k 番目の文字にアクセスできます。 *(start+k). 。優れたコンパイラは通常、アドレッシング モードでインデックス付きオフセット内の "...+k" を隠すことができます。]
適切に実行すると、このスキームは、ゼロ以外の桁ごとにおよそ 1 回の安価な積和演算、1 回の仮数の浮動小数点へのキャスト、および 1 回の浮動小数点乗算を実行して、指数と小数点の位置によって結果をスケールします。
上記は実装していません。ループを使用したバージョンを実装しましたが、それらはかなり高速です。
いくつかのデータ交換ファイルの解析中に Winforms アプリケーションのパフォーマンスが非常に遅くなったのを覚えています。私たちは皆、データベース サーバーのスラッシングが原因だと考えていました。しかし、賢明な上司は実際に、ボトルネックが解析された文字列を次の形式に変換する呼び出しにあることを発見しました。小数！
最も単純な方法は、文字列内の各桁 (文字) をループし、現在までの合計を保持し、その合計に 10 を掛けて、次の桁の値を加算することです。文字列の終わりに到達するか、ドットに遭遇するまでこれを続けます。ドットに遭遇した場合は、整数部分を小数部分から分離し、各桁を 10 で割る乗数を用意します。続けて追加してください。
例：123.456
合計= 0を実行し、1（現在は1）の合計= 1 * 10 = 10を実行します。分数パーツ乗数= 0.1、乗算4、0.4を取得し、合計を実行し、123.4の乗数= 0.1 / 10 = 0.01になり、5を乗算し、0.05を取得し、合計を実行します。 6を掛け、0.006を取得し、合計を実行し、123.456を作成します
もちろん、数値の正しさだけでなく負の数値もテストすると、テストはさらに複雑になります。しかし、入力が正しいと「想定」できれば、コードをより単純かつ高速に作成できます。
GPU にこの作業を実行させることを検討したことがありますか?文字列を GPU メモリにロードしてすべて処理させることができれば、プロセッサよりも大幅に高速に実行できる優れたアルゴリズムが見つかる可能性があります。
あるいは、FPGA で実行します。任意のコプロセッサを作成するために使用できる FPGA PCI-E ボードがあります。DMA を使用して、変換する文字列の配列を含むメモリ部分を FPGA に指定し、変換された値を残して処理を実行させます。
クアッドコアプロセッサを見たことがありますか?これらのケースのほとんどにおける本当のボトルネックは、とにかくメモリ アクセスです。
-アダム
