c - C は sin() やその他の数学関数をどのように計算しますか?

Question

私は .NET 逆アセンブルと GCC ソースコードを詳しく調べましたが、実際の実装やその他の数学関数がどこにも見つからないようですsin()...それらは常に何か他のものを参照しているようです。

誰かがそれらを見つけるのを手伝ってくれますか? C が実行されるすべてのハードウェアがハードウェアで三角関数をサポートする可能性は低いので、ソフトウェアアルゴリズムがどこかにあるはずですよね?

関数を計算するいくつかの方法を知っており、楽しみのためにテイラー級数を使用して関数を計算する独自のルーチンを作成しました。私のアルゴリズムはかなり賢いと思いますが (明らかにそうではありません)、私の実装は常に数桁遅いため、実際の製品言語がどのようにそれを行うのか興味があります。

score 248 · Accepted Answer

GNU libm では、の実装sinはシステムに依存します。したがって、各プラットフォームの実装は、 sysdepsの適切なサブディレクトリのどこかにあります。

1 つのディレクトリーには、IBM 提供の C での実装が含まれています。sin()2011 年 10 月以降、これは一般的な x86-64 Linux システムで呼び出したときに実際に実行されるコードです。明らかにfsinアセンブリ命令よりも高速です。ソースコード: sysdeps/ieee754/dbl-64/s_sin.cを探します__sin (double x)。

このコードは非常に複雑です。1 つのソフトウェアアルゴリズムで、 x値の全範囲で可能な限り高速で正確なものはありません。そのため、ライブラリはいくつかの異なるアルゴリズムを実装し、最初の仕事はxを見て、どのアルゴリズムを使用するかを決定することです。

xが0 に非常に近い場合、sin(x) == xが正解です。
少し離れて、sin(x)おなじみのテイラーシリーズを使用しています。ただ、これは0付近でしか正確ではないので…
角度が約 7° を超える場合、別のアルゴリズムが使用され、sin(x) と cos(x) の両方のテイラー級数近似が計算され、事前に計算されたテーブルの値を使用して近似が調整されます。
いつ | × | sin> 2 の場合、上記のアルゴリズムはどれも機能しないため、コードはor のcos代わりに供給できる 0 に近い値を計算することから始めます。
xが NaN または無限大であることを処理する別のブランチがあります。

このコードは、これまでに見たことのないいくつかの数値ハックを使用していますが、浮動小数点の専門家の間ではよく知られているかもしれません。場合によっては、数行のコードを説明するのに数段落かかることがあります。たとえば、次の 2 行

double t = (x * hpinv + toint);
double xn = t - toint;

π/2 の倍数、具体的には× π/2だけxと異なる 0 に近い値にxを減らすために (時々) 使用されます。分割や分岐を行わずにこれを行う方法は、かなり巧妙です。しかし、コメントはありません！xn

GCC/glibc の古い 32 ビットバージョンではこのfsin命令が使用されていましたが、これは一部の入力に対して驚くほど不正確です。わずか 2 行のコードでこれを説明する魅力的なブログ投稿があります。

純粋な Cでの fdlibm の実装はsin、glibc の実装よりもはるかに単純で、適切にコメントされています。ソースコード: fdlibm/s_sin.cおよびfdlibm/k_sin.c

score 72 · Accepted Answer

サインやコサインなどの関数は、マイクロプロセッサ内のマイクロコードに実装されています。たとえば、インテルのチップには、これらの組み立て手順があります。AC コンパイラは、これらのアセンブリ命令を呼び出すコードを生成します。(対照的に、Java コンパイラはそうしません。Java は、ハードウェアではなくソフトウェアで三角関数を評価するため、実行速度が大幅に低下します。)

チップは、テイラー級数を使用して三角関数を計算することはありません(少なくとも完全には使用しません)。まずCORDICを使用しますが、短いテイラー級数を使用して CORDIC の結果を洗練したり、非常に小さな角度に対して高い相対精度で正弦を計算するなどの特別な場合に使用したりすることもあります。詳細については、このStackOverflow answerを参照してください。

score 66 · Accepted Answer

これは、経験の浅いソフトウェアエンジニアに対する私の最大の不満の 1 つです。彼らは、人生で誰もこれらの計算をしたことがないかのように、(テイラーの級数を使用して) 超越関数をゼロから計算します。違います。これは明確に定義された問題であり、非常に賢いソフトウェアおよびハードウェアエンジニアによって何千回もアプローチされており、明確に定義された解決策があります。基本的に、ほとんどの超越関数はチェビシェフ多項式を使用して計算します。どの多項式が使用されるかは、状況によって異なります。まず、この問題に関するバイブルは、Hart と Cheney による「Computerapproximations」という本です。その本では、ハードウェアの加算器、乗算器、除算器などがあるかどうかを判断し、どの操作が最も高速かを判断できます。たとえば、非常に高速な除算器がある場合、正弦を計算する最速の方法は、P1(x)/P2(x) である可能性があります。ここで、P1、P2 はチェビシェフ多項式です。高速除算器がないと、P(x) だけになる可能性があります。ここで、P は P1 または P2 よりもはるかに多くの項を持っているため、遅くなります。したがって、最初のステップは、ハードウェアとその機能を決定することです。次に、チェビシェフ多項式の適切な組み合わせを選択します (通常は、コサインに対して cos(ax) = aP(x) の形式になります。P はチェビシェフ多項式です)。次に、必要な小数精度を決定します。たとえば、7 桁の精度が必要な場合は、私が言及した本の適切な表でそれを調べると、(精度 = 7.33 の場合) 数 N = 4 と多項式数 3502 が得られます。N は次数です。 N=4 であるため、多項式 (p4.x^4 + p3.x^3 + p2.x^2 + p1.x + p0 です)。次に、p4、p3、p2、p1 の実際の値を調べます。本の後ろにある p0 値は 3502 未満です (浮動小数点になります)。次に、次の形式でアルゴリズムをソフトウェアに実装します: (((p4.x + p3).x + p2).x + p1).x + p0 ....これは、10 進数 7 のコサインを計算する方法です。そのハードウェアに配置します。

FPU での超越演算のほとんどのハードウェア実装には、通常、このようなマイクロコードと演算が含まれることに注意してください (ハードウェアによって異なります)。チェビシェフ多項式はほとんどの超越に使用されますが、すべてではありません。たとえば、最初にルックアップテーブルを使用してニュートンラフソン法を 2 回繰り返した方が平方根の方が高速です。繰り返しになりますが、その本「Computerapproximations」はそれを教えてくれます。

これらの機能を実装する予定がある場合は、その本を入手することをお勧めします。この種のアルゴリズムのバイブルです。コーディックなど、これらの値を計算するための代替手段がたくさんあることに注意してください。ただし、これらは、低い精度のみが必要な特定のアルゴリズムに最適な傾向があります。毎回精度を保証するには、チェビシェフ多項式が最適です。私が言ったように、明確に定義された問題。50年間解決されてきました...そして、それがどのように行われたかです。

そうは言っても、チェビシェフ多項式を使用して、低次多項式で単精度の結果を取得できる手法があります (上記のコサインの例のように)。次に、「Gal's Accurate Tables Method」など、はるかに大きな多項式を使用せずに精度を上げるために値間を補間する他の手法があります。この後者の手法は、ACM 文献を参照している投稿が参照しているものです。しかし、最終的には、チェビシェフ多項式を使用して 90% を達成します。

楽しみ。

score 16 · Accepted Answer

具体的にはsin、テイラー展開を使用すると、次のようになります。

sin(x) := x - x^3/3! + x^5/5! - x^7/7! + ... (1)

それらの間の差が許容レベルよりも低くなるか、有限の量のステップだけになるまで（より高速ですが、精度が低くなります）、項を追加し続けます。例は次のようになります。

float sin(float x)
{
  float res=0, pow=x, fact=1;
  for(int i=0; i<5; ++i)
  {
    res+=pow/fact;
    pow*=-1*x*x;
    fact*=(2*(i+1))*(2*(i+1)+1);
  }

  return res;
}

注: (1) は、小さな角度の近似 sin(x)=x のために機能します。より大きな角度では、許容できる結果を得るために、より多くの項を計算する必要があります。while 引数を使用して、一定の精度で続行できます。

double sin (double x){
    int i = 1;
    double cur = x;
    double acc = 1;
    double fact= 1;
    double pow = x;
    while (fabs(acc) > .00000001 &&   i < 100){
        fact *= ((2*i)*(2*i+1));
        pow *= -1 * x*x; 
        acc =  pow / fact;
        cur += acc;
        i++;
    }
    return cur;

}

score 13 · Accepted Answer

はい、計算用のソフトウェアアルゴリズムもありsinます。基本的に、デジタルコンピューターでこれらの種類のものを計算することは、通常、関数を表すテイラー級数を近似するなどの数値的方法を使用して行われます。

数値法は関数を任意の精度で近似できます。浮動小数点数の精度は有限であるため、これらのタスクに非常に適しています。

score 12 · Accepted Answer

のような三角関数に関してはsin()、5 年後、高品質の三角関数の重要な側面である範囲削減については言及されていcos()ません。tan()

これらの関数の初期段階では、角度をラジアン単位で 2*π 間隔の範囲に減らします。x = remainder(x, 2*M_PI)しかし、π は不合理なので、誤差を導入するなどの単純な削減M_PI、またはマシン pi は π の近似値です。それで、どうすればいいx = remainder(x, 2*π)ですか？

初期のライブラリは、拡張された精度または細工されたプログラミングを使用して高品質の結果を提供しましたが、それでも範囲は限られていましたdouble。のように大きな値が要求された場合sin(pow(2,30))、結果は無意味であるか0.0、または精度の完全な損失または精度の部分的な損失のようなエラーフラグが設定されていた可能性があります。TLOSSPLOSS

-π から π のような間隔への大きな値の適切な範囲縮小は、のような基本的な三角関数sin()自体の課題に匹敵する挑戦的な問題です。

優れたレポートは、巨大な引数の引数削減: 最後のビットまで有効(1992) です。問題を十分にカバーしています。必要性と、さまざまなプラットフォーム (SPARC、PC、HP、その他 30 以上) での状況について説明し、からまでのすべて に対して質の高い結果を提供するソリューションアルゴリズムを提供します。double-DBL_MAXDBL_MAX

元の引数が度単位であるが、大きな値である可能性がある場合は、fmod()精度を向上させるために first を使用します。良いものはエラーfmod()を導入しないため、優れた範囲削減を提供します。

// sin(degrees2radians(x))
sin(degrees2radians(fmod(x, 360.0))); // -360.0 < fmod(x,360) < +360.0

さまざまなトリガー ID があり、remquo()さらに改善されます。サンプル: sind()

score 12 · Accepted Answer

テイラー級数を使用して、級数の項間の関係を見つけようとするので、何度も何度も計算する必要はありません。

コサインの例を次に示します。

double cosinus(double x, double prec)
{
    double t, s ;
    int p;
    p = 0;
    s = 1.0;
    t = 1.0;
    while(fabs(t/s) > prec)
    {
        p++;
        t = (-t * x * x) / ((2 * p - 1) * (2 * p));
        s += t;
    }
    return s;
}

これを使用して、既に使用されているものを使用して合計の新しい項を取得できます (階乗と x ^2pは避けます)

score 11 · Accepted Answer

複雑な質問です。x86 ファミリの Intel のような CPU には、このsin()機能のハードウェア実装がありますが、これは x87 FPU の一部であり、64 ビットモードでは使用されなくなりました (代わりに SSE2 レジスタが使用されます)。そのモードでは、ソフトウェア実装が使用されます。

そのような実装がいくつかあります。1 つはfdlibmにあり、Java で使用されます。私の知る限り、glibc の実装には fdlibm の一部と、IBM によって提供されたその他の部分が含まれています。

などの超越関数のソフトウェア実装では、sin()通常、テイラー級数から得られることが多い多項式による近似を使用します。

score 6 · Accepted Answer

ライブラリ関数の実際の実装は、特定のコンパイラやライブラリプロバイダ次第です。それがハードウェアで行われるかソフトウェアで行われるか、テイラー展開であるかどうかなど、さまざまです。

それは絶対に役に立たないことを理解しています。

score 5 · Accepted Answer

ハードウェアではなくソフトウェアでの実装が必要な場合、この質問に対する決定的な答えを探す場所はNumerical Recipesの第 5 章です。私のコピーは箱に入っているので、詳細を説明することはできませんが、短いバージョン (これが正しいことを覚えていれば) はtan(theta/2)、原始的な操作として取り、そこから他の操作を計算することです。計算は級数近似で行われますが、テイラー級数よりもはるかに早く収束するものです。

申し訳ありませんが、本を手に取らなければ、これ以上思い出すことはできません。

score 5 · Accepted Answer

これらは通常、ソフトウェアで実装され、ほとんどの場合、対応するハードウェア (つまり、aseembly) 呼び出しを使用しません。ただし、ジェイソンが指摘したように、これらは実装固有です。

これらのソフトウェアルーチンはコンパイラソースの一部ではなく、clib や GNU コンパイラの glibc などの対応するライブラリに含まれていることに注意してください。http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functionsを参照してください。

より細かく制御したい場合は、正確に何が必要かを慎重に評価する必要があります。典型的な方法のいくつかは、ルックアップテーブルの補間、アセンブリ呼び出し (多くの場合低速)、または平方根のニュートンラフソンなどの他の近似スキームです。

score 5 · Accepted Answer

ソースにアクセスして、一般的に使用されているライブラリで誰かが実際にどのようにそれを行っているかを見るのに勝るものはありません。特に 1 つの C ライブラリの実装を見てみましょう。私はuLibCを選びました。

sin 関数は次のとおりです。

http://git.uClibc.org/uClibc/tree/libm/s_sin.c

これは、いくつかの特殊なケースを処理しているように見えます。次に、入力を [-pi/4,pi/4] の範囲にマップするために引数の削減を実行します (引数を大きな部分と末尾の 2 つの部分に分割します)。電話する前に

http://git.uClibc.org/uClibc/tree/libm/k_sin.c

次に、これら2つの部分で動作します。裾がない場合は、次数 13 の多項式を使用しておおよその答えが生成されます。sin(x+y) = sin(x) + sin'(x')y

score 4 · Accepted Answer

sin()現在のx86プロセッサ（Intel Core 2 Duoとしましょう）でGCCのCコンパイラでコンパイルされたCプログラムの場合について答えようとします。

C 言語では、標準 C ライブラリに一般的な数学関数が含まれていますが、言語自体には含まれていません (たとえばpow、べき乗sin、cos正弦、余弦のとなど)。そのヘッダーはmath.hに含まれています。

現在、GNU/Linux システムでは、これらのライブラリ関数は glibc (GNU libc または GNU C ライブラリ) によって提供されています。しかし、GCC コンパイラは、コンパイラフラグを使用して数学ライブラリ( libm.so)にリンクし、これらの数学関数を使用できるようにする必要があります。~~標準 C ライブラリに含まれていない理由がわかりません。~~これらは、浮動小数点関数のソフトウェアバージョン、つまり「ソフトフロート」になります。-lm

余談ですが、数学関数を分離する理由は歴史的なものであり、私の知る限り、おそらく共有ライブラリが利用可能になる前に、非常に古い Unix システムで実行可能なプログラムのサイズを縮小することを目的としていたにすぎません。

これで、コンパイラは標準 C ライブラリ関数sin()(によって提供されるlibm.so) を最適化して、CPU/FPU の組み込み sin() 関数へのネイティブ命令の呼び出しに置き換えることができます。これは、FPU 命令 ( FSINx86/x87 用)として存在します。 Core 2 シリーズのような新しいプロセッサ (これは、i486DX までさかのぼってほとんど正しいです)。これは、gcc コンパイラに渡される最適化フラグに依存します。コンパイラが、i386 以降のプロセッサで実行されるコードを書くように指示された場合、そのような最適化は行われません。この-mcpu=486フラグは、そのような最適化を行っても安全であることをコンパイラに通知します。

プログラムが sin() 関数のソフトウェアバージョンを実行する場合、CORDIC (座標回転デジタルコンピューター) またはBKMアルゴリズム、または計算に現在一般的に使用されている表計算またはベキ級数計算に基づいて実行します。そんな超越機能。[ソース: http://en.wikipedia.org/wiki/Cordic#Application]

最近の (およそ 2.9x 以降の) バージョンの gcc には、組み込みバージョンの sin も用意されています__builtin_sin()。これは、最適化として、C ライブラリバージョンへの標準呼び出しを置き換えるために使用されます。

これは泥のようにはっきりしていると思いますが、期待していたよりも多くの情報が得られることを願っています。

score 4 · Accepted Answer

これらの関数の実際の GNU 実装を C で見たい場合は、glibc の最新のトランクを調べてください。GNU C ライブラリを参照してください。

score 4 · Accepted Answer

多くの人が指摘したように、それは実装依存です。しかし、私があなたの質問を理解している限り、あなたは数学関数の実際のソフトウェア実装に興味を持っていましたが、それを見つけることができませんでした. その場合は、次のようになります。

http://ftp.gnu.org/gnu/glibc/から glibc ソースコードをダウンロードします。
展開されたglibc ルート\sysdeps\ieee754\dbl-64 フォルダーdosincos.cにあるファイルを確認します
同様に、数学ライブラリの残りの実装を見つけることができます。適切な名前のファイルを探すだけです

拡張子の付いたファイルを見ることもできます。その内容は、バイナリ形式でさまざまな関数の事前計算された値の.tbl巨大なテーブルにすぎません。それが、実装が非常に高速な理由です。使用する系列のすべての係数を計算する代わりに、クイックルックアップを行うだけで、はるかに高速です。ところで、彼らはテーラー級数を使用してサインとコサインを計算します。

これが役立つことを願っています。

score 4 · Accepted Answer

そのような関数が評価されるときは常に、あるレベルで次のいずれかが存在する可能性が最も高くなります。

補間された値のテーブル (高速で不正確なアプリケーション用 - コンピュータグラフィックスなど)
目的の値に収束する級数の評価 --- おそらくテイラー級数ではなく、Clenshaw-Curtis のような派手な求積法に基づくものである可能性が高いです。

ハードウェアサポートがない場合、コンパイラはおそらく後者の方法を使用し、ac ライブラリを使用するのではなく、アセンブラコードのみ (デバッグシンボルなし) を発行します。これにより、デバッガで実際のコードを追跡するのが難しくなります。

score 2 · Accepted Answer

サイン/コサイン/タンジェントの計算は、実際にはテイラー級数を使用したコードで非常に簡単に実行できます。自分で書くのに5秒ほどかかります。

プロセス全体は、次の式で要約できます。

私が C 用に書いたいくつかのルーチンを次に示します。

double _pow(double a, double b) {
    double c = 1;
    for (int i=0; i<b; i++)
        c *= a;
    return c;
}

double _fact(double x) {
    double ret = 1;
    for (int i=1; i<=x; i++) 
        ret *= i;
    return ret;
}

double _sin(double x) {
    double y = x;
    double s = -1;
    for (int i=3; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _cos(double x) {
    double y = 1;
    double s = -1;
    for (int i=2; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _tan(double x) {
     return (_sin(x)/_cos(x));  
}

c - C は sin() やその他の数学関数をどのように計算しますか?

22 に答える 22

Related

Reference