Sine, Cosine 続

暇を利用して調査を継続。

角度が十分に小さい場合には引数のradian値とsinの結果がかなり近くなるというのは、プロットしてみたら原点付近は直線に近いので視覚的に納得出来た。degreeで10度未満ならdevmasterに投稿されていた近似方法より誤差が少なかった。
http://en.wikipedia.org/wiki/Small-angle_approximation

7世紀のインド人の数学者 Bhaskara I 西暦14世紀生まれのMadhava of Sangamagrama が見つけた近似式が単純で実装がしやすい。（power series expansion of sin x）
http://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximation_formula
http://en.wikipedia.org/wiki/Madhava_of_Sangamagrama#Trigonometry
http://en.wikipedia.org/wiki/Madhava%27s_sine_table
sin q = q – q3/3! + q5/5! – q7/7! + ...
cos q = 1 – q2/2! + q4/4! – q6/6! + ...

• • -

実際に値を出力してみると、角度が低いうちは誤差が少ないけれど角度が大きくなると誤差がどんどん増えてくる。対称性を利用して反転すればその問題は軽減出来るけれど分岐無しで演算するにはトリッキーなコードが必要かもしれない。

あと精度を上げるために次数を増やすと乗算回数が増える。5乗までしか使わない方法だと90度ぐらいになるとdevmasterの近似方法より精度が落ちる。

http://web.archive.org/web/20140528163441/http://lolengine.net/blog/2011/12/21/better-function-approximations

誤差の量が角度によって異なっているので、その分布を均等にする為に係数を調整した関数が提示されていた。最大誤差を最小化するMinimax Approximationの係数を使っていて、算出にはRemez exchange algorithmを使うとの事。

RemezアルゴリズムはFIRフィルタの設計でも使われていた。

ライブラリを使ってお手軽に係数が得られるっぽい。
http://web.archive.org/web/20140528164434/http://lolengine.net/wiki/doc/maths/remez

• • • • • -

しかしテイラー展開で精度的に問題が無いならそれで良い気がしてくる。ややこしい説明が不要だし。というか値域で精度が良い方を選択して使うのじゃ駄目なんだろうか？アンバランスってそこまで問題なのかな。。

double sine2(double x)
{
	bool bMinus = false;
	if (x < 0) {
		x = -x;
		bMinus = true;
	}
	x = fmod(x, 2 * PI);
	if (x > PI) {
		x -= PI;
		bMinus = !bMinus;
	}
	if (x > PI/2) {
		x = PI - x;
	}
#if 0
	double x2 = x * x;
	double x3 = x * x2;
	double x5 = x2 * x3;
	double x7 = x2 * x5;
	double x9 = x2 * x7;
	return
		x
		- x3/6
		+ x5/120
		- x7/5040
		+ x9/362880
	;
#else
	double x2 = x * x;
	double y = x * (1 + x2 * (-1.0/6.0 + x2 * (+1.0/120.0 + x2 * (-1.0/5040.0 + x2 * 1.0/362880))));
#endif
	return bMinus ? -y : y;
}