In [ruby-math:00936]

>正木です。

とにかく「正木法」のニュートン法の初期値をニュートン法で求める
というアイデアにはシビレました。あとは細かい話になりそうなんで
すが、、、。

>In [ruby-math:00935]
>|近岡です。
>
>| ですから、
>|おおもとのrb_int_sqrt()内では、do 〜 while ループは欠かせませんが、
>|再帰呼び出しされている側のrb_int_sqrt()内では、do 〜 while ループを
>|削除して、単純に、
>|  return rb_rshift(rb_plus(z, rb_div(x,z)), 1);
>|とすることが可能だと思います。

近岡さんの話ももっともで、ニュートン法は初期値が良ければ、特に最
後の一回の割り算は多くの場合無駄なので、省略したら速くなるかも、
という事ですよね。私の実験(後述)では若干速くなりました。

しかし、とにかく再帰のオーダーは O(log(log(N))) なので、少なくとも
スタック枯渇は心配しなくていいみたいです。20000! (十進で7万桁!)
でも深さは 80 回ぐらいでした。(不思議な事にメインループは回らず、
初期値だけで既に答えが出ていました。)


>  def isqrt
>    raise "argument is negative" if self < 0
>    return 0 if self == 0
>    return 1 if self < 4
>    j = (bitlength + 1) >> 2
>    x = ((self >> (j << 1)).isqrt + 1) << j
>    while x > y = div(x)
>      x = (x + y) >> 1
>    end
>    x
>  end

これはオリジナルの [ruby-math:00933] を少し変形していますよね。
多分オリジナルのままの方がちょっと速いと思います。つまり、

  def isqrt2
    raise "argument is negative" if self < 0
    return 0 if self == 0
    return 1 if self < 4
    return 2 if self < 9
    return 3 if self < 16

    j = (bitlength - 1) >> 2
    x = ((self >> (j << 1)).isqrt2 + 1) << j
    while x > y = div(x)
      x = (x + y) >> 1
    end
    x
  end

でいいのではないでしょうか。

ところで正木さんは、よく見るニュートン法のループ

    begin
      z = x
      x = (x + div(x)) >> 1
    end while (x < z)
    return z

も改良してますね。「m > n <==> m > (m+n).div(2)」なわけか、
なるほど。


>  def isqrt_ini
>    return 0 if self == 0
>    return 1 if self < 4
>    j = (bitlength + 1) >> 2
>    x = ((self >> (j << 1)).isqrt_ini + 1) << j
>    (x + div(x)) >> 1
>  end
>
>  def isqrtC
>    raise "argument is negative" if self < 0
>    return 0 if self == 0
>    return 1 if self < 4
>    x = isqrt_ini + 1
>    while x > y = div(x)
>      x = (x + y) >> 1
>    end
>    x
>  end
>end

こちらもこれでいいと思います。

  def isqrt1_ini
    return 0 if self < 1
    return 1 if self < 4
    return 2 if self < 9
    return 3 if self < 16

    j = (bitlength - 1) >> 2
    x = ((self >> (j << 1)).isqrt1_ini + 1) << j
    (x + div(x)) >> 1
  end

  def isqrtC1
    raise "argument is negative" if self < 0
    return 0 if self == 0
    return 1 if self < 4

    x = isqrt1_ini
    while x > y = div(x)
      x = (x + y) >> 1
    end
    x
  end


で、さらにたった今投げられた、近岡さんの [ruby-math:00939] を見たら

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

は、

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

としてますね。ううむ、なるほど!気がつかなかった。こちらの方が割
り算が軽くなりますね。これを正木さんのオリジナルに取り込むと、、、


  def isqrt2_ini
    return 0 if self < 1
    return 1 if self < 4
    return 2 if self < 9
    return 3 if self < 16

    j = (bitlength - 1) >> 2
    x = (self >> (j << 1)).isqrt2_ini + 1
    ((x<<j) + (self>>j).div(x)) >> 1
  end

  def isqrtC2
    raise "argument is negative" if self < 0
    return 0 if self == 0
    return 1 if self < 4

    x = isqrt2_ini
    while x > y = div(x)
      x = (x + y) >> 1
    end
    x
  end


というわけで、ベンチマークです。

isqrt   ... [ruby-math:00936] 
isqrt2  ... 上述
isqrtC1 ... 上述
isqrtC2 ... 上述
isqrtD  ... [ruby-math:00939] 

                          user     system      total        real
isqrt(10000!)         0.260000   0.000000   0.260000 (  0.253867)
isqrt2(10000!)        0.250000   0.000000   0.250000 (  0.251312)
isqrtC1(10000!)       0.220000   0.000000   0.220000 (  0.216147)
isqrtC2(10000!)       0.150000   0.000000   0.150000 (  0.156600)
isqrtD(10000!)        0.160000   0.000000   0.160000 (  0.155240)

                          user     system      total        real
isqrt(11000!)         0.250000   0.000000   0.250000 (  0.251974)
isqrt2(11000!)        0.260000   0.000000   0.260000 (  0.250360)
isqrtC1(11000!)       0.220000   0.000000   0.220000 (  0.219876)
isqrtC2(11000!)       0.150000   0.000000   0.150000 (  0.157385)
isqrtD(11000!)        0.150000   0.010000   0.160000 (  0.156833)

                          user     system      total        real
isqrt(12000!)         0.300000   0.000000   0.300000 (  0.305972)
isqrt2(12000!)        0.310000   0.000000   0.310000 (  0.307817)
isqrtC1(12000!)       0.370000   0.000000   0.370000 (  0.371402)
isqrtC2(12000!)       0.300000   0.000000   0.300000 (  0.300848)
isqrtD(12000!)        0.190000   0.000000   0.190000 (  0.190973)

                          user     system      total        real
isqrt(13000!)         0.610000   0.000000   0.610000 (  0.612236)
isqrt2(13000!)        0.440000   0.000000   0.440000 (  0.444748)
isqrtC1(13000!)       0.390000   0.000000   0.390000 (  0.385644)
isqrtC2(13000!)       0.280000   0.000000   0.280000 (  0.281237)
isqrtD(13000!)        0.280000   0.000000   0.280000 (  0.281503)

                          user     system      total        real
isqrt(14000!)         0.420000   0.000000   0.420000 (  0.414214)
isqrt2(14000!)        0.410000   0.000000   0.410000 (  0.413504)
isqrtC1(14000!)       0.360000   0.000000   0.360000 (  0.363345)
isqrtC2(14000!)       0.260000   0.000000   0.260000 (  0.258511)
isqrtD(14000!)        0.260000   0.000000   0.260000 (  0.260459)

                          user     system      total        real
isqrt(15000!)         0.450000   0.010000   0.460000 (  0.463047)
isqrt2(15000!)        0.460000   0.000000   0.460000 (  0.461631)
isqrtC1(15000!)       0.400000   0.000000   0.400000 (  0.404377)
isqrtC2(15000!)       0.290000   0.010000   0.300000 (  0.290540)
isqrtD(15000!)        0.280000   0.000000   0.280000 (  0.287397)

                          user     system      total        real
isqrt(16000!)         0.540000   0.000000   0.540000 (  0.545825)
isqrt2(16000!)        0.550000   0.000000   0.550000 (  0.543977)
isqrtC1(16000!)       0.470000   0.000000   0.470000 (  0.477087)
isqrtC2(16000!)       0.350000   0.000000   0.350000 (  0.344756)
isqrtD(16000!)        0.340000   0.000000   0.340000 (  0.341518)

                          user     system      total        real
isqrt(17000!)         0.640000   0.000000   0.640000 (  0.639905)
isqrt2(17000!)        0.640000   0.000000   0.640000 (  0.644164)
isqrtC1(17000!)       0.560000   0.000000   0.560000 (  0.560813)
isqrtC2(17000!)       0.400000   0.000000   0.400000 (  0.402403)
isqrtD(17000!)        0.410000   0.000000   0.410000 (  0.403905)

                          user     system      total        real
isqrt(18000!)         1.110000   0.000000   1.110000 (  1.104188)
isqrt2(18000!)        0.800000   0.000000   0.800000 (  0.804987)
isqrtC1(18000!)       0.690000   0.000000   0.690000 (  0.688353)
isqrtC2(18000!)       0.500000   0.000000   0.500000 (  0.499481)
isqrtD(18000!)        0.500000   0.000000   0.500000 (  0.498362)

                          user     system      total        real
isqrt(19000!)         0.820000   0.000000   0.820000 (  0.817477)
isqrt2(19000!)        0.820000   0.000000   0.820000 (  0.818653)
isqrtC1(19000!)       0.710000   0.000000   0.710000 (  0.715676)
isqrtC2(19000!)       0.510000   0.000000   0.510000 (  0.510921)
isqrtD(19000!)        0.510000   0.000000   0.510000 (  0.511898)

                          user     system      total        real
isqrt(20000!)         0.990000   0.000000   0.990000 (  0.994729)
isqrt2(20000!)        0.990000   0.000000   0.990000 (  0.990165)
isqrtC1(20000!)       0.870000   0.000000   0.870000 (  0.869424)
isqrtC2(20000!)       0.630000   0.000000   0.630000 (  0.638718)
isqrtD(20000!)        0.630000   0.000000   0.630000 (  0.629519)


結論としては isqrtC2 と isqrtD が速いという事でしょう。

因みに拡張ライブラリとして C で書いても*全く*速くなりません
でした。もともと rb_funcall を起動する回数がとても小さいから
でしょう。


ところで、近岡さんの isqrtD ([ruby-math:00939]) で

>  def isqrtD_ini
>    return 0 if self == 0
>    if self < 16 
>        j = bitlength >> 1
>        x = ((1<<j)+(self>>j))>>1
>        while x > y=div(x)
>          x = (x + y)>>1
>        end
>        x + 1  # 大き目の値を返す
>    else
>      j = (bitlength - 1) >> 2
>      x = ((self >> (j << 1)).isqrtD_ini)
>      ((x<<j) + (self>>j).div(x)) >> 1
>    end
>  end


となっていますが、この「大きめ」を得るにはこの方法で十分でしょ
うか?一回でもニュートン法のループ回すと「大きめ」という性質は
崩れてしまう可能性があるのでは?

また、

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

という工夫も値が小さめになるので、isqrtC2 でさえちょっと不安です。


正木さんのオリジナル isqrt2 と isqrtC1 の、再帰的に初期値を求め
て +1 してからビットシフトする方法は、正しい事が証明できました。


メモ:ポイントは次の補題です。([] はガウス記号)

【補題】非負の実数 a に対して [sqrt([a])] + 1 > sqrt(a)。
(証明)n = [sqrt(a)] とおくと a >= n**2 なので、
       [sqrt([a])] + 1 >= [sqrt([n**2])] + 1 >= n + 1 > sqrt(a).

また忘れそうなので、ほとんど自分へのメモでした。(^^;