アルゴリズム

arctanはどのように実装されていますか？



How Is Arctan Implemented



解決：

x86プロセッサでのFPATAN命令の実装は、通常、独自仕様です。 arctanまたは他の（逆）三角関数を計算するために、一般的なアルゴリズムは3つのステップのプロセスに従います。

  1. 完全な入力ドメインを狭い間隔にマッピングするための引数の削減
  2. 狭い区間（一次近似区間）でのコア近似の計算
  3. 最終結果を生成するための引数の削減に基づく中間結果の拡張

引数の削減は通常、MathWorld（http://mathworld.wolfram.com/InverseTangent.html）などのさまざまな標準参照で検索できるよく知られた三角関数公式に基づいています。 arctanの計算では、一般的に使用されるIDは次のとおりです。



  • arctan（-x）= -arctan（x）
  • arctan（1 / x）= 0.5 * pi-arctan（x）[x> 0]
  • arctan（x）= arctan（c）+ arctan（（x --c）/（1 + x * c））

最後のアイデンティティは、値のテーブルの構築に役立つことに注意してくださいarctan（i / 2NS）、i = 1 ... 2NS、これにより、追加のテーブルストレージを犠牲にして、任意に狭い一次近似間隔を使用できます。これは、空間と時間の間の古典的なプログラミングのトレードオフです。

コア間隔の近似は、通常、十分な次数のミニマックス多項式近似です。有理数近似は、浮動小数点除算のコストが高いため、通常、最新のハードウェアでは競合しません。また、2つの多項式の計算に加えて除算によって生じる誤差のために、追加の数値誤差が発生します。



ミニマックス多項式近似の係数は、通常、Remezアルゴリズム（http://en.wikipedia.org/wiki/Remez_algorithm）を使用して計算されます。 MapleやMathematicaのようなツールには、そのような近似を計算するための機能が組み込まれています。多項式近似の精度は、すべての係数が正確に表現可能なマシン番号であることを確認することで改善できます。私が知っている唯一のツールには、このための組み込み機能があります。Sollya（http://sollya.gforge.inria.fr/）は、fpminimax（）関数。

多項式の評価では、通常、効率的で正確なホーナー法（http://en.wikipedia.org/wiki/Horner%27s_method）、またはエストリン法（http://en.wikipedia.org/wiki/）を組み合わせて使用​​します。 Estrin％27s_scheme）およびHorner's。 Estrinのスキームでは、スーパースカラープロセッサが提供する命令レベルの並列性をうまく利用できますが、全体的な命令数にはわずかな影響があり、精度にはほとんど影響しません（常にではありません）。

FMA（fused-multiply add）を使用すると、丸めステップの数が減り、減算キャンセルに対する保護が提供されるため、いずれかの評価スキームの精度とパフォーマンスが向上します。 FMAは、GPUや最近のx86CPUを含む多くのプロセッサで使用されています。標準Cおよび標準C ++では、FMA操作は次のように公開されます。fma（）標準ライブラリ関数。ただし、ハードウェアサポートを提供しないプラットフォームでエミュレートする必要があるため、これらのプラットフォームでは速度が低下します。



プログラミングの観点から、近似と引数の削減に必要な浮動小数点定数をテキスト表現からマシン表現に変換するときに、変換エラーのリスクを回避したいと考えています。 ASCIIから浮動小数点への変換ルーチンは、トリッキーなバグが含まれていることで有名です（例：http：//www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/）。標準Cによって提供される1つのメカニズム（ いいえ 私が知っているC ++は、プロプライエタリ拡張としてのみ利用可能です）、浮動小数点定数を、基になるビットパターンを直接表現する16進リテラルとして指定し、複雑な変換を効果的に回避することです。

以下は、上記の設計原理と手法の多くを示す倍精度arctan（）を計算するためのCコードです。この迅速に構築されたコードは、他の回答で指摘されている実装の洗練度に欠けていますが、2 ulps未満のエラーで結果を提供するはずです。これは、さまざまなコンテキストで十分な場合があります。すべての中間ステップに1024ビット浮動小数点演算を使用するRemezアルゴリズムの簡単な実装を使用して、カスタムミニマックス近似を作成しました。 Sollyaまたは同様のツールを使用すると、数値的に優れた近似値が得られると思います。

double my_atan（double x）{double a、z、p、r、s、q、o; / *引数の削減：arctan（-x）= -arctan（x）; arctan（1 / x）= 1/2 * pi-arctan（x）、x> 0の場合* / z = fabs（x）; a =（z> 1.0）？ 1.0 / z：z; / *ミニマックス多項式近似を評価します* / s = a * a; // a ** 2 q = s * s; // a ** 4 o = q * q; // a ** 8 / *低次の項にEstrinのスキームを使用* / p = fma（fma（fma（-0x1.53e1d2a25ff34p-16、s、0x1.d3b63dbb65af4p-13）、q、fma（-0x1.312788dde0801p -10、s、0x1.f9690c82492dbp-9））、o、fma（fma（-0x1.2cf5aabc7cef3p-7、s、0x1.162b0b2a3bfcep-6）、q、fma（-0x1.a7256feb6fc5cp-6、s、0x1。 171560ce4a483p-5）））; / *高次項にホーナー法を使用* / p = fma（fma（fma（fma（fma（fma（fma（fma（fma（fma（fma（fma（p、s、-0x1.4f44d841450e1p-5）、 s、0x1.7ee3d3f36bb94p-5）、s、-0x1.ad32ae04a9fd1p-5）、s、0x1.e17813d66954fp-5）、s、-0x1.11089ca9a5bcdp-4）、s、0x1.3b12b2db51738p-4）、s、- 0x1.745d022f8dc5cp-4）、s、0x1.c71c709dfe927p-4）、s、-0x1.2492491fa1744p-3）、s、0x1.99999999840d2p-3）、s、-0x1.555555555544cp-2）* s、a、a ）; / *引数の削減に基づく逆置換* / r =（z> 1.0）？ （0x1.921fb54442d18p + 0-p）：p; copysign（r、x）;を返します。 } 
三角関数にはかなり醜い実装があり、ハッキーでビットをいじくり回します。実際に使われているアルゴリズムを説明できる人をここで見つけるのはかなり難しいと思います。
 
atan2の実装は次のとおりです：https：//sourceware.org/git/？p = glibc.git; a = blob; f = sysdeps / ieee754 / dbl-64 / e_atan2.c; h = a287ca6656b210c77367eec3c46d72f18476d61d; hb = HEAD
 
編集：実際に私はこれを見つけました：http：//www.netlib.org/fdlibm/e_atan2.cこれははるかに簡単ですが、おそらくそのために遅くなります（？）。
 
FPUは一部の回路でこれらすべてを実行するため、CPUがこのすべての作業を実行する必要はありません。
  
要約：それは難しい。また、時々SOをぶらぶらしているEricPostpischilとStephenCanonは、それが非常に得意です。
 
多くの特殊機能の通常のアプローチは次のとおりです。
 
     
  • NaN、無限大、および符号付きゼロを特殊なケースとして処理します。
    •  
  • 数が多すぎて結果が次のように丸められる場合M_PI、リターンM_PI。このしきい値を呼び出すNS。
    •  
  • 何らかの種類の引数削減IDがある場合は、それを使用して引数をより適切な範囲にします。 （（ これは注意が必要です ： にとって罪とcos、これはあなたが倍数を選ぶことを意味します ちょうど 正しい範囲に着陸するように2piの値。）
    •  
  • 別れる[0、M）を有限の数の間隔に。各区間でかなり高次のアークタンにチェビシェフ近似を使用します。 （これはオフラインで行われ、通常、これらの実装で見られるすべての魔法の数のソースです。また、Remezの交換アルゴリズムを使用してチェビシェフ近似をわずかに厳密にすることができますが、これが大いに役立つケースはわかりません。 。）
    •  
  • 引数がどの間隔にあるかを把握します（を使用してifs and stuffまたはテーブルインデックスのトリック）、その間隔でChebyshevシリーズを評価します。
    •  
 
ここでは、いくつかのプロパティが特に望ましいです。
 
     
  • NSarctanの実装は単調である必要があります。つまり、NSarctan（x）<= arctan(y).
    •  
  • NSarctanの実装では、常に正しい答えから1ulp以内の答えを返す必要があります。これは相対誤差限界であることに注意してください。
    •  
 
これらの2つの特性が成り立つように、チェビシェフ級数を評価することは完全に簡単ではありません。 2つのトリックダブルスは、単一の値のさまざまな部分を表すために使用されます。ここでは一般的です。次に、実装が単調であることを示すいくつかのケースワークがおそらくあります。また、ゼロに近い、テイラー近似チェビシェフ近似の代わりにarctan ---相対誤差限界を求めており、ホーナーの法則を使用して級数を評価すると機能するはずです。
 
あなたが探しているならatanの実装を読むと、fdlibmは現在glibcにあるものよりも厄介ではないようです。引数の削減は、三角法のアイデンティティに基づいているようですtan（a + b）=（tan（a）+ tan（b））/（1-tan（a）tan（b））、使用0.5、1、または1.5の場合必要に応じてtan（a）。