私はこれが繰り返しの質問であることを知っていますが、私はまだ実際に有用な答えを見つけていません。私は基本的acosにC++の関数の高速近似を探していますが、標準の関数を大幅に上回ることができるかどうかを知りたいです。
しかし、あなた方の何人かは私の特定の問題について洞察を持っているかもしれません：私は非常に速くなければならない科学的なプログラムを書いています。メインアルゴリズムの複雑さは、次の式を計算することに要約されます（多くの場合、異なるパラメーターを使用）。
sin( acos(t_1) + acos(t_2) + ... + acos(t_n) )
ここで、t_iは既知の実数（double）であり、n非常に小さい（6より小さいなど）。少なくとも1e-10の精度が必要です。現在、標準関数sinとacosC++関数を使用しています。
どういうわけか私はかなりスピードを上げることができると思いますか？t_i数学を知っている人にとっては、 （平方根のみを含む）代数式を取得するために、その正弦を拡張するのが賢明だと思いますか？
回答ありがとうございます。
以下のコードは、精度要件を満たし、試してみたいと思われる、sin()の簡単な実装を提供します。acos()プラットフォームでの数学ライブラリの実装は、そのプラットフォームの特定のハードウェア機能に合わせて高度に調整されている可能性が高く、効率を最大化するためにアセンブリでコード化されている可能性が高いため、ハードウェアの詳細に対応していない単純なコンパイル済みCコードでは提供されない可能性があります。精度要件が完全な倍精度からいくらか緩和されている場合でも、より高いパフォーマンス。Viktor Latypovが指摘しているように、超越数学関数への高価な呼び出しを必要としないアルゴリズムの代替案を探すことも価値があるかもしれません。
以下のコードでは、シンプルでポータブルな構造に固執しようとしました。コンパイラがrint()[C99およびC++11で指定された]関数をサポートしている場合は、の代わりにそれを使用することをお勧めしますmy_rint()。一部のプラットフォームではfloor()、マシンの状態を動的に変更する必要があるため、への呼び出しにコストがかかる場合があります。関数my_rint()、、、およびはsin_core()、最高のパフォーマンスを得るためにインライン化する必要があります。コンパイラは、これを高い最適化レベルで自動的に実行する場合があります（たとえば、でコンパイルする場合）。または、ツールチェーンに応じて、これらの関数に適切なインライン属性を追加できます（たとえば、inlineまたは__inline）。cos_core()asin_core()-O3
プラットフォームについて何も知らないので、単純な多項式近似を選択しました。これは、エストリンのスキームとホーナーのスキームを使用して評価されます。これらの評価スキームの説明については、ウィキペディアを参照してください。
http://en.wikipedia.org/wiki/Estrin%27s_scheme、http://en.wikipedia.org/wiki/Horner_scheme _ _
近似自体はミニマックスタイプであり、Remezアルゴリズムを使用してこの回答用にカスタム生成されました。
http://en.wikipedia.org/wiki/Minimax_approximation_algorithm、http://en.wikipedia.org/wiki/Remez_algorithm _ _
次の本で説明されているように、私はCody / Waiteスタイルの引数削減を使用したため、の引数削減で使用されたIDacos()はコメントに記載されています。sin()
WJ Cody、W。Waite、初等機能のソフトウェアマニュアル。プレンティスホール、1980年
コメントに記載されている誤差範囲は概算であり、厳密にテストまたは証明されていません。
/* not quite rint(), i.e. results not properly rounded to nearest-or-even */
double my_rint (double x)
{
double t = floor (fabs(x) + 0.5);
return (x < 0.0) ? -t : t;
}
/* minimax approximation to cos on [-pi/4, pi/4] with rel. err. ~= 7.5e-13 */
double cos_core (double x)
{
double x8, x4, x2;
x2 = x * x;
x4 = x2 * x2;
x8 = x4 * x4;
/* evaluate polynomial using Estrin's scheme */
return (-2.7236370439787708e-7 * x2 + 2.4799852696610628e-5) * x8 +
(-1.3888885054799695e-3 * x2 + 4.1666666636943683e-2) * x4 +
(-4.9999999999963024e-1 * x2 + 1.0000000000000000e+0);
}
/* minimax approximation to sin on [-pi/4, pi/4] with rel. err. ~= 5.5e-12 */
double sin_core (double x)
{
double x4, x2, t;
x2 = x * x;
x4 = x2 * x2;
/* evaluate polynomial using a mix of Estrin's and Horner's scheme */
return ((2.7181216275479732e-6 * x2 - 1.9839312269456257e-4) * x4 +
(8.3333293048425631e-3 * x2 - 1.6666666640797048e-1)) * x2 * x + x;
}
/* minimax approximation to arcsin on [0, 0.5625] with rel. err. ~= 1.5e-11 */
double asin_core (double x)
{
double x8, x4, x2;
x2 = x * x;
x4 = x2 * x2;
x8 = x4 * x4;
/* evaluate polynomial using a mix of Estrin's and Horner's scheme */
return (((4.5334220547132049e-2 * x2 - 1.1226216762576600e-2) * x4 +
(2.6334281471361822e-2 * x2 + 2.0596336163223834e-2)) * x8 +
(3.0582043602875735e-2 * x2 + 4.4630538556294605e-2) * x4 +
(7.5000364034134126e-2 * x2 + 1.6666666300567365e-1)) * x2 * x + x;
}
/* relative error < 7e-12 on [-50000, 50000] */
double my_sin (double x)
{
double q, t;
int quadrant;
/* Cody-Waite style argument reduction */
q = my_rint (x * 6.3661977236758138e-1);
quadrant = (int)q;
t = x - q * 1.5707963267923333e+00;
t = t - q * 2.5633441515945189e-12;
if (quadrant & 1) {
t = cos_core(t);
} else {
t = sin_core(t);
}
return (quadrant & 2) ? -t : t;
}
/* relative error < 2e-11 on [-1, 1] */
double my_acos (double x)
{
double xa, t;
xa = fabs (x);
/* arcsin(x) = pi/2 - 2 * arcsin (sqrt ((1-x) / 2))
* arccos(x) = pi/2 - arcsin(x)
* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2))
*/
if (xa > 0.5625) {
t = 2.0 * asin_core (sqrt (0.5 * (1.0 - xa)));
} else {
t = 1.5707963267948966 - asin_core (xa);
}
/* arccos (-x) = pi - arccos(x) */
return (x < 0.0) ? (3.1415926535897932 - t) : t;
}
sin( acos(t1) + acos(t2) + ... + acos(tn) )
要約すると、
sin( acos(x) ) and cos(acos(x))=x
なぜなら
sin(a+b) = cos(a)sin(b)+sin(a)cos(b).
まずは
sin( acos(x) ) = sqrt(1-x*x)
平方根のテイラー級数展開は、問題を多項式計算に減らします。
明確にするために、これがn = 2、n=3への拡張です。
sin( acos(t1) + acos(t2) ) = sin(acos(t1))cos(acos(t2)) + sin(acos(t2))cos(acos(t1) = sqrt(1-t1*t1) * t2 + sqrt(1-t2*t2) * t1
cos( acos(t2) + acos(t3) ) = cos(acos(t2)) cos(acos(t3)) - sin(acos(t2))sin(acos(t3)) = t2*t3 - sqrt(1-t2*t2)*sqrt(1-t3*t3)
sin( acos(t1) + acos(t2) + acos(t3)) =
sin(acos(t1))cos(acos(t2) + acos(t3)) + sin(acos(t2)+acos(t3) )cos(acos(t1)=
sqrt(1-t1*t1) * (t2*t3 - sqrt(1-t2*t2)*sqrt(1-t3*t3)) + (sqrt(1-t2*t2) * t3 + sqrt(1-t3*t3) * t2 ) * t1
等々。
（-1,1）のxのsqrt（）は、次を使用して計算できます。
x_0 is some approximation, say, zero
x_(n+1) = 0.5 * (x_n + S/x_n) where S is the argument.
編集：私は「バビロニア法」を意味します。詳細については、ウィキペディアの記事を参照してください。（0,1）のxで1e-10を達成するには、5〜6回の反復が必要です。
Jonas Wielickiがコメントで述べているように、正確なトレードオフはそれほど多くありません。
最善の策は、関数のプロセッサ組み込み関数を試して使用し（コンパイラがまだこれを行っていない場合）、必要な計算の量を減らすためにいくつかの数学を使用することです。
また、すべてをCPUに適した形式に保つこと、キャッシュミスが少ないことなどを確認することも非常に重要です。
おそらくGPUに移行するなど、大量の関数を計算している場合acosは、オプションがありますか？
ルックアップテーブルを作成して、標準のc ++関数の代わりに使用して、パフォーマンスが向上するかどうかを確認できます。
メモリを調整し、データをカーネルにストリーミングすることで、大幅な向上を実現できます。ほとんどの場合、これは数学関数を再作成することによって得られる利益を小さくします。カーネルオペレーターとの間のメモリアクセスを改善する方法を考えてください。
バッファリング技術を使用することにより、メモリアクセスを改善できます。これは、ハードウェアプラットフォームによって異なります。これをDSPで実行している場合は、データをL2キャッシュにDMAし、乗算器ユニットが完全に占有されるように命令をスケジュールできます。
汎用CPUを使用している場合、できることのほとんどは、整列されたデータを使用し、プリフェッチによってキャッシュラインにフィードすることです。ネストされたループがある場合は、キャッシュラインが使用されるように、最も内側のループを前後に移動する必要があります（つまり、前方に反復してから後方に反復します）。
複数のコアを使用して計算を並列化する方法も考えられます。GPUを使用できる場合、これによりパフォーマンスが大幅に向上する可能性があります（精度は低下しますが）。
他の人が言ったことに加えて、速度最適化のいくつかのテクニックがあります：
コードのどこでほとんどの時間が費やされているかを調べます。モースのメリットを得るには、その領域のみを最適化してください。
プロセッサは、実行パスの分岐やジャンプ、変更を好みません。一般に、プロセッサは命令パイプラインをリロードする必要があり、計算に費やすことができる時間を使い果たします。これには関数呼び出しが含まれます。
この手法は、ループ内により多くの「セット」の操作を配置し、反復回数を減らすことです。
頻繁に使用される変数は、として宣言する必要がありますregister。SOの多くのメンバーは、コンパイラーがこの提案を無視すると述べていますが、私はそうではないことを知りました。最悪の場合、入力に時間を浪費しました。
多くのプロセッサには、小さなforループを保持するのに十分なスペースが命令パイプラインにあります。これにより、命令パイプラインのリロードにかかる時間が短縮されます。
大きな計算ループを多くの小さなループに分散させます。
多くのプロセッサにはデータキャッシュがあり、これはプロセッサに非常に近い超高速メモリです。プロセッサは、オフプロセッサメモリからデータキャッシュを1回ロードするのが好きです。負荷が増えると、計算に費やすことができる時間が必要になります。Webで「データ指向のデザインキャッシュ」を検索します。
計算の設計を変更して、複数のプロセッサでの使用に簡単に適応できるようにします。多くのCPUには、命令を並行して実行できる複数のコアがあります。一部のプロセッサには、命令を複数のコアに自動的に委任するのに十分なインテリジェンスがあります。
一部のコンパイラは、並列処理用にコードを最適化できます（コンパイラのコンパイラオプションを検索してください）。並列処理用のコードを設計すると、コンパイラーにとってこの最適化が容易になります。
関数のアセンブリ言語リストを印刷します。関数の設計をアセンブリ言語の設計と一致するように、またはコンパイラがより最適なアセンブリ言語を生成できるように変更します。
本当に効率を上げる必要がある場合は、アセンブリ言語を最適化し、インラインアセンブリコードまたは別のモジュールとして挿入します。私は一般的に後者を好みます。
あなたの状況では、テイラー展開の最初の10項を取り、それらを別々に計算して、個々の変数に配置します。
double term1, term2, term3, term4;
double n, n1, n2, n3, n4;
n = 1.0;
for (i = 0; i < 100; ++i)
{
n1 = n + 2;
n2 = n + 4;
n3 = n + 6;
n4 = n + 8;
term1 = 4.0/n;
term2 = 4.0/n1;
term3 = 4.0/n2;
term4 = 4.0/n3;
次に、すべての用語を要約します。
result = term1 - term2 + term3 - term4;
// Or try sorting by operation, if possible:
// result = term1 + term3;
// result -= term2 + term4;
n = n4 + 2;
}
最初に2つの用語を考えてみましょう。
cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)
またcos(a+b) = cos(a)*cos(b) - sqrt(1-cos(a)*cos(a))*sqrt(1-cos(b)*cos(b))
cosをRHSに連れて行く
a+b = acos( cos(a)*cos(b) - sqrt(1-cos(a)*cos(a))*sqrt(1-cos(b)*cos(b)) )... 1
ここで、cos（a）= t_1およびcos（b）= t_2 a = acos（t_1）およびb = acos（t_2）
式（1）に代入すると、次のようになります。
acos(t_1) + acos(t_2) = acos(t_1*t_2 - sqrt(1 - t_1*t_1) * sqrt(1 - t_2*t_2))
ここでは、2つのacosを1つに結合したことがわかります。したがって、すべてのacosを再帰的にペアにして、バイナリツリーを形成できます。最後に、にsin(acos(x))等しい形式の式が残りますsqrt(1 - x*x)。
これにより、時間の複雑さが改善されます。
ただし、sqrt（）の計算の複雑さについてはよくわかりません。