Xn+1 = 2 * Xn - A * Xn2で与えられる。有効桁は、繰り返しごとに2倍になる。多桁の計算では乗算の計算時間が支配的であるので、目的とする桁数の計算結果を得るための乗算回数を少なくするために高次の収束をする漸化式を利用するのが望ましい。 ここでは 1 / (1 - h) の級数展開を使用している。
hn = 1 - A * xn1 / A の4次収束の漸化式
xn+1 = xn * (1 + hn + hn2)
繰り返しごとに有効桁は3倍になる。hn は、xn の誤差と同じである。
hn = 1 - A * xn1 / A の5次収束の漸化式
xn+1 = xn * (1 + hn) * (1 + hn2)
hn = 1 - A * xn
xn+1 = xn * (1 + (1 + hn2) * (hn + hn2))
Xn+1 = (Xn + A / Xn) / 2で与えられるが、この式には A / Xn の除算を行わなければならないという不利な点がある。
Xn+1 = Xn * (3 - A * Xn2) / 2つまり、上式は(3/2, 1/2 は定数であるので)、乗算が必要とされるだけである。
hn = 1 - A * xn21 / √A の 4次収束の漸化式
xn+1 = xn * (1 + hn * (4 + 3 * hn) / 8)
hn は、xn の誤差と同じである。繰り返しごとに有効桁は3倍になる。
hn = 1 - A * xn21 / √A の 5次収束の漸化式
xn+1 = xn * (1 + hn * (8 + hn * (6 + 5 * hn)) / 16)
hn = 1 - A * xn21 / √A の 6次収束の漸化式
xn+1 = xn * (1 + (64 * hn + hn2 * (48 + 40 * hn + 35 * hn2)) / 128)
hn = 1 - A * xn2
xn+1 = xn * (1 + hn * (1/2 + hn * (3/8 + hn * (5/16 + hn * (35/128 + hn * 63/256)))))
hn = 1 - A / xn2
xn+1 = xn * (1 - hn * (1/2 + hn * (1/8 + hn * (1/16 + hn * (5/128 + hn * 7/256)))))
hn = 1 - A * xn3
xn+1 = xn + xn * hn * (1/3 + hn * (2/9 + hn * (14/81 + hn * (35/243 + hn * 91/729))))
hn = 1 - A * xn4
xn+1 = xn + xn * hn * (1/4 + hn * (5/32 + hn * (15/128 + hn * (195/2048 + hn * 663/8192))))
*Note: 今時のプロセッサではループの条件判断をしている間に中身の計算くらい出来るので、 ループを展開した方が良い。機械ε が10-16くらいならば3回分で十分である。#include <stdio.h> #include <math.h> #ifndef DBL_EPSILON #define DBL_EPSILON 2.2204460492503131E-16 #endif double _sqrtinv(a) double a; { double x, h, g; int e; /* enhalf the exponent for a half digit initial accuracy */ frexp(a, &e); x = ldexp(1.0, -e >> 1); /* 1/sqrt(a) 4th order convergence */ g = 1.0; while(fabs( h = 1.0 - a * x * x) < fabs(g)) { x += x * (h * (8.0 + h * (6.0 + 5.0 * h)) / 16.0); g = h; } return(x); } double sqrt(a) double(a); { if(a < 0){ fprintf(stderr, "sqrt: domain error\n"); return(0); } return(a * _sqrtinv(a)); }