近 岡 です。

4年前の話題で申し訳ありません。

[ruby-math:00909]にあるプログラム
> /* 2^30 <= x < 2^32 について root を計算する */
> unsigned long isqrt(unsigned long x)
> {
>     unsigned long a, b;
>     a = 143257 - 202860544/(1582+(x>>22));
>     a = (a + x/a)/2;
>     b = (a + x/a)/2;
>     return a > b ? b : a;
> }

についての問題です。平方根をニュートン法で求めることについて詳しい人はよ
くご存知でしょうが、このプログラム中の"a = (a + x/a)/2;"の部分は、ニュー
トン法でのループの中身そのものです。
[ruby-math:00909]で提起された問題は、このニュートン法の部分を1回に減ら
すことができるがその方法は? というものでした。

結論からすれば、次のように掛算1回余分に行うことで、ニュートン法のループ
部分を1回にできました。

/* 2^30 <= x < 2^32 について root を計算する */
unsigned long isqrt(unsigned long x)
{
    unsigned long a;
    a = 143257 - 202860544 /( 1582 + (x>>22));
    a = (a + x/a)/2;
    return ((a*a-1 >= x) ? a-1 : a);
}

定義域が2^2k≦x<2^(2k+2)で与えられている場合に、一次分数関数
   f(x) = A−C/(B+x/D)
で√xを近似するとき、係数をどんなにうまく選んでも、最大相対誤差は0.25%
を超えてしまうようである。
# 別の言い方をすると、有効数字が(2進法で)8ケタ程度。
ニュートン法の前に一次分数関数を使用した際の効果は、ニュートン法のループ
の2〜3回分に相当すると思われる。32ビット整数の範囲内では素晴らしい効
果があるが、RubyでのBIGNUMでは、目に見えるような効果は得られない。

/* ===============================================
    ilog2(x) := bitlength(x)-1 ≡ floor[log_2(x)]
    2進常用対数の整数部
 ※ 0≦x≦15の範囲で、一次分数関数を利用
=============================================== */
int ilog2(unsigned long x)
{
    int ans=0+4;  // 0<=x<=15のとき, ilog2(x)=4-[22/(x+4)]
  #if  ULONG_MAX > 0xFFFFFFFFUL
    while (x > 0xFFFFFFFFUL) { ans += 32; x >>= 32; }
  #endif
    if (x >= (1UL<<16)) { ans += 16; x >>= 16; }
    if (x >= (1UL<<8))  { ans +=  8; x >>=  8; }
    if (x >= (1UL<<4))  { ans +=  4; x >>=  4; }
    return ans - 22/(x+4);
}

/* ===============================================
    isqrt(n) := floor[√n] 平方根の整数部分
    √n ≒ a - [b / (c + [n>>d])] 1次分数関数による近似の利用
=============================================== */
static struct {
        unsigned long a,b,c; unsigned d;
} isqrt_coef[] = { /*  e=[log2(n)]  A=最大絶対誤差  R=最大相対誤差 */
 {     0x8UL,      0xc3UL,    0x18UL, 0x0,}, /* e=0 A=0 R=0.0 */
 {     0x8UL,      0xc3UL,    0x18UL, 0x0,}, /* e=1 A=0 R=0.0 */
 {     0x8UL,      0xc3UL,    0x18UL, 0x0,}, /* e=2 A=0 R=0.0 */
 {     0x8UL,      0xc3UL,    0x18UL, 0x0,}, /* e=3 A=0 R=0.0 */
 {     0x8UL,      0xc3UL,    0x18UL, 0x0,}, /* e=4 A=0 R=0.0 */
 {    0x10UL,     0x546UL,    0x57UL, 0x0,}, /* e=5 A=0 R=0.0 */
 {    0x1eUL,    0x224aUL,   0x13fUL, 0x0,}, /* e=6 A=1 R=0.100000 */
 {    0x28UL,    0x4f4cUL,   0x22dUL, 0x0,}, /* e=7 A=1 R=0.076923 */
 {    0x3bUL,    0xf8e2UL,   0x4aaUL, 0x0,}, /* e=8 A=1 R=0.062500 */
 {    0x52UL,   0x2929eUL,   0x8ebUL, 0x0,}, /* e=9 A=1 R=0.045455 */
 {    0x76UL,   0x794c6UL,  0x1253UL, 0x0,}, /* e=10 A=1 R=0.031250 */
 {    0xa6UL,  0x14f744UL,  0x241dUL, 0x0,}, /* e=11 A=1 R=0.022222 */
 {    0xedUL,  0x3cb25cUL,  0x495cUL, 0x0,}, /* e=12 A=1 R=0.015625 */
 {   0x14eUL,  0xa94628UL,  0x9159UL, 0x0,}, /* e=13 A=1 R=0.011111 */
 {   0x1ddUL, 0x1ebdba0UL, 0x12801UL, 0x0,}, /* e=14 A=1 R=0.007812 */
 {   0x29eUL, 0x55066d0UL, 0x24731UL, 0x0,}, /* e=15 A=1 R=0.005525 */
 {   0x3b7UL, 0xf2d518aUL, 0x495faUL, 0x0,}, /* e=16 A=1 R=0.003906 */
 {   0x53dUL,0x2a8339beUL, 0x91cc2UL, 0x0,}, /* e=17 A=1 R=0.002762 */
 {   0x76bUL,0x78a5a1e0UL,0x1243e8UL, 0x0,}, /* e=18 A=1 R=0.001953 */
 {   0xa31UL,0x9bf15560UL,0x112ea1UL, 0x1,}, /* e=19 A=2 R=0.002762 */
 {   0xe6aUL,0xdc97ef9cUL,0x112f65UL, 0x2,}, /* e=20 A=2 R=0.001953 */
 {  0x1462UL,0x9bf14ed9UL, 0x89750UL, 0x4,}, /* e=21 A=2 R=0.001381 */
 {  0x1cd4UL,0xdc97f9cdUL, 0x897b3UL, 0x5,}, /* e=22 A=2 R=0.000977 */
 {  0x28c5UL,0x9bfce628UL, 0x44bdfUL, 0x7,}, /* e=23 A=3 R=0.001036 */
 {  0x39a8UL,0xdc97e5f3UL, 0x44bd9UL, 0x8,}, /* e=24 A=3 R=0.000732 */
 {  0x5189UL,0x9bf727d4UL, 0x225e2UL, 0xa,}, /* e=25 A=4 R=0.000691 */
 {  0x734fUL,0xdc922bbcUL, 0x225e3UL, 0xb,}, /* e=26 A=5 R=0.000488 */
 {  0xa312UL,0x9bf727d4UL, 0x112f1UL, 0xd,}, /* e=27 A=6 R=0.000518 */
 {  0xe69fUL,0xdc951c26UL, 0x112f4UL, 0xe,}, /* e=28 A=8 R=0.000427 */
 { 0x14627UL,0x9bfb6f51UL,  0x897bUL,0x10,}, /* e=29 A=11 R=0.000432 */
 { 0x1cd3eUL,0xdc951c26UL,  0x897aUL,0x11,}, /* e=30 A=15 R=0.000397 */
 { 0x28c47UL,0x9bf65f0bUL,  0x44bcUL,0x13,}, /* e=31 A=22 R=0.000436 */
#if ULONG_MAX > 0xFFFFFFFFUL
 { 0x39a7cUL,0xdc951c26UL,  0x44bdUL,0x14,}, /* e=32 A=32 R=0.000436 */
 { 0x5188eUL,0x9bf65f0bUL,  0x225eUL,0x16,}, /* e=33 A=54 R=0.000531 */
 { 0x734eaUL,0xdc902f99UL,  0x225eUL,0x17,}, /* e=34 A=77 R=0.000536 */
 { 0xa311bUL,0x9bf64811UL,  0x112fUL,0x19,}, /* e=35 A=150 R=0.000743 */
 { 0xe69d5UL,0xdc9047ffUL,  0x112fUL,0x1a,}, /* e=36 A=213 R=0.000742 */
 {0x1462d2UL,0x9c041587UL,   0x898UL,0x1c,}, /* e=37 A=467 R=0.001160 */
 {0x1cd487UL,0xdca3d16eUL,   0x898UL,0x1d,}, /* e=38 A=661 R=0.001162 */
#endif
};
#define	_A(e)	isqrt_coef[e].a
#define	_B(e)	isqrt_coef[e].b
#define	_C(e)	isqrt_coef[e].c
#define	_D(e)	isqrt_coef[e].d

unsigned long isqrt(unsigned long n)
{
    int e,e2; unsigned long s,t;

    if (n<=35) return 8-195/(24+n);

    e = ilog2(n);  /* 2^e≦ n <2^(e+1) */

    if (e < 19) {    // 36≦n<2^19
        t = _A(e) - _B(e) / (_C(e) + n);
        return ( (t*t > n) ? t-1 : t);
    }
    else {
  #if  (ULONG_MAX>>39) < 1    /* このとき n<2^39 */
        t = _A(e) - _B(e) / (_C(e) + (n>>_D(e)));
        t = (n/t + t)>>1;
        return ( (t*t-1 >= n) ? t-1 : t );
  #else
        if (e < 39) {    // 2^19≦n<2^39
            t = _A(e) - _B(e) / (_C(e) + (n>>_D(e)));
            t = (n/t + t)>>1;
            return ( (t*t-1 >= n) ? t-1 : t );
        }
        else {    // 2^39≦n
            e2 = (e - 30)>>1;
            t = 143211UL-3245822441UL/(25315UL+(n>>(e2+e2+18)));
            t = ( (n>>e2)/t + (t<<e2) )>>1;
            do { s = t; t = ((n/t)+t)>>1; } while (t<s);
            return s;
        }
  #endif
    }
}
#undef 	_A()
#undef 	_B()
#undef 	_C()
#undef 	_D()