近岡です。

正木さんの[ruby-math:00945] の証明を拝見しました。
まさしく、コロンブスの卵を見た気がします。

この証明を検討してみると、
[ruby-math:00943]のプログラムisqrt_iniの
  >     j = (bitlength - 1) >> 2
は
  >     j = (bitlength + 1) >> 2
でも、構わないんじゃないかと思うようになりました。

# 最初は、この証明を応用して、
#  >     j = (bitlength + 1) >> 2
# では上手くいかない事を説明できると考えていたのですが、
# 何時の間にか逆の結論が得られました。

>2 ** (4*n) <= a < 2**(4*n + 4)   , (n > 1)
>
>と仮定して
>
>b = a >> 2*n
>
>x = b.isqrt_ini << n
>
>と置く。
>
>b.isqrt_ini に関する帰納法の仮定で
>
>|b.isqrt_ini - sqrt(b)| < 1    (sqrt は本来の平方根:実数)
>
>a の平方根の範囲は
>
>sqrt(b)*2**n <= sqrt(a) < sqrt(b + 1) * 2**n 
>
>n > 1 ならば b >= 4 なので
>
>sqrt(b+1) - sqrt(b) < 1/4
>
>|b.isqrt_ini - sqrt(b + 1)| < 5/4
>
>となり、 x の誤差は
>
>| x - sqrt(a)| < 5/4 * 2 ** n  
>
>とおさえられる。

上の不等式は、もう少し制限を強めることができます。

今、仮に s=b.isqrt=[sqrt(b)](ただし[]はガウス記号)
とおくと、s ≦ sqrt(b) < s+1 から、
  s**2≦b<(s+1)**2
ここで、b も (s+1)**2 も整数なので
  s**2≦b<b+1≦(s+1)**2
従って、
  s≦sqrt(b)<sqrt(b+1)≦s+1
 
そうすると、前述の a の平方根の範囲の式
>sqrt(b)*2**n <= sqrt(a) < sqrt(b + 1) * 2**n 
から
  s*2**n ≦sqrt(a)<(s+1)*2**n …(*)
が得られます。

(ア) b.isqrt_ini=b.isqrt=sのとき、
 x=s*2**n で、(*)式より
  x ≦ sqrt(a) < x+2**n
(イ) b.isqrt_ini=b.isqrt+1=s+1のとき、
 x=(s+1)*2**nで、(*)式より
  x−2**n ≦ sqrt(a) < x

(ア)(イ)いずれの場合も、
  | x - sqrt(a)| ≦ 2 ** n  …(米)

# (米)式の等号成立は、aが平方数でしかも
#  sqrt(a)=s*2**n、x=(s+1)*2**n
# のとき。

>有理数での Newton 法の近似
>
>y = (x + a/x)/2      (y は有理数)
>
>に対して
>
>0 <=  y - sqrt(a) == (x - sqrt(a))**2/(2*x) < ((5/4)**2)/2 < 1
>
>sqrt(a) <= y < sqrt(a) + 1
>
>したがって
>
>a.isqrt <= y.floor < a.isqrt + 1
>
>Q.E.D.

ここで、| x - sqrt(a)| ≦ 2**n ですから
  (x - sqrt(a))**2/(2*x) ≦(2**(2*n))/(2*x)
従って、
  2**(2*n)≦2*x
ならば、
  y ≦ sqrt(a)+1
が成り立つはず。

元に返って考えて見ると、
 2**(4*n-2)≦a
でありさえすれば、
 2**(2*n-2)≦[a/(2**(2n))]=b
ですから、
 2**(n-1)≦[sqrt(b)]=s
となり、
 2**(2*n-1)≦s*2**n≦x
が成り立ちます。

もしかして、[ruby-math:00936]内のプログラム isqrtC の
効率が悪かったのは、
>     j = (bitlength + 1) >> 2
の影響じゃなくて、
>    x = ((self >> (j << 1)).isqrt_ini + 1) << j
の+1の影響だったんでしょうか?

# 確かに上の証明を見る限りにおいては、
# [ruby-math:00936]内のプログラム isqrtC において
#>    j = (bitlength + 1) >> 2
#>    x = ((self >> (j << 1)).isqrt_ini + 1) << j
# のとき、最悪の条件が続けば、ループを回す毎に
# 誤差が2→3→8→63→…のように
# 幾何級数より速く増大する可能性があります。

===(オマケ)============

Fixnum用のルーチンを別に用意する場合は関係ないことですが、
isqrt_ini中の
    return 0 if self < 1
    return 1 if self < 4
    return 2 if self < 9
    return 3 if self < 16
の部分を次のように置き換えても構わない。
(これらはいずれも[sqrt(self)]か[sqrt(self)]+1を返す)

(例1)
    return 0 if self < 1 
    return ((self >> 2) + 1 ) if self <= 19 
    # 1≦self≦7,9≦self≦11で正しい値を返す
(例2)
    return 0 if self < 1 
    return (self * 3  + 21) >> 4 if self <= 30 
    # 1≦self≦14で正しい値を返す
(例3)
    return (self * 5 + 12) >> 4 if self <= 7
    return (self * 9  + 112) >> 6 if self <= 44 
    # 前半は0≦self≦7で正しい値を返す
    # 後半は4≦self≦23で正しい値を返す
    # 「j = (bitlength - 1)>> 2」の代わりに
    # 「j = bitlength >> 2」「j = (bitlength + 1) >> 2」
    # とする場合は、後半を除く

0----+----1----+----2----+----3----+----4----+----5----+----6
近岡 宣吉  Chikaoka, Nobuyoshi
富山県立高岡西高等学校(数楽科)
 E-mail : chikaoka-nobuyoshi / tym.ed.jp