√6 - 計算用紙

某所に $\sqrt{6}$ を計算しろという問題があったので、前に書いた円周率のプログラムをいじって書いてみた。
$\Large \sqrt{6} = \frac{5}{2} \sqrt{1-\frac{1}{25}}\hspace{50}(1)$
として √ のところをテーラー展開して
$\Large \sqrt{6} = \frac{5}{2} \left(1-\frac{1}{50}-\frac{1!!}{2!\cdot 50^2}-\frac{3!!}{3!\cdot 50^3}-\frac{5!!}{4!\cdot 50^4}-\frac{7!!}{5!\cdot 50^5}-\cdots\right)$
これで10進2桁ずつ求めていくというやり方。
速くはないけど、プログラムは簡単に書ける。

平方根をこの方法で計算するとき、 $|b|\ll 1$ として、(1) のように
$\Large a \sqrt{1+b}$
という形のものを探すんだけど、このとき、a, b の分母が10の冪乗の約数、b の分子が 1、b＜0 となっていると都合がいい(b＞0 だと級数が交代級数になって少し面倒)。
面白いことに、 $\sqrt{6}$ はこれらの条件を全て満たした形にできる。

下がプログラム。言語は C。
M の 1.4( $\sim \log_{10}{25}$ )倍程度の桁数まで計算して止まるから、欲しい桁数に応じて M を調節すればいい。
120万桁くらいまで正しく計算できるはず。
The first 1 million digits of the square root of 6. - NASA と比べてみると100万桁までの結果は正しく出力されているようだ。

#include <stdio.h>
#define M 100000
#define N 10
unsigned a[M], b[N], j, k, m, s, x;
main(i)
{
    puts("2.");
    for (i = 1 - N, x = 5, s = 55; m < M - 1; i++, m = k, s = 0) {
        for (k = 2; k < m || x; k++) {
            x = 100 * a[k] + x * (2 * k - 3);
            a[k] = x % (50 * k);
            s += x /= 50 * k;
        }
        for (j = 0; j < N; j++) {
            s += 100 * b[j];
            b[j] = s % 100;
            s /= 100;
        }
        if (i > 0) {
            printf("%02d", 99 - s);
            if (i % 500 == 0) printf("\n[%d]", 2 * i);
            if (i % 25 == 0) puts("");
            else if (i % 5 == 0) printf(" ");
        }
    }
}