🌐English

概要

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

ハードウェア化を念頭に、2-bit 単位、すなわち radix-4 での漸近近似を行っています。 全ビットの近似が完了した時点で LSB を調整し、誤差を ±1/2 LSB に収めています。

64-bit の符号なし整数の最小誤差丸め平方根の最大値は 2³² になるため、戻り値には 33-bit が必要となります。

アルゴリズム

ここで A_n は radix-4 において、 その時点で評価に用いる被開平数の prefix（上位ビット群）を表します。
コード中では an = a >> (4 * i - 4) が A_n に対応します。

(xs << 4) + ((x << 3) + k) * k

(xs << 4) + ((x * k) << 3) + (k * k)

「上位 2bit ずつ取り込みながら、常に平方が prefix 以下になる最大の桁を選ぶ」

LSB 補正

if (a - x > xs)
    x++;

if (a - x > xs)

floor(sqrt(a) + 0.5)

誤差分布

誤差のヒストグラムは、例えば 38-bit システムにおける 0 〜 2³⁸ - 1（全 2³⁸ 個、最大値 274877906943 = 0x3fffffffff） の全数検査では、次のようになる。

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

総サンプル数 2³⁸ 個に対して、最近接整数丸めにおける誤差区間幅が 1/4, 1/2, 1/4 LSB であることに対応し、1:2:1 の分布になっていることが確認できる。

C 実装

sqrt64q -- 与えられた 64-bit 符号なし整数の平方根に対して誤差が最小となる整数を返す。

/* sqrt64q - calculate the square root of a given unsigned integer
 * Rev.1.0 (2026-01-12) (c) 2026 Takayuki HOSODA
 * SPDX-License-Identifier: BSD-3-Clause
 */
#include <inttypes.h>
uint64_t sqrt64(uint64_t a);
uint64_t sqrt64(uint64_t a) {
    uint64_t x       = 0;        // current approximation
    uint64_t xs      = 0;        // square of the current approximation
    uint64_t k          ;        // trial digit in radix-4 (0..3)
    uint64_t an         ;        // A_n: current prefix of the radicand
    for (int i = 2 * sizeof(a); i > 0; i--) {
        if ((an = a >> (4 * i - 4))) {
            k    = (((xs << 4) +   x * 24  + 9) <= an) ? 3 :
                   (((xs << 4) +   x * 16  + 4) <= an) ? 2 :
                   (((xs << 4) +   x *  8  + 1) <= an) ? 1 : 0;
            xs   =   (xs << 4) + ((x << 3) + k) * k;
            x    =   (x  << 2) + k;
        }                        // skip if upper digits are zero
    }
    if (a - x > xs)              // residual = a - x^2 is in [0, 2x]
        x++;                     // round to the nearest integer sqrt
    return x;
}

Note:
UINT64_MAX (0xffffffffffffffff) の平方根は、0x100000000 (= 4294967296) となって、 UINT32_MAX (0xffffffff) よりも大きくなるため、戻り値を uint32_t にすることは出来ません。

ハードウェア実装観点での最適化

k の決定ロジック

k = (… ? 3 : … ? 2 : … ? 1 : 0);

• 16x2 は共通
• 8x は共通
• k2 は LUT

乗算器の扱い

• x * 3, x * 2, x * 1
• k * k（k は 2-bit）

• シフト + 加算
• 2-bit × 2-bit の小さな LUT

xs の保持

平方を再計算せずに xs を保持しているのは、 このクラスのアルゴリズムでは最も重要な高速化です。

1bit-radix だと減算型が多いですが、 radix-4 で xs を持つ方式は比較器中心の設計になり、 高速化（高周波数化）しやすい構成です。

まとめ

このアルゴリズムは「数値演算」ではなく「比較問題」です。

16x² + 8xk + k² ≤ A

を 3 回比較しているだけで、これは CPU でも極めて計算コストの低い命令列になります。

Newton–Raphson 法は FPU や乗除算器を使用しますが、こちらは整数 ALU と比較器だけで完結します。

ハードウェア化においても、各ステップで必要なのは：

• 64-bit 加算器（xs << 4 + …）
• 64-bit 比較器 × 3
• シフタ
• k 用の小 LUT（2-bit）

しかも 32 step で必ず終了（64-bit 版の場合）。
これは、典型的な CORDIC 型・digit recurrence 型 の回路で、

• 高速・低遅延
• 低消費電力
• 小面積
• 完全に deterministic

と、実にハードウェア化に適した構造となっています。

REFERENCES

• M. D. Ercegovac, T. Lang, Digital Arithmetic, Morgan Kaufmann, 2004. (Chapter 10: Digit-Recurrence Square Root)
• P. Markstein, IA-64 and Elementary Functions, Prentice Hall, 2000. (Square Root and Rounding Analysis)

SEE ALSO

• 整数の平方根アルゴリズム（sqrt64, Newton–Raphson method）
• 整数の自乗アルゴリズム（sqr32, Karatsuba-like method）