開平法（二進数）による平方根

筆算の開平法の数値を、１０進数から２進数に置き換えただけです。二進数の開平法では、乗除算が不要です。

現代では、なかなかお目にかかれないであろう非力なプロセッサ（演算命令は整数のみで除算がないとか…）で平方根を計算する場合、二進数の開平法が選択肢の１つになります。

サンプル プログラムは整数部で計算をやめていますが、少し弄れば小数部も求められます。

sample.cpp

#include <iostream>

using namespace std;

typedef int32_t I32;
typedef uint32_t U32;
typedef uint64_t U64;

struct sqrt_res
{
    U32 root; // 平方根
    U32 rem;  // 余り
};

// 開平法による符号なし整数(32bit)の平方根
sqrt_res sqrt(U32 x)
{
    sqrt_res res;

    int shift = 32/2; // 残り回数
    U32 root = 0; // 平方根
    U32 sub = 0;  // 副算用
    U32 rem = 0;  // 余り

    // 上位の "0" をスキップして詰める
    if ((x >> 0x10) == 0) { x <<= 0x10; shift -= 16/2; }
    if ((x >> 0x18) == 0) { x <<= 0x08; shift -=  8/2; }
    if ((x >> 0x1c) == 0) { x <<= 0x04; shift -=  4/2; }
    if ((x >> 0x1e) == 0) { x <<= 0x02; shift -=  2/2; }

    // 開平法で1ビットずつ求める
    do
    {
        // x から2ビット取り出す
        rem = (rem << 2) | (x >> (32 - 2));
        x <<= 2;

        // "0" の場合で初期化
        root <<= 1;
        sub  <<= 1;

        // 判定処理
        I32 chk = rem - sub - 1;
        if (chk >= 0)
        {
            // "1" の場合で更新
            root += 1;
            sub  += 2;
            rem = chk;
        }
    }
    while (--shift > 0);

    res.root = root;
    res.rem  = rem;
    return res;
}

// メイン関数
int main(int argc, char **argv)
{
    if (argc != 2)
    {
        cout << "Usage: " << argv[0] << " N" << endl;
        return 1;
    }

    U64 x = strtoul(argv[1], nullptr, 10);
    if (errno == ERANGE)
    {
        cout << "error: range: " << argv[1] << endl;
        return 2;
    }
    if (x >= (1ULL<<32))
    {
        cout << "error: range: " << argv[1]
             << " >= " << (1ULL<<32)
             << endl;
        return 2;
    }

    sqrt_res res = sqrt(x);
    U32 y = res.root * res.root + res.rem;

    if (x != y)
        cout << "Bug:" << x
             << " != " << y
             << endl;
    cout << "N = " << x
         << " = " << res.root
         << " * " << res.root
         << " + " << res.rem
         << endl;
    return 0;
}

実行結果

$ ./sample 127806
N = 127806 = 357 * 357 + 357

例 127806 を２進数の筆算（開平法）で平方根を求めると

127806 : 1_1111_0011_0011_1110
   357 :           1_0110_0101

   1  0  1  1  0  0  1  0  1
----------------------------
) 01 11 11 00 11 00 11 11 10     1
   1                             1
  -----                         ---
     11                         100
      0                           0
     -----                      ----
     11 11                      1001
     10 01                         1
     --------                   -----
      1 10 00                   10101
      1 01 01                       1
      ----------                ------
           11 11                101100
               0                     0
           --------             -------
           11 11 00             1011000
                  0                   0
           -----------          --------
           11 11 00 11          10110001
           10 11 00 01                 1
           --------------       ---------
            1 00 00 10 11       101100100
                        0               0
            ----------------    ----------
            1 00 00 10 11 10    1011001001
              10 11 00 10 01             1
            ----------------    ----------
               1 01 10 01 01    1011001010