備忘録：除算命令を使わない除算

除算命令がないプロセッサ向け。

一般的に筆算の方法を２進数で行う。

試しにニュートン法も実装してみる。

割り算の除数を $D$ として、関数 $f \left( x \right)$ を

f \left( x \right) = \frac{1}{x} - D

とすると、$f \left( x \right) = 0$ では

\begin{eqnarray}
f \left( x \right) & = & \frac{1}{x} - D = 0 \\
x & = & \frac{1}{D}
\end{eqnarray}

に近づく漸化式

\begin{eqnarray}
x_{i+1} & = & x_i - \frac{ f \left( x_i \right) }{ f' \left( x_i \right) } \\
& = & x_i - \frac{ x_i^{-1} - D }{ - x_i^{-2} } \\
& = & x_i + \left( x_i - D\ x_i^2 \right) \\
& = & x_i \left( 2 - D\ x_i \right)
\end{eqnarray}

を使う。

サンプル プログラム

sample.cpp

#include <iostream>
#include <iomanip>

using namespace std;

// 分岐予測用
#ifndef likely
#if defined(__GNUC__)
#define likely(b) __builtin_expect(!!(b), 1)
#else
#define likely(b) b
#endif
#endif
#ifndef unlikely
#if defined(__GNUC__)
#define unlikely(b) __builtin_expect(!!(b), 0)
#else
#define unlikely(b) b
#endif
#endif

//　型宣言

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

struct udiv_t
{
    U32 quot;
    U32 rem;
};

// フラグ

static bool debug = false;

// 上位ビットから連続する "0" を数える
#if defined(__GNUC__)
#define CLZ(x) __builtin_clz(x)
#else
static int CLZ(U32 N)
{
    int r = 0;
    if ((N >> 0x10) == 0) { r += 16; N <<= 16; }
    if ((N >> 0x18) == 0) { r +=  8; N <<=  8; }
    if ((N >> 0x1c) == 0) { r +=  4; N <<=  4; }
    if ((N >> 0x1e) == 0) { r +=  2; N <<=  2; }
    if ((N >> 0x1f) == 0) { r +=  1; N <<=  1; }
    if (N == 0) r += 1;
    return r;
}
#endif

// 筆算（二進数）の除算 : 符号なし整数(32bit)
udiv_t udiv_bn(U32 N, U32 D)
{
    if (unlikely(N < D))
    {// 被除数が除数より小さい
        udiv_t res;
        res.quot = 0;
        res.rem  = N;
        return res;
    }

    if (unlikely(N == D))
    {// 被除数と除数が等しい
        udiv_t res;
        res.quot = 1;
        res.rem  = 0;
        return res;
    }

    int Zd = CLZ(D);
    int Zn = CLZ(N);

    if (likely(D != (D & -D)))
    {// 除数が 2 の乗数ではない
        int S = Zd - Zn;
        U32 Q = 0;
        U32 R = N;

        D <<= S;
        do
        {
            Q <<= 1;
            if (R >= D)
            {
                R -= D;
                Q += 1;
            }
            D >>= 1;
        }
        while (--S >= 0);

        udiv_t res;
        res.quot = Q;
        res.rem  = R;
        return res;
    }

    if (likely(D != 0))
    {// 除数が 0 ではない
        udiv_t res;
        res.quot = N >> (31 - Zd);
        res.rem  = N & (D - 1);
        return res;
    }

    {// 除数が 0
        udiv_t res;
        res.quot = 0;
        res.rem  = N;
        return res;
    }
}

// ニュートン法による除算 : 符号なし整数(32bit)
udiv_t udiv_nr(U32 N, U32 D)
{
    if (unlikely(N < D))
    {// 被除数が除数より小さい
        udiv_t res;
        res.quot = 0;
        res.rem  = N;
        return res;
    }

    if (unlikely(N == D))
    {// 被除数と除数が等しい
        udiv_t res;
        res.quot = 1;
        res.rem  = 0;
        return res;
    }

    int Zd = CLZ(D);

    if (likely(D != (D & -D)))
    {// 除数が 2 の乗数ではない

        // 整数部 32bit、小数部 32bit の固定小数点演算を使って
        // ニュートン法による除数の逆数を求める

        U64 TWO = 2ULL << 32; // 2.0

        // 初期値: x[0]
#if NO_TABLE
        // テーブル参照なし
        // ループは 5,6 回程度
        U32 xn = 1 << Zd;
#else
        // テーブル参照あり
        // ループは 1〜4 回程度
        // 除数の先頭 4bit を使用
        static U32 table[8] = {
            (U32)((8ULL << 32) / (16 +  1)),
            (U32)((8ULL << 32) / (16 +  3)),
            (U32)((8ULL << 32) / (16 +  5)),
            (U32)((8ULL << 32) / (16 +  7)),
            (U32)((8ULL << 32) / (16 +  9)),
            (U32)((8ULL << 32) / (16 + 11)),
            (U32)((8ULL << 32) / (16 + 13)),
            (U32)((8ULL << 32) / (16 + 15)),
        };
        int Tn = (D << Zd) >> (31 - 3);
        U32 xn = table[Tn - 8] >> (30 - Zd);
#endif
        U32 xp = 0;

        do
        {// x[n+1] = x[n] (2 - D x[n])
            xp = xn;
            xn = ((TWO - ((U64)D * xp)) * xp) >> 32;
#if DEBUG
            if (debug)
                cout << showbase
                     << "xp=" << hex << xp
                     << ",xn=" << hex << xn
                     << dec << endl;
#endif
        }
        while (xp != xn);

        // 商と余を求める
        U32 Q = ((U64)N * xn) >> 32;
        U32 R = N - (Q * D);

        // 演算誤差の補正
        if (R >= D)
        {
            Q += 1;
            R -= D;
        }

        udiv_t res;
        res.quot = Q;
        res.rem  = R;
        return res;
    }

    if (likely(D != 0))
    {// 除数が 0 ではない
        udiv_t res;
        res.quot = N >> (31 - Zd);
        res.rem  = N & (D - 1);
        return res;
    }

    {// 除数が 0
        udiv_t res;
        res.quot = 0;
        res.rem  = N;
        return res;
    }
}

// メイン関数
int main(int argc, char **argv)
{
    const char *program = argv[0];

    if (argc >= 2)
    {
        if (argv[1][0] == '-' &&
            argv[1][1] == 'd' &&
            argv[1][2] == 0)
        {
            debug = true;
            --argc;
            ++argv;
        }
    }

    if (argc != 3)
    {
        cout << "Usage: " << program << " N D" << endl;
        return 1;
    }

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

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

    U32 N = (U32)Na;
    U32 D = (U32)Da;
    U32 Q = (D != 0) ? (N / D) : 0;
    U32 R = (D != 0) ? (N % D) : N;

    udiv_t res;

    res = udiv_bn(N, D);
    if (res.quot != Q || res.rem != R)
        cout << "Bug! udiv_bn: " << Q << "," << R << endl;
    cout << "udiv_bn: " << N
         << " = " << D
         << " * " << res.quot
         << " + " << res.rem
         << endl;

    res = udiv_nr(N, D);
    if (res.quot != Q || res.rem != R)
        cout << "Bug! udiv_nr : " << Q << "," << R << endl;
    cout << "udiv_nr: " << N
         << " = " << D
         << " * " << res.quot
         << " + " << res.rem
         << endl;

    return 0;
}

gcc version 8.3.0 (Raspbian 8.3.0-6+rpi1) でコンパイルして udiv_nr 部分を抜粋

udiv_nr(unsigned, unsigned):
        cmp     r1, r2
        bcc     .L36
        beq     .L37
        rsb     ip, r2, #0
        and     ip, ip, r2
        cmp     ip, r2
        push    {r4, r5, r6, r7, r8, r9, lr}
        clz     r4, r2
        beq     .L24
        lsl     r3, r2, r4
        ldr     ip, .L38
        lsr     r3, r3, #28
        sub     r3, r3, #8
        rsb     r5, r4, #30
        ldr     r4, [ip, r3, lsl #2]
        mov     lr, #2
        mov     ip, #0
        lsr     r4, r4, r5
        b       .L25
.L28:
        mov     r4, r3            ; ループ先頭
.L25:
        umull   r8, r9, r2, r4
        subs    r6, ip, r8
        sbc     r3, lr, r9
        umull   r6, r7, r6, r4
        mla     r3, r4, r3, r7
        cmp     r4, r3
        bne     .L28              ; ループ終端

        umull   r4, r5, r1, r4
        mov     ip, r5
        mul     r3, r5, r2
        sub     r1, r1, r3
        cmp     r2, r1
        addls   ip, r5, #1
        subls   r1, r1, r2
        str     ip, [r0]
        str     r1, [r0, #4]
        pop     {r4, r5, r6, r7, r8, r9, pc}
.L37:
        mov     r2, #1
        mov     r3, #0
        strd    r2, [r0]
        bx      lr
.L24:
        cmp     ip, #0
        beq     .L27
        rsb     r3, r4, #31
        sub     ip, ip, #1
        lsr     r3, r1, r3
        and     ip, ip, r1
        stm     r0, {r3, ip}
        pop     {r4, r5, r6, r7, r8, r9, pc}
.L27:
        mov     r3, #0
        str     r1, [r0, #4]
        str     r3, [r0]
        pop     {r4, r5, r6, r7, r8, r9, pc}
.L36:
        mov     r3, #0
        str     r1, [r0, #4]
        str     r3, [r0]
        bx      lr
.L39:
        .align  2
.L38:
        .word   .LANCHOR0


        .section        .rodata
        .align  2
        .set    .LANCHOR0,. + 0
        .type   _ZZ7udiv_nrjjE5table, %object
        .size   _ZZ7udiv_nrjjE5table, 32
_ZZ7udiv_nrjjE5table:
        .word   2021161080
        .word   1808407282
        .word   1636178017
        .word   1493901668
        .word   1374389534
        .word   1272582902
        .word   1184818564
        .word   1108378657
        .bss