三角関数を自前実装してみる（前半）

はじめに

三角関数は、大体のプログラミング言語で標準的なライブラリに含まれているため、自前で実装することはあまりないでしょう。

普段何の考えなしに使っている関数ですが、三角関数をプログラム的に計算するにはどうやったらよいでしょうか。

ちゃんと調べればもっと効率や精度の良い方法が出てくるでしょうが、とりあえず今思いつく方法で実装してみます。

前提

Python3
あらかじめ十分な精度で円周率 $\pi$ を計算している
それっぽい値が計算できることを確認するだけで、結果の精度については考えない

アイデア

単位円上の点の座標を求める

単位円上の点 $p=(x, y)$ について、 $p$ と $x$ 軸正方向とがなす角を $\theta$ （単位はラジアン）とすると、

$\cos(\theta)=x, \sin(\theta)=y$

と表すことができます。

したがって、任意の $\theta$ に対して $x$ 軸正方向となす角が $\theta$ となるような単位円上の点 $p$ の座標を計算することが、三角関数を計算することになるのです。

f:id:ttk1:20190528015749p:plain

二分法

とはいえ、どうやってそのような点 $p$ を計算すればよいでしょうか。

二分法によって計算することができます。

$\theta=0$ に対応する点が $p=(1, 0)$ になることはわかります。 $\theta=\pi$ に対応する点が $p=(-1, 0)$ ということもわかります。

ある $\theta_1$ と $\theta_2$ について対応する点 $p_1=(x_1, y_1), p_2=(x_2, y_2)$ を知っているとき、 $\theta_3=(\theta_1+\theta_2)/2$ に対応する点を $p_3=(x'/k, y'/k)$ （ただし $x'=(x_1+x_2)/2, y'=(y_1+y_2)/2, k=||(x', y')||_2$ ）と表すことができます。

f:id:ttk1:20190528021633p:plain

したがって、 $\theta_1=0, \theta_2=\pi$ から始めて、 $\theta_3 - \theta$ が十分小さくなるまで、 $\theta_3 \lt \theta$ ならば $\theta_1 \gets \theta_3$ , $\theta_3 \geq \theta$ ならば $\theta_2 \gets \theta_3$ として繰り返し $\theta_3$ を計算することで目的の点 $p$ を得ることができます。

実装

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