文字列で 10 進任意精度整数の四則演算 その 2 符号なし除算

任意精度整数の除算の筆算を2進数以外の位取り記数法の除算ルーチンとして書くには面倒な理論を理解できてなければなりません。なぜ、面倒になるかというと、商の着目している桁の数の候補が 10 進数なら 10 個あり、それらのうちからどれを選べば良いかを少ない計算手順で求める必要があるからです。ネタ元の D. KnuthThe Art of Computer Programming Vol.2 asin:4756145434 (以下 Knuth1998) は除算のアルゴリズムの元になる理論の根拠を説明するのにページを割いており、練習問題合わせてすべて理解して始めて、プログラミングできるようになります。

結論から言えば、割られる数の着目している桁の上位 3 桁と、割る数の上位 2 桁を使えば、商の候補になる数を求めることができます。この数を求めてしまえば、2 進数のときと同じように、とりあえず掛けて引いてみて、負にならなければその数を採用し、負のときは引きすぎているので割る数を足して一つ戻してその数より一つ小さな数を採用することによって除算を進めることができます。ただし、候補の数がこのようになるためには条件があり、条件を満たすようにあらかじめ前処理で割られる数と割る数に因子をかけておかなければならない場合があります。そのため、2進数以外の任意精度整数の除算ルーチンは、桁ごとの候補の数の絞り込みと、前処理と後処理のつじつまあわせがくっついて長くなります。

この商の各桁の候補を求めるアルゴリズムは、倍精度浮動小数点除算の実装の試みがおこなわれていた 1950 年代に遡り、1960 年に Pope と Stein が発表したものを、Knuth が改良して 1969 年に発表したもののようです。他にもいくつかバリエーションがあり、下の Perl のコーディングの書き方にならうと $u01 / ($v0 + 1) とするやりかたも目にします。

ftp://ftp.cs.wisc.edu/pub/techreports/1969/TR55.pdf

さて、「割る数の上位 2 桁を使えば」とあるからには 2 桁以上でなければならないわけでして、割る数が 1 桁のときはどうするか別に考えないといけません。幸いこちらは簡単な方法で計算をおこなえます。この方法は BASE62 を求めるなどの応用に利用できます。

# 割る数が 1 桁のときの除算
# see [Knuth1998] page.625 4.3.1 Excercise 16 (Short division)
sub _divmod_unsigned_short {
    my($r0, $r1) = @_;
    # ok 0 < $r1 && $r1 < $base;
    use integer;
    my $q = '';
    my $r = 0;
    for my $i (0 .. (length $r0) - 1) {
        $r = $r * $base + (substr $r0, $i, 1);
        my $q0 = $r / $r1;
        $r -= $q0 * $r1;
        $q .= chr $q0 + $ord0;
    }
    return ($q, "$r");    
}

大物は次のとおりです。ネタ元のアルゴリズムでは、被除数の着目している範囲の最上位桁も求めていますが、どうせゼロになるのはわかりきっているので、借りがあるかどうか(その桁の引き算が負になるか)だけを判定して範囲の最上位桁には問答無用でゼロを書き込むようにアレンジしています。

2014年4月17日: r0 と q への値のセットする場所を間違えていたので修正

# 割る数が 2 桁以上のときの除算
# see [Knuth1998] page.272 Algorithm D (Division of nonnegative integers)
sub _divmod_unsigned_large {
    my($r0, $r1) = @_;
    use integer;
    my $m = (length $r0) - (length $r1);
    my $n = length $r1;
    # 商の各桁の候補の数 $qhat を絞り込めるように、$r0 と $r1 を正規化する。
    $r0 = '0' . $r0; # now (length $r0) == $m + $n + 1
    my $d = $base / ((substr $r1, 0, 1) + 1);
    if ($d > 1) {
        # $r0 = _mul($r0, $d);
        my $i = $m + $n + 1;
        my $k = 0;
        while (--$i >= 0) {
            my $a = (substr $r0, $i, 1) * $d + $k;
            $k = $a / $base;
            $a -= $k * $base;
            substr $r0, $i, 1, chr $a + $ord0;
        }
        # $r1 = _mul($r1, $d);
        $i = $n;
        $k = 0;
        while (--$i >= 0) {
            my $a = (substr $r1, $i, 1) * $d + $k;
            $k = $a / $base;
            $a -= $k * $base;
            substr $r1, $i, 1, chr $a + $ord0;
        }
    }
    # 筆算の手法で、桁をずらしながら商と余りを求める。
    my $q = '0' x ($m + 1);
    #  q: 0 1 ..       m
    # r0: 0 1 2 .. n   m  m+1 .. m+n
    # r1:   0 1 .. n-1
    #            --> 桁をずらしていく
    #                    0   .. n-1
    my $v0 = 0 + (substr $r1, 0, 1);
    my $v1 = 0 + (substr $r1, 1, 1);
    for my $j (0 .. $m) {
        # 商の候補の数 $qhat を絞り込む
        my $u0 = 0 + (substr $r0, $j, 1);
        my $u01 = $u0 * $base + (substr $r0, $j + 1, 1);
        my $u2 = $j + 2 >= (length $r0) ? 0 : 0 + (substr $r0, $j + 2, 1);
        # $qhat = min $u01 / $v0, $base - 1; が初期値
        my $qhat = $u01 / $v0;
        if ($qhat > $base - 1) {
            $qhat = $base - 1;
        }
        # アルゴリズム中の $rhat は $u01 - $qhat * $v0 で、このコードではその都度求めることにした
        #      $qhat * $v1 > $rhat * $base + $u2  の比較をしている
        while ($qhat * $v1 > ($u01 - $qhat * $v0) * $base + $u2) {
            --$qhat;
        }
        if ($qhat > 0) {
            # これで $q[$j] は $qhat - 1 か $qhat のどちらかになる
            # とりあえず、 $qhat をかけて引いてみる。
            # (substr $r0, $j+1, $n) = _sub((substr $r0, $j+1, $n), _mul($r1, $qhat));
            my $k1 = 0;
            my $k2 = 1;
            my $i = $n;
            while (--$i >= 0) {
                my $a1 = $qhat * (substr $r1, $i, 1) + $k1;
                $k1 = $a1 / $base;
                $a1 -= $k1 * $base;
                my $a2 = (substr $r0, $j + $i + 1, 1) - $a1 + $k2 + $base - 1;
                $k2 = $a2 >= $base ? 1 : 0;
                $a2 -= $k2 * $base;
                substr $r0, $j + $i + 1, 1, chr $a2 + $ord0;
            }
            if ((substr $r0, $j, 1) - $k1 + $k2 < 1) { # 最上桁の引き算が負になるならば、
                # 引きすぎたので、一つ戻す。
                --$qhat;
                # (substr $r0, $j+1, $n) = _add((substr $r0, $j+1, $n), $r1);
                my $k = 0;
                my $i = $n;
                while (--$i >= 0) {
                    my $a = (substr $r0, $i + $j + 1, 1) + (substr $r1, $i, 1) + $k;
                    $k = $a >= $base ? 1 : 0;
                    $a -= $k * $base;
                    substr $r0, $i + $j + 1, 1, chr $a + $ord0;
                }
            }
        }
        substr $r0, $j, 1, '0'; # 余りで使うため計算済みの桁をゼロで潰す。
        substr $q, $j, 1, chr $qhat + $ord0; # 商に求めた数をセットする。
    }
    # 余りを、正規化のときにかけた因子で割ってつじつまをあわせをする。
    if ($d > 1) {
        # (substr $r0, $m+1, $n) = _div((substr $r0, $m+1, $n), $d);
        my $k = 0;
        for my $i (0 .. $n - 1) {
            my $a = $k * $base + (substr $r0, $m + $i + 1, 1);
            my $x = $a / $d;
            $k = $a - $x * $d;
            substr $r0, $m + $i + 1, 1, chr $x + $ord0;
        }
    }
    return ($q, $r0);
}