符号なし 256 ビット整数の被除数左ビットシフト除算

C言語で記述されたビットシフト筆算法の除算のコードの大半は、除数を右シフトして被除数から引き算していく書き方になっているようです。例えば、古い GNU C のランタイムの整数除算などがそうなっています。これは筆算のやりかたをそのまま手順化してあるわけですが、私が高校生の時分に 8 ビットマイコン機械語で除算を書いていた頃は、除数ではなく、余り・被除数・商を合わせて左シフトしながら除数を引いていき、被除数が左シフトして空いたビットスペースに商を埋めていく書き方をしていたものです。

単純な例として、32 ビットマシンで、256 ビットの整数除算を被除数左ビットシフトでやってみます。ゼロで割ろうとすると、商と余りを最大値で埋めて false を返します。それ以外では、商と余りを求めて true を返します。

なお、この手のビットシフト除算は、論理回路を作製するのが容易で、ハードウェアで使われることが多いのですが、最も遅い除算アルゴリズムなので、ソフトウェアによる 256 ビットの整数の除算に使うには向いていません。筆算除算の系統でも、桁ごとに掛け算をして引き算する スクールブック除算 の方がワードのビット長さ分、速度を稼ぐことができます。ちなみに、新しい GNU C のランタイム倍精度整数除算では、高速除算法の一つ Burnikel の再帰除算を使っています。

#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>

/*       (numer[0]+nemer[1]*B+..+numer[7]*B**7)
 *     = (denom[0]+denom[1]*B+..+denom[7]*B**7)
 *     * (q[0]+q[1]*B+..+q[7]*B**7) + (r[0]+r[1]*B+..+r[7]*B**7)
 *     where B = 2**32
 */
bool
uint256_divmod (uint32_t const * const numer, uint32_t const * const denom,
                  uint32_t * const q, uint32_t * const r)
{
    int i;
    uint32_t c;
    uint32_t x;
    uint32_t y;
    uint32_t w[8];

    if (denom[0] == 0 && denom[1] == 0 && denom[2] == 0 && denom[3] == 0
            && denom[4] == 0 && denom[5] == 0 && denom[6] == 0 && denom[7] == 0) {
        q[0] = q[1] = q[2] = q[3] = q[4] = q[5] = q[6] = q[7] = UINT32_MAX;
        r[0] = r[1] = r[2] = r[3] = r[4] = r[5] = r[6] = r[7] = UINT32_MAX;
        return false;
    }
    r[0] = r[1] = r[2] = r[3] = r[4] = r[5] = r[6] = r[7] = 0;
    q[0] = numer[0]; q[1] = numer[1]; q[2] = numer[2]; q[3] = numer[3];
    q[4] = numer[4]; q[5] = numer[5]; q[6] = numer[6]; q[7] = numer[7];
    for (i = 0; i < 256; i++) {
        /* r[7..0]:q[7..0] 合わせて 512 ビットを 1 ビット左シフトする */
        r[7] = (r[7] << 1) + (r[6] >> 31);
        r[6] = (r[6] << 1) + (r[5] >> 31);
        r[5] = (r[5] << 1) + (r[4] >> 31);
        r[4] = (r[4] << 1) + (r[3] >> 31);
        r[3] = (r[3] << 1) + (r[2] >> 31);
        r[2] = (r[2] << 1) + (r[1] >> 31);
        r[1] = (r[1] << 1) + (r[0] >> 31);
        r[0] = (r[0] << 1) + (q[7] >> 31);
        q[7] = (q[7] << 1) + (q[6] >> 31);
        q[6] = (q[6] << 1) + (q[5] >> 31);
        q[5] = (q[5] << 1) + (q[4] >> 31);
        q[4] = (q[4] << 1) + (q[3] >> 31);
        q[3] = (q[3] << 1) + (q[2] >> 31);
        q[2] = (q[2] << 1) + (q[1] >> 31);
        q[1] = (q[1] << 1) + (q[0] >> 31);
        q[0] = (q[0] << 1);
        /* c:w[7..0] = r[7..0] - denom[7..0] */
        x = r[0]; y = denom[0];                  c  = (x < y); w[0] = x - y;
        x = r[1]; y = denom[1] + c; c = (y < c); c += (x < y); w[1] = x - y;
        x = r[2]; y = denom[2] + c; c = (y < c); c += (x < y); w[2] = x - y;
        x = r[3]; y = denom[3] + c; c = (y < c); c += (x < y); w[3] = x - y;
        x = r[4]; y = denom[4] + c; c = (y < c); c += (x < y); w[4] = x - y;
        x = r[5]; y = denom[5] + c; c = (y < c); c += (x < y); w[5] = x - y;
        x = r[6]; y = denom[6] + c; c = (y < c); c += (x < y); w[6] = x - y;
        x = r[7]; y = denom[7] + c; c = (y < c); c += (x < y); w[7] = x - y;
        if (c == 0) {
            /* r[7..0] >= denom[7..0] */
            r[0] = w[0]; r[1] = w[1]; r[2] = w[2]; r[3] = w[3];
            r[4] = w[4]; r[5] = w[5]; r[6] = w[6]; r[7] = w[7];
            q[0] += 1;
        }
    }
    return true;
}