Go to Qiita Advent Calendar 2024 Top
LoginSignup
4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@ikiuo

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

Posted at

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

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

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

割り算の除数を $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
4
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?