ニュートン法による逆数の計算

【 逆数とは 】

 1~9までの整数の逆数について調べてみると...
 1の逆数は 1/1 = 1
 2の逆数は 1/2 = 0.5
 4の逆数は 1/4 = 0.25
 5の逆数は 1/5 = 0.2
 8の逆数は 1/8 = 0.125
これらの数値はその逆数も有限少数として求めることができる。しかし
 3の逆数は 1/3 = 0.333333...
 6の逆数は 1/6 = 0.166666...
 7の逆数は 1/7 = 0.142857...
 9の逆数は 1/9 = 0.111111...
これらは無限小数となるため求めることができない。

 ここでは任意の数値 a に対する逆数 1/a を高精度※1で求める。

 ※1 ここで言う「高精度」とは小数点以下15桁以上のことを指す。15桁程度を
  求めたければ...普通に変数をdoubleで定義して 1/a を計算すればいい。
  またWindows付属の電卓でもその程度は求められる


【 目的 】

 そもそも逆数を計算してどうするのか?

 ある任意の数値 a と b があるとき、 b * a という乗算と b / a という除算では、 乗算の方が計算量が少ない※2

 そのため b / a = b * (1/a) とし、「 b に対し a の逆数を乗算する」とする。こ こで使うための逆数を求める※3のが目的だ。


 ※2 どっちも1回じゃないか!乗算と除算って違いはあるけど、なんて突っ込みは禁止
 ※3 逆数を求めるための計算量が増えてんじゃん、ってのもなしね


【 漸化式の取得 】

 ニュートン法に適用するための方程式 f(x) とその1次導関数 f'(x) を求める。
 このとき、f(x) = 0 の解がaの逆数 1/a になるようにする。また無限に微分できるこ とが好ましい。


 次に上の方程式に対するニュートン法の漸化式を求める。

 これで逆数 1/a を求めるための漸化式として xn+1 = xn(2 - a * xn) が求まった。


【 初期値と収束条件 】

 ここで漸化式の初期値 x0 は求める逆数 1/a に近い値であればあるほど 収束が速いそのため漸化式の初期値は
          x0 = (double) 1/a ※4
とする。

 また求めた数値が収束したかどうかは 2 - a * xn の項がどれだけ1.00に 近づいたかを見ることで判定する。

 ※4 ここでの 1/a はdouble型の変数を使うため有効桁数16桁程度


【 計算例 】

 ここでは 123456789 の逆数 1/123456789 を計算してみた。
 有効桁数 (1800-9) 桁を 7 回のループで求めることができた。
 求めた数値の検算は同様の方法で求めた2500桁までの数値に対し 12345679 を乗算することで行い、誤差部分は赤く示した。

初期値:
0.000000008100000073710001300

1回目:
0.000000008100000073710000670761006103925106664026590000

2回目:
0.000000008100000073710000670761006103925155545718915466042130740983389742653856559247121258551761229100000000

3回目:
0.0000000081000000737100006707610061039251555457189154660421307409833897429488466608345046135939919837053270517184761706381331528070116905
33063277869297304365939249417034960140894927399673682220404746910000000

4回目:
0.0000000081000000737100006707610061039251555457189154660421307409833897429488466608345046135939919837053270517184761706381331528070116905
43806383948638093932606654786720558559157082888329454283798033982562109241315194095968266273311223087132130092902383845411692993232157556
174185826998993612371281756292720486968179344424947852272492244479407878399986975061066005930572402184167032286582040258827949100000


5回目:
0.0000000081000000737100006707610061039251555457189154660421307409833897429488466608345046135939919837053270517184761706381331528070116905
438063839486380939326066547867205585591570828883294542837980339825621092413151940959682662733112230871321300929023838454116929932464062385
422967707349006136875955845571198194697903571750922502933394776693892467914421458021235267993240938738492542520282136934567446104563759551
530211918924928462216849006173325956177266201213122431039414122458668514373883480802339675301291045241748511699660808386009329606009987014
143478809074187844714000491183416435881061299444977925164630313578509934408348200114824921218210099155333722308122477914827991571636880115
51285208307607815850317627963377701121233501237228895360046553758183966026431320875626530303637583697510042090910


6回目:
0.0000000081000000737100006707610061039251555457189154660421307409833897429488466608345046135939919837053270517184761706381331528070116905
438063839486380939326066547867205585591570828883294542837980339825621092413151940959682662733112230871321300929023838454116929932464062385
422967707349006136875955845571198194697903571750922502933394776693892467914421458021235267993240938738492542520282136934567446104563759551
530211918924928462216849006173325956177266201213122431039414122458668514373883480802339675301291045241748511699911456469194253869667710213
976162947183082819366053656231088271702903272496419779717419995428521958399549821435903375066720713107158489275142252403794496874529921558
222286179822804236387518551126418815250411218778742090886553027067632546315456171470651160382925559484622591310065580921596786386530756117
429880668611914084368418167752605326548708471593247091498548532636791646994803987652716287639718217521435779445065592950096895845881752197
523944997467899476957885240316755686882476750630538430737899719714887449404700178959797906096990028925667496102833354906178189432783902205
458137138883843917510794160533364055622731458999944493917872808881682519493387663990827197582189260053323386937274760415631189890665408423
750506435541324302977816087218314849415674535225152615161821232559459292564838862789023857365598306864073718337837000122148735181428882485
101783605704934933164770041641837883507037930482822634779656823885495068410045629597415916071562544642489832826747471894891075403121233248
57347506532060449083781084559976209045363087872855652010485806539160453707961334829100


7回目:
0.0000000081000000737100006707610061039251555457189154660421307409833897429488466608345046135939919837053270517184761706381331528070116905
438063839486380939326066547867205585591570828883294542837980339825621092413151940959682662733112230871321300929023838454116929932464062385
422967707349006136875955845571198194697903571750922502933394776693892467914421458021235267993240938738492542520282136934567446104563759551
530211918924928462216849006173325956177266201213122431039414122458668514373883480802339675301291045241748511699911456469194253869667710213
976162947183082819366053656231088271702903272496419779717419995428521958399549821435903375066720713107158489275142252403794496874529921558
222286179822804236387518551126418815250411218778742090886553027067632546315456171470651160382925559484622591310065580921596786386530756117
429880668611914084368418167752605326548708471593247091498548532636791646994803987652716287639718217521435779445065592950096895845881752197
523944997467899476957885240316755686882476750630538430737899719714887449405475789589829685267450135933796236997545756677666385766764110477
553405345735988646197496680397219791614700103693770943613315586881171840618663749629840121631545106847060472308250298005077711846207177800
485317984416393658189182289521558834646185395280287097050612583160574506761228011527174904897291634565353874544720258356954351048284594539
389810308447273806870191642518743946920569916977186244492394824880792906415215448378460580243991280220320650004917915044753026907252544855
998158189583239525207479679388065082431392250125669476143592232906689319450872807002942543726777147913672046014415618731182130453757387129
192222875649228168407976332512584625864520095367132867840909097352272785905682351741709400849555547730955484351694907600423659163855298391
08321535


【 さらなる収束のために 】

 この方法は2次収束をする。
 3次、4次と高次に収束する方法については...そのうちできるかなぁ。


カテゴリー「VC++ TIPS」 のエントリー