原です。

Rational のスレッドを切ろうと思ったのですが、始まったものはしかたが
ない。テクニカルな話はこのまま行きましょうか。(^^;

>近岡です。

>[ruby-math:00907]のプログラムではオーダーが
>        O(log(i)^2×log(i))=O(log(i)^3)、
>ニュートン法ではオーダーが
>        O(log(i)×log(i))+O(log(i)^2×log(log(i)))=O(log(i)^2×log(log(i)))
>になります。

近岡さんの見積もりの正当性は私にはちょっと分かりません。例えば
かけ算のコストについても、各ステップで桁数が大きく変動している
ので、評価はなかなかむずかしいと思います。確かに乗法、除法、な
ど重さが分からないと、どうしようもないですよね。

Bignum のビットシフトについても無視できないという噂もあります。


># 整数の2乗の部分をビットずらしと加法のみで実現すると、
># オーダーがO(log(i))で実現できるのでは???
>#
># 例えば、rb_int_sqrt(x)の中の((n<<1)+1)の2乗を求めている
>#>    else {
>#>        n = rb_big_lshift(y,INT2FIX(1));
>#>      m = rb_big_or(n,INT2FIX(1));
>#>      m2 = rb_big_mul(m,m);
>#>    }
># の部分は、すでに前回のループで n*n を求めているのですから、
># これを変数 n2 として保存してあれば、
>#   m2=((n2+n)<<2) + 1
># と、加法2回とビットずらし1回ですむのでは???

つまり、

  def kaiheiho
    return 0 if self == 0
    s = ilog2 / 2
    t = a = 0
    (s + 1).times do |i|
      x = self >> ((s - i) << 1)
      a <<= 1
      t <<= 2
      u = t + (a << 1) + 1
      if u <= x
        a |= 1
        t = u
      end
    end
    a
  end

こんな感じですよね。これをほぼそのまま C で書くと

static VALUE
kaiheiho(x)
    VALUE x;
{
  long i, s = int_ilog2(x) / 2;
  VALUE t = INT2FIX(0), a = INT2FIX(0), y, u;
  for (i = 0; i <= s; i++) {
    y = rb_rshift(x, LONG2FIX((s - i) << 1));
    a = rb_lshift(a, INT2FIX(1));
    t = rb_lshift(t, INT2FIX(2));
    u = rb_plus(rb_plus(t, rb_lshift(a, INT2FIX(1))), INT2FIX(1));
    if (RTEST(rb_le(u, y))) { // u <= y
      a = rb_funcall(a,  '|', 1, INT2FIX(1));
      t = u;
    }
  }
  return a;
}

みたいになって、実際計測すると格段にはやくなり、[ruby-math:00912]
の正木さんのニュートン法と 10000! までほぼ同じオーダーになりまし
た。これは乗算の重複を省いただけでなく、再帰呼び出しをループにし
た効果が大きいのだと思います。それでもニュートン法に比べれば大分
遅いのですが。

そして、どうも走り始め開平法が強く、ある程度近似値が得られるとニュ
ートン法が強いという印象があったので、両方をミックスしたらどうなる
か実験してみました。すると全般に速く、10000! では純ニュートン法の
3倍程度速いプログラムが得られました。

えーと、参考までに一応載せておきます。

static VALUE
isqrt_ini(x, magic)
     VALUE x;
     long magic;
{
  long i, s, k;
  VALUE t, a, y, u;
  s = int_ilog2(x) / 2;
  if (magic > 0) {
    k = s / magic;
  }
  else {
    k = s + 1;
  }
  t = INT2FIX(0);
  a = INT2FIX(0);
  k = k > s + 1 ? s + 1 : k;
  for (i = 0; i < k; i++) {
    y = rb_rshift(x, LONG2FIX((s - i) << 1));
    a = rb_lshift(a, INT2FIX(1));
    t = rb_lshift(t, INT2FIX(2));
    u = rb_plus(rb_plus(t, rb_lshift(a, INT2FIX(1))), INT2FIX(1));
    if (RTEST(rb_le(u, y))) {
      a = rb_funcall(a,  '|', 1, INT2FIX(1));
      t = u;
    }
  }
  a = rb_lshift(rb_plus(a, INT2FIX(1)), INT2FIX(s + 1 - k));
  return a;
}

static VALUE
isqrt_haraK(x)
    VALUE x;
{
    VALUE y, z;
    if (FIXNUM_P(x)) return INT2FIX(int_isqrt_hara(FIX2INT(x))); 
    if (!RBIGNUM(x)->sign)  rb_raise(rb_eTypeError,"argument is negative");
    z = isqrt_ini(x, 16);
    do {
        y = z;
        z = rb_rshift(rb_plus(y, rb_div(x,y)), Unity);
    } while(rb_greater(y, z));
    return y;
}