Newton 法について

倍精度浮動小数計算では精度が不足であっという間に収束するので, 任意精度浮動小数計算 (bigfloat) を使う. プログラムは次のように なる.

/* A は小数点が付かない値を指定 */
/* Digit では必要な桁数を指定する. */

def sqrt(A,Digit) {
    setprec(Digit);
    Eps = eval(10^(-Digit)*exp(0));
    X = eval(A*exp(0));
    while ( 1 ) {
        Y = (X+A/X)/2;
        if ( abs(X-Y) < Eps )
            return Y;
        else
            X = Y;
    }
}

def abs(X) { return X >= 0 ? X : -X; }
end$

ここで, eval(X) は, X が超越関数の評価を含む場合には bigfloat の値を返すので, exp(0) ( = 1) という値を かけることで強制的に bigfloat にしている. このプログラムで sqrt(2,1000) とすると $\sqrt{2}$ が 1000 桁程度求まる. eval(2^(1/2),1000) と比較してみよ.