マクローリン展開によるsin関数(Python)



マクローリン展開によるsin関数

に触発されて、Pythonでも同じ事をやってみました。


・・・と言っても、Haskellの直訳はしません（できません）。




書いてみたら、イテレータ・ジェネレータだらけに・・・。



#sin.py
from __future__ import (
    with_statement, 
    division,
    print_function,
)

from itertools import (
    count,
    takewhile,
    ifilter
)
from functools import (
    reduce,
    partial
)
from operator import mul
from fractions import Fraction

def fact(n):
    return reduce(mul, xrange(1, n + 1), 1)
    
def is_odd(n):
    return bool(int(n) % 2)
    
def iter_sin_coefficients():
    for n in count(0):
        yield 2 * n + 1, (-1) ** n / fact(2 * n + 1)

def iter_cos_coefficients():
    for n in count(0):
        yield 2 * n, (-1) ** n / fact(2 * n)

def iter_asin_coefficients():
    for n in count(0):
        yield 2 * n + 1, fact(2 * n) / ((4 ** n) * fact(n) ** 2 * (2 * n + 1))

def iter_exp_coefficients():
    for n in count(0):
        yield n, 1 / fact(n)
    
def iter_term(x, icoefficients):
    for n, a in (icoefficients):
        yield a * (x ** n)
    
def takewhile_greater_than(epsilon, iterable):
    epsilon = abs(epsilon)
    assert epsilon > 0
    
    def abs_is_greater_than_epsilon(y):
        return epsilon < abs(y)
    
    return takewhile(abs_is_greater_than_epsilon, iterable) 

def maclaurin(x, icoefficients, epsilon&#061;10 ** -4):
    return sum(
        takewhile_greater_than(
            epsilon, 
            iter_term(
                x, 
                icoefficients)))

def sin(x):
    return maclaurin(x, iter_sin_coefficients())

def cos(x):
    return maclaurin(x, iter_cos_coefficients())

def exp(x):
    return maclaurin(x, iter_exp_coefficients())

def asin(x):
    return maclaurin(x, iter_asin_coefficients())

    
def main():
    import math
    for myfunc, truefunc in [(sin, math.sin), (cos, math.cos), 
                             (exp, math.exp), (asin, math.asin)]:
        print(myfunc.__name__)
        for i in xrange(1, 12 + 1):
            s &#061; i / 12
            
            y1, y2 &#061; truefunc(s), myfunc(s)
            d &#061; abs(y1 - y2)
            try:
                dl &#061; math.log10(d) - math.log10(y1)
            except ValueError:
                dl &#061; 0
            
            print("%e, %e, %e, %e"%(y1, y2, d, dl))
        print()
    
    import timeit
    print("sin",    timeit.Timer(partial(sin, math.pi / 6)).timeit())
    print("math.sin", timeit.Timer(partial(math.sin, math.pi / 6)).timeit())
    
    
if __name__ &#061;&#061; "__main__":
    main()


ついでに、コサイン・指数関数・アークサインも作ってみました。

iter_***_coefficientsを差し替えれば、良い様になっています。



それぞれ、mathモジュールのやつと値を比較し、
誤差と、(誤差 / 値)の対数を計算して表示しています。



そして最後に、timeitで、math.sinとパフォーマンスを比較しています。

まあ、math.sinはCで書かれてるんでしょうから、こんなもんですかね？













sin 20.9441755866
math.sin 0.378902346527