はじめに
三角関数は、大体のプログラミング言語で標準的なライブラリに含まれているため、自前で実装することはあまりないでしょう。
普段何の考えなしに使っている関数ですが、三角関数をプログラム的に計算するにはどうやったらよいでしょうか。
ちゃんと調べればもっと効率や精度の良い方法が出てくるでしょうが、とりあえず今思いつく方法で実装してみます。
前提
- Python3
- あらかじめ十分な精度で円周率
を計算している
- それっぽい値が計算できることを確認するだけで、結果の精度については考えない
アイデア
単位円上の点の座標を求める
単位円上の点 について、
と
軸正方向とがなす角を
(単位はラジアン)とすると、
と表すことができます。
したがって、任意の に対して
軸正方向となす角が
となるような単位円上の点
の座標を計算することが、
三角関数を計算することになるのです。
二分法
とはいえ、どうやってそのような点 を計算すればよいでしょうか。
二分法によって計算することができます。
に対応する点が
になることはわかります。
に対応する点が
ということもわかります。
ある と
について対応する点
を知っているとき、
に対応する点を
(ただし
)と表すことができます。
したがって、 から始めて、
が十分小さくなるまで、
ならば
,
ならば
として繰り返し
を計算することで目的の点
を得ることができます。
実装
import math pi = 3.141592653589793 def g(p1, p2): if (p1[0] + p2[0]) == 0: return [0.0, 1.0] p = [(p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2] k = math.sqrt(p[0] ** 2 + p[1] ** 2) return [p[0] / k, p[1] / k] def f(t): t1 = 0 t2 = pi p1 = [1.0, 0.0] p2 = [-1.0, 0.0] while t2 - t1 > 0.000000000000001: t3 = (t1 + t2) / 2 if t3 < t: t1 = t3 p1 = g(p1, p2) else: t2 = t3 p2 = g(p1, p2) return p1 def sin(t): return f(t)[1] def cos(t): return f(t)[0] def tan(t): p = f(t) return p[1] / p[0]
実行結果
動作確認
if __name__ == '__main__': print(sin(1), math.sin(1)) print(cos(1), math.cos(1)) # 結果 0.8414709848078963 0.8414709848078965 0.5403023058681402 0.5403023058681398
math.sin/math.cosと比較して、それっぽい値になっていることがわかります。
ベンチマーク
Benchmarker を使って実行時間を計測してみます。
from benchmarker import Benchmarker import numpy as np import random import math if __name__ == '__main__': random.seed(0) loop = [random.uniform(0, pi) for _ in range(100 * 1000)] with Benchmarker(width=20) as bench: @bench('tama') def _(bm): for r in loop: sin(r) @bench('numpy') def _(bm): for r in loop: np.sin(r) @bench('math') def _(bm): for r in loop: math.sin(r) # 結果 ## real (total = user + sys) tama 5.0977 5.1094 5.1094 0.0000 numpy 0.0560 0.0469 0.0469 0.0000 math 0.0100 0.0156 0.0156 0.0000 ## Ranking real math 0.0100 (100.0) ******************** numpy 0.0560 ( 17.8) **** tama 5.0977 ( 0.2) ## Matrix real [01] [02] [03] [01] math 0.0100 100.0 561.5 51094.4 [02] numpy 0.0560 17.8 100.0 9099.3 [03] tama 5.0977 0.2 1.1 100.0
パフォーマンスはmathやnumpyに遠く及ばない...厳しい...
まとめ
- 自力で何とか出来た
- ただしパフォーマンスは残念なことに...
圧倒的に差が出てしまいましたが、mathやnumpyの三角関数はどういった方法で実装されているのでしょうか。 時間があるときに調べてみようと思います(後編へ続く)。