原です。

>正木です。

>もしかして何か勘違いしていて、間違いを指摘してもらえるかも知れ
>ないので一応概略を書いておきます。

うまく行きましたね。というわけで、もう、[ruby-math:00943] の isqrt 
は世界の教科書を書き換えるようなすばらしい式かも。

-----

ちょっと証明の補強します。プログラムでは

    x = (self >> (j << 1)).isqrt_ini
    ((x<<j) + (self>>j).div(x)) >> 1

となっていて、整数演算が入るので微妙なところがあります。つまり

  [(x*(2**j) + [[a/2**j]/x])/2]

なので、ここを厳密に評価します。

-----

今、正の実数 a に対して、k**4 <= a である自然数 k を取り、整数 
x > 1 に対し、問題の sqrt(a) 近似を作る式を N(x) とします。つま
り、ガウス記号を使うと

  N(x) = [(x*k + [[a/k]/x])/2]

です。このとき次が成り立ちます。

【命題】 -1 < x - sqrt(a)/k <= 1 ならば -1 < N(x) - sqrt(a) <= 1.

【証明】 r を整数とすると

  N(x) - r = [(x*k + [[a/k]/x])/2] - r
           = [[[x**2*k + a/k - 2*r*x]/x]/2]
           = [[[k*(x - sqrt(a)/k)**2 + 2*x*(sqrt(a) - r)]/x]/2]

ここで r = [sqrt(a)] とおく。上の [] の中の各項は >=0 なので

  N(x) - [sqrt(a)] >= 0. 

また、仮定より (x - sqrt(a)/k)**2 <= 1 であり、また
sqrt(a) - [sqrt(a)] < 1 なので、

  N(x) - [sqrt(a)] <= [[[k + 2*x]/x]/2].

再び仮定より k <= sqrt(a)/k < x+1 だから

  右辺 <= [[x+1 + 2*x]/x]/2] = [(3*x + 1)/x]/2] = 1.

以上より 0 <= N(x) - [sqrt(a)] <= 1 よって -1 < N(x) - sqrt(a) <= 1。

-----

今回のケースでは a.isqrt_ini のために j =(a.bitlength - 1)/4,
k = 2**j, x = (a/k**2).isqrt_ini として再帰させています。最後に
-1 < N - sqrt(a) <= 1 から N**2 > a ? N-1 : N で [sqrt(a)] が得
られるわけですね。因みに k は k**4 <= a を満たす計算しやすい値
ならなんでもいい。