記憶力が無い

プログラミングとランニングとカメラと何か

三角関数を自前実装してみる(前半)

はじめに

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

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

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

前提

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

アイデア

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

単位円上の点 p=(x, y) について、 px 軸正方向とがなす角を \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

パフォーマンスはmathやnumpyに遠く及ばない...厳しい...

まとめ

  • 自力で何とか出来た
  • ただしパフォーマンスは残念なことに...

圧倒的に差が出てしまいましたが、mathやnumpyの三角関数はどういった方法で実装されているのでしょうか。 時間があるときに調べてみようと思います(後編へ続く)。

Copyright © 2017 ttk1