平方根( sqrt )の数値的解法③
既に平方根の計算方法として最も有用なニュートン法を紹介していますが、その他にも興味深い方法がありますので紹介しておきます。
＜テーラー展開＞
平方根の計算方法としては難がありますが、計算機科学はテーラー展開の発見なしには語れません。\( sin(x), cos(x), exp(x), log(x)\) など、平方根以外の初等関数はそのほとんどがテーラー展開によるものですし、力学シミュレーションの漸化式を与えているのもテーラー展開です。
平方根の計算においては次の式の形に直します。(０近傍で無限回微分が可能な形)
\begin{align*} \sqrt{\alpha} = f(x) = (x+1)^{\frac{1}{2}} \tag{1}\end{align*}
(１)を0近傍で展開(マクローリン展開)してみましょう。
\begin{align*} f(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(0)}{n!}x^n = 1 + \frac{1}{2}x-\frac{1}{8}x^2+\frac{1}{16}x^3 - \frac{5}{128}x^4+\dots \end{align*}
係数の一般化をしておくと
\begin{align*} \binom{z}{m}:=\frac{z(z-1)(z-2)\dots(z-m+1)}{m!} \end{align*}
と定義して
\begin{align*} f(x) = 1 + \binom{1/2}{1}x+\binom{1/2}{2}x^2+\binom{1/2}{3}x^3 + \binom{1/2}{4}x^4+\dots \tag{2}\end{align*}
と書けます。この形はテーラー展開が発見される以前からニュートンの一般二項定理として知られた形でした。\(x\)が０に近いほど収束が早く\(x < 1\)である必要があります。ですので残念なことに例えば \(\sqrt{1+1}\)として２の平方根を取ろうとすると収束が悪いです。この制約からテーラー展開による平方根は入力値によって形を変えてやる必要があります。例えば、３、５、６、７の平方根であれば、
\begin{align*} \sqrt{3} = 2(1-\frac{1}{4})^{\frac{1}{2}} \end{align*}
\begin{align*} \sqrt{5} = 2(1+\frac{1}{4})^{\frac{1}{2}} \end{align*}
\begin{align*} \sqrt{6} = 2(1+\frac{1}{2})^{\frac{1}{2}} \end{align*}
\begin{align*} \sqrt{7} = 2(1+\frac{3}{4})^{\frac{1}{2}} \end{align*}
といったように式を直します。(実際はもっと分数部分を小さくする工夫ができます) これらの式の手直しを数値計算で対応するのは難しく、テーラー展開は平方根の計算には向いていません。では代表として\(\sqrt{2}\)のプログラムを見てみましょう。係数の性質からホーナー法を用いています。
#include "stdafx.h"
#include <iostream>
double sqrt2(void);
int _tmain(int argc, _TCHAR* argv[ ])
{
std::cout << "sqrt(2) =" << sqrt2() << std::endl;
return 0;
}
double sqrt2(void)
{
double x1=0 ,x2, a;
int n;
a = 1.0;
x2 = a ;
for(n=1; ( (a>=0)? a: -a) > 0.000001; n++){
x1 = x2;
a = a * ( (1.5 -(double)n )/(double)n);
x2 = x1 + a ;
}
std::cout << "n =" << n-1 << std::endl;
return x2;
}
このプログラムの結果は下記になります。
n =4302 sqrt(2) =1.41421
収束が悪いですが、テーラー展開の問題ではなく初期値として\(x=1\)としてしまっているからです。テーラー展開は収束半径に気を付けて使う必要があります。
＜連分数＞
連分数は通常の学習指導要領には含まれないかと思いますので学校では習はなかったかと思いますが、分数では表せないはずの無理数を分数の形で表現する事が出来ます。黄金比を表現する方法として知られています。数値計算では割り算が続くので計算量が多く好まれないのかと思われます。
さてここでは例として\(\sqrt{2}\)と\(\sqrt{5}\)を連分数の形で表して一般化してみます。
\begin{align*} \sqrt{2} = 1 + \sqrt{2} -1 \tag{1}\end{align*}
と表せます。外に出した1 は\( a^2 \le \alpha: \alpha = 2 \) である最大の\(a\)を選びます。ここで\(\sqrt{2}-1\)の項に対して乗法の公式を使って有理化っぽい事(本来は分母を有理数にする意味)をしてみます(分母に無理数を作り出すので無理化かもしれません)。
\begin{align*} \sqrt{2}-1 = \frac{\sqrt{2}+1}{\sqrt{2}+1} \times (\sqrt{2} -1) = \frac{1}{\sqrt{2} + 1} \end{align*}
となります。式①はこれにより次のように表せます。
\begin{align*} \sqrt{2} = 1 + \frac{1}{\sqrt{2} + 1} \tag{2}\end{align*}
ここで右辺の\(\sqrt{2}\) にも式②を代入し繰り返してみましょう。(再起呼び出しと同じこと)
\begin{align*} \sqrt{2} = 1 + \frac{1}{ 2 + \frac{1}{ 2 + \frac{1}{2 + \frac{1}{2 + \dots}}}} \tag{3}\end{align*}
このように連分数として表現する事が出来ました。無限に分母が続いてしまいますが、それでは計算できませんので適当なところで\(\sqrt{2}\)には１を入れます。 同様に\(\sqrt{5}\)でもやってみましょう。
\begin{align*} \sqrt{5} = 2 + \sqrt{5} -2 = 2 + \frac{3}{ 4 + \frac{3}{ 4 + \frac{3}{ 4 + \frac{3}{ 4 + \dots}}}} \tag{4}\end{align*}
こちらの場合の最終的な分母の\(\sqrt{5}\)には\( a^2 \le \alpha: \alpha = 5 \) である最大の\(a\)である２を入れます。 通常連分数では分子に１を持ってくる正則連分数が好まれるようですが、ここではこの形でOKとさせて下さい。さて、それでは代数的に一般化した形で表現してみたいと思います。
\begin{align*} \sqrt{\alpha} = a + \frac{\alpha - a^2 }{ 2a + \frac{\alpha - a^2}{ 2a + \frac{\alpha - a^2}{2a + \frac{\alpha - a^2}{2a + \frac{\alpha - a^2}{a+a}}}}} \tag{5}\end{align*}
\(\alpha\)が入力値、\(a\)は\( a^2 \le \alpha\)である最大の\(a\)を選びます。数値計算においてはこのaを選ぶのに直観というわけには行かないので、iterationで探す手間があります。それではプログラムを見てみましょう。
#include <iostream>
double sqrt(double alpha);
int _tmain(int argc, _TCHAR* argv[])
{
std::cout << "sqrt(2)=" << sqrt(2.0) << std::endl;
return 0;
}
double sqrt( double alpha ){
double numerator, x1=0, x2;
int a=0, n;
while(alpha > (a+1)*(a+1) ) a++;
numerator = alpha - (double)(a*a);
x2 = numerator / (double)(a+a);
for(n=0;x1!=x2;++n){
x1 = x2;
x2 = numerator / ( (double)2.0 * a + x1 );
}
x2 = x2 + (double)a;
std::cout<<"n="<< n << std::endl;
return(x2);
}
このプログラムの結果は下記になります。
n=21 sqrt(2)=1.41421
doubleで計算できるマックスまでiterationさせて２１回で収束しているので、二分法より早く収束しています。決して悪くない方法に思えます。
＜開平法＞
ルートを取る事を開平するといったり、ルートを開くと言ったりします。昔は学校でも必ず習う計算手法だったみたいでルートの計算といえば開平法という印象があります。幾何的な特徴からルートを求めます。基本的には筆算用の方法なので、プログラム化すると複雑化してしまいわざわざ開平法を再現するプログラムというのは見かけません(2進だとすっきりしたプログラムになりますがここでは10進でのプログラムを紹介します)。まず前提として10進の数は分解すると次のようになります。
\begin{align*} r = r_n 10^n + r_{n-1} 10^{n-1} + \dots + r_0 10^0 + r_{-1} 10^{-1}\dots = \sum_{k}^{n}r_k 10^k\tag{1}\end{align*}
\(n\)は整数、\(r_k\)は\(0 \le r_k \le 9 \)である整数です。何進数であろうとこの考え方は同じです。もう一つ前提としてルートを取りたい数\(\alpha\)が\(10^{2n} \le \alpha \le 10^ {2(n-1)}\) とすると、\( (10^{i})^2 \le \alpha \)である最大の\(i\)が\(n\)となります。
さて、では例として144の開平を見ていきます。144は\( 10^{2 \times 1} < 144 < 10^{2 \times (1+1)} \)なので開平後の最大の\(n\)は 1 となります。開平後を①式のように表してみると、
\begin{align*} \sqrt{144} = r_1 10^1 + r_0 10^0 \tag{2}\end{align*}
ルートの解として知りたいのが\(r_1\)と\(r_0\)という問題にすり替わります。②式を二乗すると、
\begin{align*} 144 = r_1^2 10^2 + 20 r_1 r_0 + r_0^2 \tag{3}\end{align*}
です。二乗というのは面積計算なので\(144＝S_k\)と置くと③式は\(S_k = S_1 + S_0 \)
\begin{align*} S_1 = r_1^2 10^2 \tag{4}\end{align*}
\begin{align*} S_0 = 20 r_1 r_0 + r_0^2 \tag{5}\end{align*}
となります。この後は段階的に\(r\)を検算して求めます。結局人間の手の終えるレベルにして検算で求めるという方法なのですが、試しに\(r_1\) に2を入れると、S1=200になって144より大きくなってしまうので\(r_1\)は１になります。とすると、\(S_0=44\)となりますので、⑤式を書き直すと、
\begin{align*} 44 = 20 r_0 + r_0^2 \tag{6} \end{align*}
です。\(r_0\)に合いそうな数を類推して当てはめます。３だと大きくなりすぎるので２です。よって答えは１２となります。開平法というと割り算みたいな計算を思い出しますが、実際やっている事はこの手続きになります。
それではプログラムを見てみましょう。
#include <iostream>
#include <deque>
#include <cmath>
double kaihei(double alpha);
int _tmain(int argc, _TCHAR* argv[ ])
{
std::cout << kaihei(2) << std::endl;
return 0;
}
double kaihei(double alpha)
{
int i, j, k, t1, t2;
std::deque<unsigned char> v_p;
std::deque<unsigned char> v_a;
double ans = 0;
i=0;
while( pow(10, 2*(i+1)) <= alpha ) i++;
j=0;
while(alpha != 0 ){
v_p.push_back( (unsigned char)(alpha / pow(10, 2*(i-j))));
alpha = alpha - v_p[j]*pow(10, 2*(i-j));
j++;
}
k=0;
while( (k+1)*(k+1) <= v_p[0]) k++;
v_a.push_back( k );
t2 = ( v_p[0] - v_a[0] * v_a[0] ) * 100;
if( j-1 > 0 ) t2 = t2 + v_p[1];
t1 = v_a[0] + v_a[0];
for(int l=1; l<10; l++){
k=0;
while( ( ( t1 * 10 ) + (k+1) ) * (k+1) <= t2 ) k++;
t2 = ( t2 - ( ( t1 * 10 ) + k ) * k) * 100;
if( j-1 > l ) t2 = t2 + v_p[l];
v_a.push_back( k );
t1 = ( ( t1 * 10 ) + k ) + k;
}
for(int m=0; m<10; m++) ans = ans + (double)( v_a[m] * pow(10, i-m) );
return (ans);
}
このプログラムの結果は下記になります。
sqrt(2)=1.41421
cmathをincludeしてしまっていますし、アロケートして配列を増やしています。さらに今までとは毛色が違いiterationの数を10と決めています。iterationするだけ下位の桁まで計算できます。基本的には標準関数には向かないと思いますが頭の体操と思ってプログラム内容を考えてみてください。
まとめ
５種類の方法を見てきたが平方の計算にはニュートン法が向いている。
