2005/11/01 数値解析学と演習 福井 先生 関数近似 初等関数をどう計算するか ? ( 計算量は、四則の 50 倍程度かかった.. ) 単純には.. テーラー展開 多項式 / 連分数 しかし、テーラー展開は、h が 0 の近くでは誤差が小さいが、 h が大きくなると、誤差が大きくなる。 望ましい誤差曲線 近似区間の誤差の最大が、最小になるように係数を調整する min-max 近似 ( 最良近似 ) が望ましい。 この為には、誤差カーブが決らないといけない しかし、誤差カーブが解るということは、真値が解るということ これは、矛盾 : 真値が解らないから、関数近似をしている.. 例: 24 bit の真値を計算するには ( まるめ誤差のことを考えると.. ) 倍精度の 64 bit で計算し、テーラーの項目をふやせば、時間がかかるが、ある程度求めることができる。 じゃあ、64bit は ? 128bit ? じゃあ 128bit は.. ? ハードウェアの bit 数は、限界がある ソフトウェアでの多倍長ルーチンを作っている 真値がわかっても.. 与えられた真値と、近似関数から、min-max を計算するのも 非線形な計算なので大変 # 関数近似の大家、日大の故山下先生、小野先生、多田先生、浜田先生 三角関数の近似 sin が理解れば、cos も解る。 # できるだけ計算をサボりたいので.. 周期関数なので [0,2π] 形が対称なので [0,π] 同じ対称なので [0,π/2] 加法定理を使えば、更に π/4, π/6 まで狭めることができる。 # 人の好き嫌いでπ/4, π/6を選ぶ # π/4 : 計算する区間は広いが、折畳む手間が少い # π/6 : 計算する区間は狭いが、折畳む手間が多い ## 計算量的には、どっちもどっち \sin{x} のテーラ展開 f(x)=\sin{x} f(x+h)= f(x)+\frac{h}{1!}f'(x)+\frac{h^2}{2!}f''(x)+\frac{h^3}{3!}f'''(x)+\frac{h^4}{4!}f^{(4)}(x) + .. より、 \sin{x+h}= \sin{x}+\frac{h}{1!}\cos{x}-\frac{h^2}{2!}\sin{x}-\frac{h^3}{3!}\cos{x}+\frac{h^4}{4!}\sin{x} + .. ここで x = 0 を代入すれば、 \sin{h}= \frac{x}{1!}-\frac{x^3}{3!}+\frac{h^5}{5!}-\frac{h^7}{7!} .. これが、近似式 0 <= h <= π/4 < 1 だし、階乗の計算は速く大きくなるので収束は速い さらに、これを min-max にすると... \sin{h}=h-0.1.666h^3 + 0.083h^5 .. のようになる。 # ここらへんの係数の調整は、それだけで、ノウハウが沢山ある !! 引数の誤差の問題 周期関数なので、 x の計算をするときには、 x mod 2π を計算すれば 良いが.. x の値がすごく大きいと、 x mod 2π の精度が悪くなってしまう.. => 桁落が生れる # 名古屋の二宮先生 # x を a + b x 2π の形で、計算する !! # => a の所に誤差は出ない。 πの値に含まれる誤差も問題 !! # 計算機の中で、何を使っているかは良くわからない !! # \arctan{1} で π/2 を計算する.. 位相の誤差がしょうじている # 丁度、符号が変るあたりが、怪しい。 # ここらへんは、sin 関数固有の問題 # tan も別な点で気を付ける必要がある sqrt は、多項式近似でもよいが、ニュートン法が速くてよい。 f(x) = x^2 - a とし、 f(x) = 0 の解αを求めれば、 α = \sqrt{a} となる。 ニュートン法 f'(x) = 2x より、 x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} なので、 x_{k+1} = x_k - \frac{(x_k^2-a)}{2x_k} = \frac{1}{2}( x_k + \frac{a}{x_k} ) 初期値の x_0 が問題だが、一般には、テーブルをもっている。 単精度の値がでれば、もう１回で、倍精度の精度が出る # cf 二次収束 a が 0 に近いと収束が速い ( 大きいと遅くなる )。 a が大きい時は、 a = A * 4^n 1 \le A \lt 4 とすれば、 \sqrt{a} = \sqrt{A*4^n} = \sqrt{A}*\sqrt{4^n} = \sqrt{A}*2^n となる。 こうすれば、x_0 のテーブルも小くできる。 4 でわる、2 をかけるは、非常に単純な計算 指数部を調整するだけ。 cf. 一松信、「関数近似」教育出版 培風館の「関数近似」の三冊 ( 最後に、定数表があるので厚い ) ハーツの本 # 関数近似を行っている先生は少い == 関数近似は、計算機の前から研究されている 30 年前に、関数近似に関する衝撃があった。 HP 35 : プログラム電卓 50 万円 : 一松先生 日本で最初 30 万円 : 平野先生 多項式近似ではない近似 計算機が得意な、論理演算、シフト等で関数計算をする。 多項式近似は、浮動小数点計算ができることを前提としている 浮動小数点計算ではなく、論理演算やシフトで計算したい => CORDIC x を二進数展開し、各桁を x_i で表す x_i = 0 / 1 x = \sum{i= 0}^n{x_i 2^{-i}} sin の加法公式を使うと.. \sin{A+B} = \sin{A}\cos{B}+\cos{B}\sin{A} \sin{x} = \sin{\sum{i= 0}^n{x_i 2^{-i}}} = ... 展開した先の値にかんしては関数表をを作る 関数表は、まばらなので、全部記憶してなくて、シフトで計算すればよいので、それほど大きな表にならない シフトと和で計算ので掛け算と良くにた計算でできる。 # 手計算の時代には、考えられなかった手法 # 仮におもいついても、実行できなかっただろう.. STL ( Successiable, Table, Lookup ) Intel の Chip の関数表に誤りがあった CODIC のよい点 bit 数がふえても、loop 数が増えるだけで実現できる == 折畳み 計算する空間を縮小する手法 関数の性質を巧く利用して行う 初等関数はだいたい、性質がよい \tan だけはだめ 結果が無限に飛ぶのは性質が悪い \log は、 \log{x*2^n} = \log{x} + n\log{2} なので、なんとかなる。 == 電卓で遊んでみよう sqrt 1 より大きな数の sqrt は、1 になる。 1 より小さい数の sqrt は、1 にならない.. 0.99998 位まで.. sin/cos はともかく tan は結構変な値ができる == [今日の演習] # Text にはのっていない。 \sin{x} のテーラ展開を行い関数近似を行う。 x^7 の項目までの、展開した結果は、 \sin{x}= \frac{x}{1!}-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!} なので、これを近似式としてみる。 真値は、ライブラリの sin 関数で求め、計算結果との間の誤差を出力する。 x の範囲は、0 から 1.0 まで、0.02 刻みで行ってみる。