原です。

大分間があいてしまいました。

In [ruby-math:00965] 

>近岡です。

>再帰呼び出しを繰り返さない isqrt です。

すばらしいです。でも、確かに難しいですねえ。はじめぜんっぜんわ
かりませんでした。しかし、検証していくうちに証明ができたので、
載せておきます。


次のコードは近岡さんの [ruby-math:00965] のものを若干短くした
だけのもので、ループ中の x の値は完全に同じです。後、l が 1 
と区別しにくかったんで e に変えました。(^^;


【コード C】 以下は [sqrt(self)] を求めるアルゴリズムである。

class Integer
  def isqrtF1
    raise RangeError if self < 0
    return [0,1,1,1,2,2,2,2][self] if self < 8

    e = (bitlength - 2) & ~1
    x = [0,1,1,1,2,2,2,2][self >> e]

    s = 0
    (e.bitlength - 1).downto 1 do |i|
      t = s
      s = e >> i
      x = ((x << (t + (s & 1))) + (self >> (e - s - t)).div(x)) >> 1
    end

    (x * x > self) ? (x - 1) : x
  end
end

で、これが正しい isqrt を返す事を示します。まず、ニュートン近似 
1 回分を N(x, a) と表しておきます。(以下、a を正の実数、log は 
2 を底とする対数、/ は実数の割り算とします。)


【定義1】 N(x, a) = (x + a/x)/2


次に log(log(a)) 回の N(x, a) による近似数列の一般的な定理をいっ
ておきます。


【定理 2】 n を自然数、u(0), u(1), ..., u(n) を単調増大な非負整
  数列で次を満たすとする:

    (u1)  u(0) = 0.
    (u2)  u(i) - u(i - 1) は 偶数  (0 < i <= n).
    (u3)  u(i)*2 - u(i - 1) <= log(a) + 2  (0 < i <= n).

  このとき、非負整数列 x(0), x(1), ... x(n) を x(n) から次のように定める:

    (xi) 0 < i <= n となる i について

         k(i) = 2**((u(i) - u(i - 1))/2).

      とおき、

         x(i - 1) = [N(k(i) * x(i), a * 2**(-u(i - 1)))].

       と帰納的に定義する。
       

  このときもし

    (xn) |x(n) - sqrt(a * 2**(-u(n)))| < 1.

  が成り立っているとすると、全ての i  (0 <= i < n) についても

    |x(i) - sqrt(a * 2**(-u(i)))| < 1.

  が成り立つ。特に

    |x(0) - sqrt(a)| < 1.

  である。


後でこの定理がコード C に適用できることを証明します。まずこの定
理の証明のために、前に [ruby-math:00957] で証明した正木さんのア
ルゴリズムの核心部分を定理 2 に適用しやすいように書き換えて証明
しておきます。内容は同じです。


【定理 3】x, k を自然数で k = 1 または偶数とする。このとき
  任意の実数 b >= k/2 に対し、次が成り立つ。

    |x - b| < 1 ならば |[N(k*x, (k*b)**2)] - k*b| < 1.

(証明)
  N(k*x, (k*b)**2) - k*b = k/(2*x) * (x - b)**2
                         < k/(2*x)

  一方、 k/2 <= b < x + 1 よって k/2 < x + 1. 仮定より k/2 = 1/2 
  または整数だから k/2 <= x. よって k/(2*x) <= 1. よって

    0 <= N(k*x, (k*b)**2) - k*b < 1.

  よって

    |[N(k*x, (k*b)**2)] - k*b| < 1.
(証明終)


(定理 2 の証明)全ての i  (0 < i <= e) について、

    |x(i) - sqrt(a * 2**(-u(i)))| < 1
                        ==> |x(i - 1) - sqrt(a * 2**(-u(i - 1)))| < 1

  がいえればよい。ここで

    x = x(i), k = k(i), b = sqrt(a * 2**(-u(i)))

  とおくと

     k*b = 2**((u(i) - u(i - 1))/2) * sqrt(a * 2**(-u(i)))
         = sqrt(a * 2**(-u(i - 1)))

  なので、(xi) より

    x(i - 1) = [N(k*x, (k*b)**2)]

  である。定理 3 の条件の成立を調べると、k は条件 (u2) より整数で
  あり、1 または偶数であるのはあきらか。また
  
    log(b) = (log(a) - u(i))/2
    log(k/2) = (u(i) - u(i - 1))/2 - 1

  なので、条件 (u3) より log(b) >= log(k/2)、よって  b >= k/2.
  よって定理 3 が適用できて証明された。
(定理 2 の証明終)


【コード C の正当性の証明】

  コード C に、定理 2 を適用することを考える。以下 a = self として
  表現する。


  《初期値について》

    a >= 8 としてよい。最初に

      e = [([log(a)] - 1)/2]*2  (= (a.bitlength - 2) & ~1)
      n = [log(e)]              (= e.bitlenght - 1)

    としているので、

      log(a) - e =  log(a) - (([log(a)] - 1) & ~1)
                 <= log(a) - ([log(a)] - 2)  <  3.

    よって (a >> e) < 2**3 = 8 だから
 
       x = [0,1,1,1,2,2,2,2][a >> e] 

    は整数値を返し、配列の選び方から

       x = [sqrt(a >> e)]
         = [sqrt(a * 2**(-e))]

    である(一般に非負の実数 b について [sqrt([b])] = [sqrt(b)] )。


  《u(i) の設定と条件 (u1), (u2), (u3) について》

    数列 u(0), u(1), ..., u(n) を

      (Cu)  u(i) = e - ((e >> i) & ~1)

    で定義する。このとき u(i) が条件 (u1), (u2) を満たすのは明らか。
    (u3) は、s = e >> (i - 1) とおくと、
  
      u(i)*2 - u(i - 1) = (e - ((s >> 1) & ~1))*2 - (e - (s & ~1)),
                        = e + (s & ~1) - ((s >> 1) & ~1) * 2,
                        = e + ((s >> 1) & 1) * 2,
                        <= e + 2,
                        <= log(a) + 2.

    なので、成り立っている。


  (注意)上の変形では、整数 z に対する次の恒等式を使っている。
    (i)  z = (z & 1) + (z & ~1).
    (ii) z & ~1 = (z >> 1) + (z >> 1).


  《条件 (xn) について》

      u(n) = e - (e >> n) & ~1 = e

    なので、x(n) = x = [sqrt(a * 2**(-e))] とおけば、

      |x(n) - sqrt(a * 2**(-u(n)))| = |x - sqrt(a * 2**(-e))|,
                                    < 1.

  《条件 (xi) について》

    while ループにおける x の更新:

      x = ((x << (t + (s & 1))) + (a >> (e - s - t)).div(x)) >> 1

    を調べる。

      s = e >> i
      t = e >> (i + 1) = s >> 1
      (右辺の x) = x(i) 

    とおくと

      u(i) = e - ((e >> i) & ~1) = e - 2*t.
      u(i - 1) = e - ((e >> (i - 1)) & ~1) = e - 2*s.

    から(j(i) = (u(i) - u(i - 1))/2 とおいて)、

      t + (s & 1) = s - t = j(i).
      e - s - t = (e - 2*s) + (s - t),
                = u(i - 1) + j(i).

    がいえるので、これらを x の右辺に代入すると

      x = ((x(i) << j(i)) + (a >> (u(i - 1) + j(i))).div(x(i))) >> 1

    であり、 [ruby-math:00957] の命題の変形を用いて更に計算すると、
     2**j(i) = k(i) だから、

      x = [(x(i) * k(i) + a * 2**(-u(i - 1)) / (k(i) * x(i)))/2],
	= [N(x(i) * k(i), a * 2**(-u(i - 1)))],
        = x(i - 1).

    である。すなわち、この部分の x の更新は、帰納的に数列 x(i) から 
    x(i - 1) を定義している。


  《結論》

     以上より定理 2 の全ての条件が満たされることがいえた。よって、
     x = x(0) に対して

       |x - sqrt(a)| < 1.

     がいえる。従って

       (x ** 2 > a) ? : x - 1 : x

     は [sqrt(a)] に等しくなる。

(証明終)