近 岡 です。
４年前の話題で申し訳ありません。
[ruby-math:00909]にあるプログラム
> /* ２^30 <= x < ２^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]で提起された問題は、このニュートン法の部分を１回に減ら
すことができるがその方法は？ というものでした。
結論からすれば、次のように掛算１回余分に行うことで、ニュートン法のループ
部分を１回にできました。
/* ２^30 <= x < ２^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);
}
定義域が２^2k≦ｘ＜２^(2k+2)で与えられている場合に、一次分数関数
ｆ(ｘ) ＝ Ａ−Ｃ／(Ｂ＋ｘ／Ｄ)
で√ｘを近似するとき、係数をどんなにうまく選んでも、最大相対誤差は0.25％
を超えてしまうようである。
＃ 別の言い方をすると、有効数字が（２進法で）８ケタ程度。
ニュートン法の前に一次分数関数を使用した際の効果は、ニュートン法のループ
の２〜３回分に相当すると思われる。３２ビット整数の範囲内では素晴らしい効
果があるが、RubyでのBIGNUMでは、目に見えるような効果は得られない。
/* ===============================================
ilog2(x) := bitlength(x)-1 ≡ floor[log_2(x)]
２進常用対数の整数部
※ ０≦x≦15の範囲で、一次分数関数を利用
=============================================== */
int ilog2(unsigned long x)
{
int ans=0+4; // 0<=x<=15のとき, ilog2(x)=４-[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])] １次分数関数による近似の利用
=============================================== */
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); /* ２^e≦ n ＜２^(e+1) */
if (e < 19) { // 36≦ｎ＜２^19
t = _A(e) - _B(e) / (_C(e) + n);
return ( (t*t > n) ? t-1 : t);
}
else {
#if (ULONG_MAX>>39) < 1 /* このとき ｎ＜２^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) { // ２^19≦ｎ＜２^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 { // ２^39≦ｎ
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()
