整数の平方根アルゴリズム（sqrt64）

2025-12-31
Takayuki HOSODA

概要

これは、18〜24-bit の信号処理システム用に設計した、最小誤差丸めの 64-bit 平方根（square root）アルゴリズムと C での実装例です。

64-bit の整数の平方根は 32-bit になるので、係数の乗算や右シフトでのオーバーフローやアンダーフローを防ぐため、 MSB の 8-bit 下位付近で演算を行っています。
0x8000 未満の数は 16-bit 上位にスケーリングしてから計算しています。
初期値を １次の Min-Max 的近似で求め、 次いでニュートン法を３回施し、1 LSB 以内の誤差が得られた後に、LSB を調整し、誤差を ±1/2 LSB に収めています。

アルゴリズム

x0 = (22 a + 11) / 32 | a ∈ (0.25, 1) … (1)
xi+1 ← (xi + a / xi ) / 2 … (2)

本アルゴリズムは 式 (1) と (2) を中核とする。

初期近似

区間 x ∈ (0.25, 1) における平方根の初期近似として、 定数の除算が 2 の冪になるように調整した 1 次の近似式 (1) を用いる。
この初期近似で、± 0.31 % 以内の精度が得られる。

誤差低減

誤差の低減には、(2) 式のニュートン法を用いている。
重解を持たない場合、ニュートン法では解の近傍では 2 次の収束を示す。
繰り返しごとの理論精度は次のようになる。

1st: 4.8×10-4
2nd: 1.2×10-7
3rd: 6.5×10-15

1 / 232 ≈ 2.33×10-10 ≫ 6.5×10-15 であるので、 ニュートン法の 3 回の繰り返しで十分と言える。

ところで、ニュートン法の (2) 式は f (x) = xi2 - a = 0 として、 次のようにして導かれる：

xi+1 = xi - f (xi) / f '(xi)
   = xi - (xi2 - a) / (2 xi)
   = (xi + a / xi) / 2
と言えるのであるが、次のようにも捉えることができる。
誤差項を h = 1 - a / xi2 として、 √(1 - h) の、例えば 4 次までの級数展開は次のように表される。
h = 1 - a / xi2 … (3)
series(√(1 - h), h, 4) = 1 - h / 2 - h2 / 8 - h3 / 16 + O(h4) … (4)
この級数展開を 1 次の項まで用いて、xi を補正すると、
xi+1 = xi (1 - h / 2)
   = (xi + a / xi) / 2
となって、(2) 式と同じになる。

ここで、(4) 式の 2 次の項に着目すると、-h / 8 は負であるから、 h が十分小さい場合、これ以降の補正項の和は負であると言える。
従って、(2) 式は、その分だけ常に補正不足であり、近似誤差は常に正となる。

ところで、(4) 式を 2 次まで用いると、Halley 法と等価であるし、3 次まで用いると Housholder 法と等価である。

±1/2 LSB 誤差への調整

与えられた整数を x 、計算結果の整数を y とすると

y2 < x < (y + 1)2 … (5)
となる y が存在し、y に対して誤差が最大の 1/2 LSB となるときの x は
x == (y + 1/2)2 … (6)
であるが、 x は整数なので x == y^2 + y のときである。 従って、 誤差を 1/2 LSB に納めるためには、 y2 + y < x の場合に y + 1 を答えとすれば良い。

誤差分布

誤差のヒストグラムは、例えば 38-bit システムにおける 0〜238 - 1 ( = 274877906943 (0x3fffffffff)) の全数検査では、次のようになる。 

error(-1/4..-1/2 LSB) =  68719476736
error(-1/4..+1/4 LSB) = 137438953472
error(+1/4..+1/2 LSB) =  68719476736

C 実装

sqrt64 -- 与えられた 64-bit 整数の平方根に対して誤差が最小となる整数を返す。 負数の数に対しては、例外を発生させずそのまま返す。 

/* sqrt64 - calculate squre root of given integer value
 * Rev.1.4 (2025-12-31) (c) 2003-2025 Takayuki HOSODA
 * SPDX-License-Identifier: BSD-3-Clause
 */
#include <inttypes.h>
int64_t sqrt64(int64_t x);
int64_t sqrt64(int64_t x) {
    int64_t a;                          // approximation register
    int64_t t;                          // scaled '1'
    int     sr;                         // scaling (for relative domain of (0.25, 1))
    int     sx;                         // scaling (for small values)
    int     n_bit = sizeof(x) * 8;      // number of digits of given x (e.g. 64 for int64_t)
    int     n_room = 8;                 // head room
    if (x <= 0) return x;               // netative values are returned unchanged.
    sx = 0;
    a  = x;
    if (x < (1LL << (n_bit / 2 - 1))) { // If x is in the lower half of the digits,
        sx = n_bit / 4;                 // note that square root halves the digits.
        a  = x = (x << (n_bit / 2));    // enlarge to the upper half of the didits.
    }
    a >>= n_room;                       // scale down to get head-room
    sr  = n_bit / 2 - n_room;
    t = 1LL << (n_bit - n_room - 2);
    while (a < t) {                     // find operating point
        t >>=2;
        sr--;
    }
    a >>= sr;                           // now, a / t is in (0.25, 1)
    t = 1LL << (sr + n_room);           // scaled '1'
    a = (22 * a + 11 * t) / 32;         // -3.1e-2 < err < +2.9e-2, a / t in (0.25, 1)
    a = (x / a + a) >> 1;               // 0 <= err < 4.8e-4
    a = (x / a + a) >> 1;               // 0 <= err < 1.2e-7
    a = (x / a + a) >> 1;               // 0 <= err < 6.5e-15
    if ((x - a) > (a * a))              // the approximation error is always zero or negative
        a++;                            // adjust LSB so that the error is within 1/2 LSB
    if (sx) {
        a += (1LL << (sx - 1)) - 1;     // for round half
        a >>= sx;
    }
    return a;
}
#ifdef DEBUG
#include <stdio.h>
#include <math.h>
int test_sqrt(int64_t x) {
    int64_t y;
    double  ref;
    ref = round(sqrt((double)x));
    y = sqrt64(x);
    fprintf(stdout, "sqrt64(%"PRId64") = %"PRId64" %s %.f\n",
        x, y, y != ref ? "!=" : "==", ref);
    return (y != ref ? 1 : 0);
}
int main(int argc, char *argv[]) {
    int     i;
    int64_t x, y;
    if (argc == 2) {
        sscanf(&argv[1][0], "%" SCNd64 "\n", &x);
        y = sqrt64(x);
        fprintf(stdout, "sqrt64(%"PRId64") = %"PRId64"\n", x, y);
    } else {
        x = 1LL;
        for (i = 0; i < 63; i++) {
            test_sqrt(1LL << i);
            test_sqrt(x);
            x = (x << 1) | 1;
        }
    }
    return 0;
}
#endif

※ INT64_MAX (0x7fffffffffffffff) の平方根は、0xb504f334 (= 3037000500 ) となって、 INT32_MAX (0x7fffffff) よりも大きくなるため、戻り値を int32_t にすることは出来ません。

Ｑ＆Ａでみる整数の平方根の設計と実装

Q1: Halley 法は 3 次収束で速いのに、なぜ使わないの？

A1: 整数演算では、途中で値が 大きくなりすぎる からです。

平方根に対する Halley 法の更新式の一例は次の通りです。

x ← x (x + 3 a / x ) / ( 3 x + a / x )
この式では、反復の途中で x (x + 3 a / x) ≈ 4a のように、一時的に「元の値 a と同程度」の大きさの数が現れます。
浮動小数点演算では気になりませんが、固定ビット幅の整数演算では、この一時的な膨張が致命的です。

⛔ Halley 法は速いが、整数ではダイナミックレンジを消費しすぎる。

Q2. 初期近似を 2 次にすれば、Newton の回数を減らせませんか？

A2. 初期誤差は減りますが、代わりに中間量が肥大化します。

例えば、区間 x ∈ (0.25, 1) における 2 次の Min-Max 近似として、

g(x) = (66 + 270x - 81x2) / 256
を用いると、初期近似誤差は約 0.6 % まで下げられます。 これは 1 次近似より 約 0.7 桁良い精度です。
しかし、この式には x2 の項が含まれています。 整数演算では、これが一時的に 倍以上のビット幅を要求します。

⛔ 近似次数を上げると、ダイナミックレンジが一気に増大する。

Q3. そういうの、ワークラウンドで回避できませんか？

A3. 可能ですが、その時点で負け。本末転倒です。

といった回避策はあります。 しかし、 条件分岐は計算のパイプラインを乱して遅くなるなどの見えにくい弊害もありますし、 システマティックな誤差解析も困難になります。つまり、それらを入れるくらいなら
Newton 法をもう 1 回回したほうが、速くて安全
になります。

Newton 法の更新式は、

x ← (x + a/x) / 2
のように、反復中ずっと √a オーダーの値だけを扱います。

✅ ニュートン法は速くはないが、スリムで強くて挙動も読める。

Q4. 結局、何を優先したのですか？

A4. 値域全体を、安全に、説明可能な形でカバーすることです。

本実装では、

という構成を採っています。

結果として、

という、
✅ 整数演算として最も堅牢な設計になっています。

まとめ

浮動小数点では「速い方法」が、
整数演算では「危ない方法」になることがある。
この実装では、ダイナミックレンジを最優先に設計判断を行った。
✅ 付随的ですが、数値的信頼性も高く、シンプルで移植性も高いものとなっています。
www.finetune.co.jp [Mail] © 2000 Takayuki HOSODA.