マクローリン展開による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=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 = i / 12 y1, y2 = truefunc(s), myfunc(s) d = abs(y1 - y2) try: dl = math.log10(d) - math.log10(y1) except ValueError: dl = 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__ == "__main__": main()
ついでに、コサイン・指数関数・アークサインも作ってみました。
iter_***_coefficientsを差し替えれば、良い様になっています。
それぞれ、mathモジュールのやつと値を比較し、
誤差と、(誤差 / 値)の対数を計算して表示しています。
そして最後に、timeitで、math.sinとパフォーマンスを比較しています。
まあ、math.sinはCで書かれてるんでしょうから、こんなもんですかね?
sin | 20.9441755866 |
math.sin | 0.378902346527 |