float sinx(float x)
{
static const float a[] = {-.1666666664,.0083333315,-.0001984090,.0000027526,-.0000000239};
float xsq = x*x;
float temp = x*(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq);
return temp;
}
これらの定数はどのように計算されますか? この方法を計算cosしtanて使用する方法は？これを拡張して精度を上げることはできますか? もっと定数を追加する必要があると思いますか？
同じ次数のテイラー多項式に対する上記の「高速」サインの誤差のプロット。
この記事の執筆時点でのほぼすべての回答は、関数 sin のテイラー展開に言及していますが、関数の作成者が真剣であれば、テイラー係数を使用しないでしょう。テイラー係数は、ゼロに近づくと必要以上に良くなり、ゼロから離れるほど悪くなる多項式近似を生成する傾向があります。通常、目標は、-π/2…π/2 などの範囲で一様に適切な近似を取得することです。多項式近似の場合、これは Remez アルゴリズムを適用することで取得できます。現実的な説明はこの投稿です。
両方の多項式が同じ関数を近似しようとしているため、その方法で得られた多項式係数はテイラー係数に近くなりますが、多項式は同じ数の操作に対してより正確であるか、同じ (均一な) 品質に対してより少ない操作を必要とする場合があります。近似の。
あなたの質問の係数が正確にテイラー係数であるか、Remez アルゴリズムによって得られたわずかに異なる係数であるかは、それらを見るだけではわかりませんが、そうでなくても、おそらく使用されるべきものでした。
最後に、書いた人は誰でも、(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq)より良い多項式評価スキームについて読む必要があります。たとえば、Horner の:
1 + xsq*(a[0]+ xsq*(a[1] + xsq*(a[2] + xsq*(a[3] + xsq*a[4]))))N 2 /2の代わりに N 乗算を使用します。
それらは-1/6、、、1/120など-1/5040です。
というか: -1/3!, 1/5!, -1/7!, 1/9!... など
ここで、sin x のテイラー級数を見てください。
そのすぐ下に cos x があります。
cos x の場合、上の図からわかるように、定数は-1/2!, 1/4!, -1/6!, 1/8!...
tan x はわずかに異なります。
これを cosx 用に調整するには、次のようにします。
float cosx(float x)
{
static const float a[] = {-.5, .0416666667,-.0013888889,.0000248016,-.0000002756};
float xsq = x*x;
float temp = (1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq);
return temp;
}
係数は、 Handbook of Mathematical Functions , edに記載されているものと同じです。Abramowitz と Stegan (1964)、76 ページ、および Carlson と Goldstein による関数の合理的近似、Los Alamos Scientific Laboratory (1955) によるものです。
最初のものはhttp://www.jonsson.eu/resources/hmf/pdfwrite_600dpi/hmf_600dpi_page_76.pdfにあります。
2 つ目はhttp://www.osti.gov/bridge/servlets/purl/4374577-0deJO9/4374577.pdfです。37ページには次のように記載されています。
3 番目の質問である「精度を上げるためにこれを拡張できますか?」については、 http: //lol.zoy.org/wiki/doc/maths/remez に Remez アルゴリズムのダウンロード可能な C++ 実装があります。の6次多項式の係数を提供します（私はチェックしていません）sin：
error: 3.9e-14
9.99999999999624e-1
-1.66666666660981e-1
8.33333330841468e-3
-1.98412650240363e-4
2.75568408741356e-6
-2.50266363478673e-8
1.53659375573646e-10
またはもちろん、改善を実現するには float から double に変更する必要があります。cosそして、これはとに関する 2 番目の質問にも答えるかもしれませんtan。
また、コメントには、最終的に定点の回答が必要であることがわかります。私は約 26 年前に 8031 アセンブラーで 32 ビット固定小数点バージョンを実装しました。掘り下げて、有用なものがないかどうかを確認します。
更新: 32 ビットの double に固執している場合、「1 桁または 2 桁」の精度を上げる唯一の方法は、浮動小数点を忘れて固定小数点を使用することです。驚いたことに、Google は何も表示しないようです。次のコードは、標準の Linux マシンで実行される概念実証を提供します。
#include <stdio.h>
#include <math.h>
#include <stdint.h>
// multiply two 32-bit fixed-point fractions (no rounding)
#define MUL32(a, b) ((uint64_t)(a) * (b) >> 32)
// sin32: Fixed-point sin calculation for first octant, coefficients from
// Handbook for Computing Elementary Functions, by Lyusternik et al, p. 89.
// input: 0 to 0xFFFFFFFF, giving fraction of octant 0 to PI/8, relative to 2**32
// output: 0 to 0.7071, relative to 2**32
static uint32_t sin32(uint32_t x) { // x in 1st octant, = radians/PI*8*2**32
uint32_t y, x2 = MUL32(x, x); // x2 = x * x
y = 0x000259EB; // a7 = 0.000 035 877 1
y = 0x00A32D1E - MUL32(x2, y); // a5 = 0.002 489 871 8
y = 0x14ABBA77 - MUL32(x2, y); // a3 = 0.080 745 367 2
y = 0xC90FDA73u - MUL32(x2, y); // a1 = 0.785 398 152 4
return MUL32(x, y);
}
int main(void) {
int i;
for (i = 0; i < 45; i += 2) { // 0 to 44 degrees
const double two32 = 1LL << 32;
const double radians = i * M_PI / 180;
const uint32_t octant = i / 45. * two32; // fraction of 1st octant
printf("%2d %+.10f %+.10f %+.10f %+.0f\n", i,
sin(radians) - sin32(octant) / two32,
sin(radians) - sinf(radians),
sin(radians) - (float)sin(radians),
sin(radians) * two32 - sin32(octant));
}
return 0;
}
係数は、Lyusternik et alによるComputing Elementary Functions のハンドブックからのものです。89、こちら：
この特定の関数を選択した唯一の理由は、元のシリーズより項が 1 つ少ないためです。
結果は次のとおりです。
0 +0.0000000000 +0.0000000000 +0.0000000000 +0
2 +0.0000000007 +0.0000000003 +0.0000000012 +3
4 +0.0000000010 +0.0000000005 +0.0000000031 +4
6 +0.0000000012 -0.0000000029 -0.0000000011 +5
8 +0.0000000014 +0.0000000011 -0.0000000044 +6
10 +0.0000000014 +0.0000000050 -0.0000000009 +6
12 +0.0000000011 -0.0000000057 +0.0000000057 +5
14 +0.0000000006 -0.0000000018 -0.0000000061 +3
16 -0.0000000000 +0.0000000021 -0.0000000026 -0
18 -0.0000000005 -0.0000000083 -0.0000000082 -2
20 -0.0000000009 +0.0000000095 -0.0000000107 -4
22 -0.0000000010 -0.0000000007 +0.0000000139 -4
24 -0.0000000009 -0.0000000106 +0.0000000010 -4
26 -0.0000000005 +0.0000000065 -0.0000000049 -2
28 -0.0000000001 -0.0000000032 -0.0000000110 -0
30 +0.0000000005 -0.0000000126 -0.0000000000 +2
32 +0.0000000010 +0.0000000037 -0.0000000025 +4
34 +0.0000000015 +0.0000000193 +0.0000000076 +7
36 +0.0000000013 -0.0000000141 +0.0000000083 +6
38 +0.0000000007 +0.0000000011 -0.0000000266 +3
40 -0.0000000005 +0.0000000156 -0.0000000256 -2
42 -0.0000000009 -0.0000000152 -0.0000000170 -4
44 -0.0000000005 -0.0000000011 -0.0000000282 -2
したがって、この固定小数点計算はまたはよりも約10 倍正確であり、29 ビットまで正確であることがわかります。切り捨てではなく丸めを使用しても、わずかな改善しかありませんでした。sinf()(float)sin()MUL32()
その関数はsin、テイラー展開を使用しての値を計算しています。
これらの定数は、さまざまな -1/3!、1/5! です。など（テイラー級数の他の関数については、ここを参照）。
さて、級数のすべての項を指定した場合、sin(x) のテイラー展開はすべての x に対して正確ですが、ソフトウェアで三角関数を決定するためのより高速で正確な方法があります。
また、多くのプロセッサは、プロセッサに直接実装されたそのような機能を提供します (たとえば、x86 には既製のオペコードがあります)。
cos(x) = sqrt(1 - sin^2(x))
tan(x) = sin(x)/cos(x)
Sin(x) = x -x^3/3! + x^5/5! + (-1)^k*x^(2k+1)/(2k+1)! , k = 1, 2, ...
これは無限関数です