ネムネコ、桁落ちに会う！！

数日後の記事に使おうと、Fotranでパタパタとプログラムを作ってみた。
関数f(x)の微分f'(x)を使う箇所があるので、「微分なんていちいちしてられるか！！」と思い、この部分は中心差分を使って、微分係数f'(x)を求めることにした。
そうしたら、真面目に微分を計算したものと計算結果が違いやがった！！

　　

中心差分というのは、f(x)の微分を上のように計算する方法。h=1.0×10⁻⁶くらいにとれば、真面目に微分を計算したものとほとんど一致するはずだ（註を参照）と考えたのだが――打切誤差という観点からすると、この予想はまったく正しい――、これが合わないんだケロ。h=1.0×10⁻⁴あたりで精度がもっとも良くなり、それ以上hを小さくとると、精度が逆に落ちてくるんだケロ。

「何故だろう？」としばし考えたところ、どうやら、桁落ちや情報落ち、積み残し――以下、これらを総称して桁落ちと呼ぶことにする――が発生し、それで精度が落ちているみたいだね。

かれこれ１年ほど前に、数値計算のところで、丸め誤差、打切誤差、桁落ちなどを記事に書き、その時、桁落ちを発生させる計算例を示したけれど、あれは意図的に作り出した計算例で、ネムネコは、これまでに、自分の計算で、桁落ちで計算結果が大きく異なるという目に会ったことがなかった。それだけに、この計算結果を目の当たりにした際の衝撃は大きく、「何故（なぜ）に・・・」と、しばらく思考停止状態に陥ってしまうほどであった。
プログラムのどこかが間違っているというのならば、ソースコードを見なおせばすぐにわかるけれど、プログラム自体はどこも間違っていないんだから、そりゃ〜、すぐに、この不可解な現象の原因特定ができないケロ。

真面目に微分を計算した計算結果

微分を中心差分で計算した結果

計算結果が４．１１１１１ホニャララになるところが約４．０６８になっており、計算結果が１％以上も異なってしまう。微分はもっとひどくて、なんと、最大約１０％の誤差、アンダー・エスティメイトになるんだよね〜。しかも皮肉なもので、微分を真面目に計算した奴は２０回の反復で収束しないのに対して、微分を中心差分で代用したものは反復回数１８回で収束する――微分係数を過小評価するので、その結果、修正量が大きくなり、計算が早く収束する！！――。微分計算で桁落ちさせ、この部分の有効精度を落とした方が収束が早いのだから、もう笑うしかない。
果たして、どちらが、数値解法として優れているのやら（苦笑い）

う・ふ・ふ・ふ　毎日だれかに
う・ふ・ふ・ふ　見られることが
ビタミンになる

ではあるが、ネムネコの場合、これだけでは足りないにゃ。

これくらいしないと・・・。

さて、この事象の根本的解決策は、単精度の計算を倍精度にする以外にないね。計算が単純すぎるので、計算方法の工夫によって桁落ちを解消するということは、どうやら、できそうにない。
プログラムを作るのが楽なので、Fortranでプログラムを作ったのだけれど、使用する変数の全ての型宣言をするのならば、C言語の方が楽だな。しかも、精度的にも、FotranよりCの方が優れているし。
３０行に満たないオモチャみたいなプログラムだから、最初から作りなおしたところで大変じゃないけれど――プログラムを書き、コンパイル、バグ取り、計算を実行させてたとしても、この全作業は５分以内に終わるはず。１０分なんて時間は、絶対に、かからない――、プログラムを作った時間が無駄になってしまったケロよ。

でも、実際、「桁落ちで本当に精度が落ちてしまう」ということを知ることができて良い経験になった。これまで、そんなことに遭遇したことがなかったから。

微分を中心差分で計算した結果（C言語で倍精度計算）

表計算ソフトで計算した結果（真面目に微分計算）

プログラムで解いたものは、３列目が修正量、対して、表計算ソフトの３列目は推測値と解との誤差で違うのだけれど、倍精度で計算すれば、桁落ちが解消されていることがわかる。

ここしばらくC言語なんて使っていなかったから、「C言語の出力コマンドprintfの書式ってどうだったっけ」と、一瞬、考えこんでしまったケロよ。

C言語は、UNIXというOSを開発するために作られた、限りなくアセンブラに近いプログラミング言語だから、入出力系のコマンドは弱いね〜。

【補足】

情報落ち、積み残しの例　通常のコンピュータの計算（浮動小数点を用いた計算）では、百億＋１＝百億となるので、１を加えたという情報が落ちてしまう、消失する。このように、絶対値の大きな数と絶対値の数を加えた時、絶対値の小さな数が計算結果に反映されない現象を情報落ちや積み残しなどという。
金銭計算、大きな整数同士の四則演算に特化したプログラム言語、COBOLが現在も使われ続けるのは、整数同士の演算に限っては、このような情報落ち、積み残しが起きにくいため。COBOLならば、１０００兆＋１＝１０００兆１と、正しい値を出してくれる。お金が絡む計算だから、１円でも計算結果があわないと、大変なことになるからね〜。

桁落ちの例　たとえば、近い数字同士で引き算をすると、計算の有効精度が大きく落ちる、このような現象を桁落ちという。100.0001−100.0000=0.0001　７桁の精度をもっていたものが、この引き算を行うと、精度はなんと１桁に落ちてしまう！！
　――説明を容易にするために、ここでは、７桁、１桁の精度と書いた。――
それだけで済めばいいが・・・

単精度のFortranの場合、さらに、丸め誤差の影響がこれに加わり、100.0001−100.0000という引き算の結果は上のようになってしまうのであった。

【註】

　　

h=0.000001として、ネムネコのPCにある電卓で計算すると、

　　

となるんですがね〜。打切誤差が高精度な電卓の計算の丸め誤差の範囲に収まり、厳密な３６とドンピシャ一致する。

しかし、Fortranなどのプログラミング言語でこの計算をすると、この値は３６とかなり違う値が出るんだケロ。「誤差はあっても、せいぜい、0.1％程度だろう」という、ネムネコの甘い期待は木っ端微塵に打ち砕かれてしまった。

それにしても、電卓って、ホント、高精度だよな〜。