JavaのBigIntegerを知る2
前回「JavaのBigIntegerを知る1」Javaの標準ライブラリが持つ多倍長整数の実装、BigIntegerのデータ構造やコンストラクタをおさえたので、今回はBigIntegerを使用した四則演算について見てみます。
一口に四則演算といっても、内部では計算に使用する桁数等に応じて種々のアルゴリズムが選択されていたりするのですが、ここでは一番オーソドックスなものを取りあげることにしようと思います。
環境
- OS: Arch Linux (linux kernel: 4.8.13)
- Java: openjdk version “1.8.0_121”
1. BigIntegerの加算
加算に限らず四則演算全般そうですが、多倍長整数の計算を考えるときは筆算でどのように行うかを考えるのがヒントになります。
加算の場合、例えば 261 + 143 という計算を考えてみると、筆算の手順は
1 + 3 = 4で一の位は46 + 4 = 10で十の位は0かつ繰り上がり発生2 + 1 + 1 = 4で百の位は4- よって和は404
のようになります。
これと同じことをやっているのが BigInteger#add(int[], int[]) です。
private static int[] add(int[] x, int[] y) {
// If x is shorter, swap the two arrays
if (x.length < y.length) {
int[] tmp = x;
x = y;
y = tmp;
}
int xIndex = x.length;
int yIndex = y.length;
int result[] = new int[xIndex];
long sum = 0;
if (yIndex == 1) {
// yが一桁ならば下からの繰り上がりを考えずにただ和をとればよい
sum = (x[--xIndex] & LONG_MASK) + (y[0] & LONG_MASK) ;
result[xIndex] = (int)sum;
} else {
// Add common parts of both numbers
while (yIndex > 0) {
// (1) 下の位からの繰り上がりを考慮しながら二数の同じ桁の和をとる
sum = (x[--xIndex] & LONG_MASK) +
(y[--yIndex] & LONG_MASK) + (sum >>> 32);
// (2) (int) でcastした値がこの桁に収まる値となる
result[xIndex] = (int)sum;
}
}
// (3) Copy remainder of longer number while carry propagation is required
boolean carry = (sum >>> 32 != 0);
while (xIndex > 0 && carry)
carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
// Copy remainder of longer number
while (xIndex > 0)
result[--xIndex] = x[xIndex];
// (4) Grow result if necessary
if (carry) {
int bigger[] = new int[result.length + 1];
System.arraycopy(result, 0, bigger, 1, result.length);
bigger[0] = 0x01;
return bigger;
}
return result;
}
一番重要な処理はコメントの(1), (2)の部分です。ここで二数に共通な各桁の計算と繰り上がりの処理を行っています。
後半の処理では、 (3) yとの加算は終わったが繰り上がりが発生し続ける場合(e.g. 10進数の計算でいう19999+1)や、 (4) もとの値より桁数が増える場合(e.g. 10進数の計算でいう99+1)を考慮しながら最終的な和を計算します。
実際のBigInteger同士の加算ではいきなりこのメソッドが呼び出されるわけではなく、まずはじめに BigInteger#add(BigInteger val) メソッド内で二数の符号を確認し、その符号が同じ場合にこの加算処理が行われるようになっています。 もし符号が異なる場合は、後述の減算処理が行われます。
public BigInteger add(BigInteger val) {
if (val.signum == 0)
return this;
if (signum == 0)
return val;
// 符号が同じなので加算処理
if (val.signum == signum)
return new BigInteger(add(mag, val.mag), signum);
int cmp = compareMagnitude(val);
if (cmp == 0)
return ZERO;
// 符号が異なるので減算処理
int[] resultMag = (cmp > 0 ? subtract(mag, val.mag)
: subtract(val.mag, mag));
resultMag = trustedStripLeadingZeroInts(resultMag);
return new BigInteger(resultMag, cmp == signum ? 1 : -1);
}
2. BigIntegerの減算
BigInteger#subtractも行うことは加算と似ています。
加算の場合は各桁を繰り上がりを考慮しながら足し合わせましたが、減算の場合は各桁を借り入れを考慮しながら減じます。
private static int[] subtract(int[] big, int[] little) {
int bigIndex = big.length;
int result[] = new int[bigIndex];
int littleIndex = little.length;
long difference = 0;
// Subtract common parts of both numbers
while (littleIndex > 0) {
// (1) 上の位からの借り入れを考慮しながら二数の差をとる
difference = (big[--bigIndex] & LONG_MASK) -
(little[--littleIndex] & LONG_MASK) +
(difference >> 32);
// (int) でcastした値がこの桁に収まる値となる
result[bigIndex] = (int)difference;
}
// Subtract remainder of longer number while borrow propagates
boolean borrow = (difference >> 32 != 0);
while (bigIndex > 0 && borrow)
borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
// Copy remainder of longer number
while (bigIndex > 0)
result[--bigIndex] = big[bigIndex];
return result;
}
一番重要な処理(1), (2)もプラスがマイナスになった程度でかなり加算と似ていますが、減算の場合これでなぜ借り入れが実現できているのか少しイメージし辛いです。
結論を先に言うと、2の補数表現を考えればlong同士の計算で(2の補数上)正の数になろうが負の数になろうが、その下位32bitをそのまま差としてうまく扱うことができます。
正の数の場合はもちろんそのまま差として使えます (e.g. 10進数で5-3は実際の値も2の補数で考えてももちろん2)。
負の数の場合を具体的な例を使って、例えばひと桁4bit(基数が16)の場合で考えてみます。
この場合 3 - 5 を考えると、bit上では 3 = (0011), 5 = (0101) で 3 - 5 = -2 = (1110) と表されます。
(1110) は2の補数表現とせずそのまま見れば14であり、つまり 3 - 5 の計算時に上の桁から1借りて基数である16を足した結果に等しくなっています。
これは補数の性質から来るものであり、基数がいくつでも成り立つ話なので、このやり方で上の桁からの借り入れを考慮した計算を行うことができます。
(ちなみにある基数bのある数cの桁数がnのとき、b進数におけるcの補数は b^n-c と定義される。言葉でいえば、cにもう一つ桁を増やすために足す必要のある最小の数、という感じ?。詳細は補数 - Wikipedia)
また借り入れが発生した場合、その桁の計算結果は(2の補数表現上)負の数になるため、計算結果を表すdifferenceの上位32bitは 0xffffffff となっています。
一方で借り入れが発生しなければ正の数なので上位32bitは 0x00000000 です (元々はintであったもの同士の引き算なので上位32bitを使用することにはならない)。
ということで次の桁の計算時にlongとしてこの上位32bitの値をたせば、借り入れが発生した場合には -1 を足し、発生しなければ 0 を足すということになり借り入れをうまく表現することができます。
3. BigIntegerの乗算
乗算では使用する数の桁数に応じて異なるアルゴリズムが選択されるようですが、ここでは一番単純な実装のBigInteger#multiplyToLenを見てみます。
このアルゴリズムではこれまで同様筆算をイメージしたものになります。
// x, yが乗算に使用され、zには結果が入る
private static int[] multiplyToLen(int[] x, int xlen, int[] y, int ylen, int[] z) {
int xstart = xlen - 1;
int ystart = ylen - 1;
if (z == null || z.length < (xlen+ ylen))
z = new int[xlen+ylen];
long carry = 0;
for (int j=ystart, k=ystart+1+xstart; j >= 0; j--, k--) {
long product = (y[j] & LONG_MASK) *
(x[xstart] & LONG_MASK) + carry;
z[k] = (int)product;
carry = product >>> 32;
}
z[xstart] = (int)carry;
for (int i = xstart-1; i >= 0; i--) {
carry = 0;
for (int j=ystart, k=ystart+1+i; j >= 0; j--, k--) {
// (1) 下位の桁からの繰り上がりを考慮したlong値の乗算
long product = (y[j] & LONG_MASK) *
(x[i] & LONG_MASK) +
(z[k] & LONG_MASK) + carry;
// (2) 積のセット
z[k] = (int)product;
// (3) 繰り上がりの計算
carry = product >>> 32;
}
z[i] = (int)carry;
}
return z;
}
(1), (2), (3)で筆算と同様の処理を行います。 加算、減算を確認したあとならば素直に実装が理解できると思います。
ちなみにわざわざ最初に一回別にforを回しているのは、xが一桁の場合にこれのみで終わらせられてパフォーマンスが上がるからなのかと思うのですがはっきりとはわかりません。 単純に実装するだけならば最初のfor文を無くして後ろのものと統合しても結果は同じはずなので…
4. BigIntegerの除算
最後にBigIntegerの除算アルゴリズムについて見てみるのですが、その前に除算処理で頻出するMutableBigIntegerクラスを確認しておきます。
class MutableBigInteger {
int[] value;
int intLen;
int offset = 0;
...
}
MutableBigIntegerはBigIntegerと同様多倍長整数を表すクラスですが、名前の通りmutableであり、またoffsetやintLenによりvalue配列の一部を使って数を表すことができます。 パッケージプライベートなクラスなので外部に提供されているクラスではなく、主にBigIntegerやBigDecimalが行う一部の計算で性能を上げるために使用されるクラスです。
それでは改めて除算の処理を見ていきます。 考え方はこれまでと同じで、手で筆算を使って行うことを何とかコンピュータアルゴリズムに変換すればよいのですが、除算に関しては他と違い中々一筋縄ではいきません。
そこでまずはわかりやすくするために除数がひと桁である場合を見てみることにします。
例えば10進数でいう 512 / 5 のような計算になるのですが、この場合はほぼ筆算と同じノリで考えられます。
MutableBigIntegerではこの場合の処理がdivideOneWordメソッドで実装されています。
// thisが被除数、divisotrが除数、quotientに商を入れる
int divideOneWord(int divisor, MutableBigInteger quotient) {
... (省略)
// (1) 被除数の最上位の桁に関する計算
int rem = value[offset];
long remLong = rem & LONG_MASK;
if (remLong < divisorLong) {
quotient.value[0] = 0;
} else {
quotient.value[0] = (int)(remLong / divisorLong);
rem = (int) (remLong - (quotient.value[0] * divisorLong));
remLong = rem & LONG_MASK;
}
// (2) 被除数の上位から2番目の桁からひとつずつ計算
int xlen = intLen;
while (--xlen > 0) {
// (2-1) 一つ上の桁の剰余をもらい現在の桁との和を算出する
long dividendEstimate = (remLong << 32) |
(value[offset + intLen - xlen] & LONG_MASK);
int q;
// (2-2) 現在の桁に対する商、剰余の計算
if (dividendEstimate >= 0) {
// (2-2a) 単純にlong値の計算で商と剰余を計算
q = (int) (dividendEstimate / divisorLong);
rem = (int) (dividendEstimate - q * divisorLong);
} else {
// (2-2b) dividentEstimateが最上位のbitを使用している場合
// long値として単純に計算をしてしまうと結果が負になる
// そのため両辺を2で割ったあと調整して答えを出す
long tmp = divWord(dividendEstimate, divisor);
q = (int) (tmp & LONG_MASK);
rem = (int) (tmp >>> 32);
}
// qをこの桁の商としquotientにセットする
quotient.value[intLen - xlen] = q;
// remLongは次の桁の計算に使用する
remLong = rem & LONG_MASK;
}
... (省略)
}
実装では最上位(1)とそれ以外の桁(2)に関する計算で処理をわけているので少し長く見えますが、本質的なところでは同じようなことをやっています。 ここでは(1)を無視して(2)に注目します。 (2)の処理にいくつかコメントを加えたので筆算を思い出しながら見ていただくと何をやっているかイメージできるのではないかと思います。
(2-2)の分岐は最初見たときは戸惑いましたが、恐らくJavaでlong値を普通に計算してしまうと2の補数表現とみなされて負の値になってしまうために、最上位bitを使用する値が被除数になるときは別途一手間かけて除算を行っているのかなと思います。 なおここで使用しているdivWordの実装は以下のようになっています。 このあと見るアルゴリズムでもそうですが、まず商を概算して、それから正しい値となるように調整するというのが一つの除算の方針になるみたいです。
static long divWord(long n, int d) {
long dLong = d & LONG_MASK;
long r;
long q;
if (dLong == 1) {
q = (int)n;
r = 0;
return (r << 32) | (q & LONG_MASK);
}
// Approximate the quotient and remainder
// 被除数、除数それぞれに右シフト処理することで最上位の桁を0にする
// そうすることでlong値として商、剰余が求められる
q = (n >>> 1) / (dLong >>> 1);
r = n - q*dLong;
// Correct the approximation
// 上で右シフトをかましたことで本来の計算結果と若干異なる可能性がある
// その補正を以下の2処理で行う
while (r < 0) {
r += dLong;
q--;
}
while (r >= dLong) {
r -= dLong;
q++;
}
// n - q*dlong == r && 0 <= r <dLong, hence we're done.
return (r << 32) | (q & LONG_MASK);
}
簡単な除算の感覚を掴んだところで、次に多倍長整数同士の除算を見てみるために、MutableBigInteger#divideMagnitudeに注目します。 実際には被除数の桁数に応じて異なるアルゴリズムが選択されるようですが、恐らくこちらの実装が一番多倍長整数の除算としては標準的なものかなと思います。 (実装のコメントを見る限りKnuthさんの「The Art of Computer Programming Vol.2 4.3.1 D」のアルゴリズムが使用されているとのことですし)
このアルゴリズムに関しては実装だけではとても何をやっているかはわからないと思うので、アルゴリズムの詳細を知りたい場合には上記「The Art of Computer Programming Vol.2」を確認したほうがよいと思います。 少なくとも私は本を読むまでまるで何もわかりませんでした。
考え方の肝だけでもここで触れておくと、筆算のときに商の各桁を上から順番に求めていくと思いますが、その計算値が除数の1桁目と被除数の現在の2桁からほぼ確定できるということがこのアルゴリズムのポイントになっています。 これにより本来の除数がいくつだろうと、上で触れた除数が1桁のときのアルゴリズムで計算を行うことができるようになります。この証明全体はまだ追えていないのですが、はじめて読んだときはむちゃくちゃすごいなと思いました。
とはいえ上記本でアルゴリズムの動作を理解しても、実際の実装は難しいところやわからないところが多々あったので、そういったところをいくつかpickupしておこうと思います。
D1
まず以下がアルゴリズム上D1(正規化)に当たる部分のコードになります。 日本語のコメントは私が勝手に加えたものです。
private MutableBigInteger divideMagnitude(MutableBigInteger div,
MutableBigInteger quotient,
boolean needRemainder ) {
// D1 normalize the divisor
// (D1-1) 本で言うdの値を計算するためにまずshiftを求める
int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
// Copy divisor value to protect divisor
final int dlen = div.intLen;
int[] divisor;
MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
if (shift > 0) {
// (D1-2-A) shift>0 ということは除数の最上位がb/2より小さいということになる
// その場合、2^shiftを被除数、除数にそれぞれかける
divisor = new int[dlen];
copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
int[] remarr = new int[intLen + 1];
rem = new MutableBigInteger(remarr);
rem.intLen = intLen;
rem.offset = 1;
copyAndShift(value,offset,intLen,remarr,1,shift);
} else {
// この場合、単純に2^shiftをかけてしまうと被除数はoverflowしてしまう
// なので桁を一つ増やしたremarrを用意し、それをremにセットする
int[] remarr = new int[intLen + 2];
rem = new MutableBigInteger(remarr);
rem.intLen = intLen+1;
rem.offset = 1;
int rFrom = offset;
int c=0;
int n2 = 32 - shift;
for (int i=1; i < intLen+1; i++,rFrom++) {
int b = c;
c = value[rFrom];
remarr[i] = (b << shift) | (c >>> n2);
}
remarr[intLen+1] = c << shift;
}
} else {
// (D1-2-B) すでに除数の最上位がb/2より大きいので何もせず、単純にdivisorをdivからコピーする
divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
rem = new MutableBigInteger(new int[intLen + 1]);
System.arraycopy(value, offset, rem.value, 1, intLen);
rem.intLen = intLen;
rem.offset = 1;
}
int nlen = rem.intLen;
// Set the quotient size
final int limit = nlen - dlen + 1;
if (quotient.value.length < limit) {
quotient.value = new int[limit];
quotient.offset = 0;
}
quotient.intLen = limit;
int[] q = quotient.value;
// Must insert leading 0 in rem if its length did not change
if (rem.intLen == nlen) {
rem.offset = 0;
rem.value[0] = 0;
rem.intLen++;
}
int dh = divisor[0];
long dhLong = dh & LONG_MASK;
int dl = divisor[1];
...
}
正規化の目的は除数の最上位桁の値が b/2 (この場合2^31)より大きくなるように除数、被除数に適切な値dをかけることです。これはアルゴリズムを機能させるために必要な条件になります。
ここでは(D1-1)で取得したshiftを使用して 2^shift をdとして取るようにしています。
Integer.numberOfLeadingZeros は与えられたintを上位から見てはじめに1が出てくるまでの桁数をカウントしたもので、例えば8の場合、0x00000008 なので28となります。
もしshiftの値が0の場合、これは除数の最上位桁の最上位bit (わかりにくいですが) が1であることを意味し、つまりb/2 よりすでに大きい値であることを意味します。
なのでこの場合は(D1-2-B)のように正規化を行わずに処理を続けることができます。
一方でshift値が0より大きい場合、除数の最上位桁をb/2より大きくするためには 2^shift をかけるのですが、このときshiftは除数をもとに決定しているだけなので、被除数にも同じ値をかけた場合にoverflowを起こす可能性があります。
これを避けるために実装では被除数も Integer.numberOfLeadingZeros を計算し、overflowを起こすことがわかればひと桁多く取り直したMutableBigIntegerを用意し、今後の処理に使用することにしています。
なお、これ以降の処理では、divisorが元々のdivをコピーして正規化した除数、quotientが商、remが(はじめは正規化された被除数で、適宜更新されながら最終的に)剰余となる値をそれぞれ表します。
D2, D3
// D2 Initialize j
// 商の各桁を上位から計算する
for (int j=0; j < limit-1; j++) {
// D3 Calculate qhat
// estimate qhat
int qhat = 0; // 商の推測値
int qrem = 0; // 剰余の推測値
boolean skipCorrection = false;
int nh = rem.value[j+rem.offset];
int nh2 = nh + 0x80000000;
int nm = rem.value[j+1+rem.offset];
// (D3-1) qhat, qremの計算
if (nh == dh) {
// (D3-1-A) qhatの計算値が基数b-1を超えてしまう場合の計算?
// この場合qhatはb-1にする?
// qremの設定の仕方等わからないところが多い
qhat = ~0;
qrem = nh + nm;
skipCorrection = qrem + 0x80000000 < nh2;
} else {
// (D3-1-B) 通常の場合の計算
long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
if (nChunk >= 0) {
qhat = (int) (nChunk / dhLong);
qrem = (int) (nChunk - (qhat * dhLong));
} else {
long tmp = divWord(nChunk, dh);
qhat = (int) (tmp & LONG_MASK);
qrem = (int) (tmp >>> 32);
}
}
if (qhat == 0)
continue;
// (D3-2) 推測値を本来の値に近づけるためにいくつか処理が行える
// (詳細はThe Art of Computer Programming参照)
if (!skipCorrection) { // Correct qhat
long nl = rem.value[j+2+rem.offset] & LONG_MASK;
long rs = ((qrem & LONG_MASK) << 32) | nl;
long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
if (unsignedLongCompare(estProduct, rs)) {
qhat--;
qrem = (int)((qrem & LONG_MASK) + dhLong);
if ((qrem & LONG_MASK) >= dhLong) {
estProduct -= (dl & LONG_MASK);
rs = ((qrem & LONG_MASK) << 32) | nl;
if (unsignedLongCompare(estProduct, rs))
qhat--;
}
}
}
...
}
D2が定義するforループの中で、商の各桁を計算し、D3で現在の桁に関する商の推測値(qhat)、剰余の推測値(qrem)を計算します。
(D3-1)でqhat, qremの計算を行っており、通常の場合(D3-1-B)のように被除数の2桁と除数の最上位1桁の商、剰余をそれぞれ設定します。
しかし被除数の上位1桁と除数の最上位1桁が等しい場合、qhatとして基数b-1(2^32-1)となる値をとるようになっています(D3-1-A)。 ただこのときの実装はわからない点が多くいまいちピンときていないです。
(D3-2)では推測値の修正を行います。 詳細は本を参照して頂きたいのですが、この処理をはさむことでqhatをかなり高確率で本来の値にとることができるようです。
D4
D4では計算したqhatを使用してremの更新を行います。
// D4 Multiply and subtract
rem.value[j+rem.offset] = 0;
int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
実際の計算はmultsubで行われています。 この処理は実際は多倍長整数と1桁の整数の乗算と、多倍長整数同士の減算が組み合わさったものになります。
private int mulsub(int[] q, int[] a, int x, int len, int offset) {
long xLong = x & LONG_MASK;
long carry = 0;
offset += len;
for (int j=len-1; j >= 0; j--) {
long product = (a[j] & LONG_MASK) * xLong + carry;
long difference = q[offset] - product;
q[offset--] = (int)difference;
carry = (product >>> 32)
+ (((difference & LONG_MASK) >
(((~(int)product) & LONG_MASK))) ? 1:0);
}
return (int)carry;
}
D5
ごくまれにD4で計算した値が負になる場合があり、そのときに被除数divisorを除数として使っているremに足し戻します。
というのがD5の目的なはずですが、実装からそのニュアンスが伝わらずうーんという感じです。borrowは上で計算したmulsubの返り値、nh2は現在着目している被除数の1桁目と 0x80000000 の和をとったものなのですが…
// D5 Test remainder
if (borrow + 0x80000000 > nh2) {
// D6 Add back
divadd(divisor, rem.value, j+1+rem.offset);
qhat--;
}
divaddは単純な多倍長整数の加算です。
private int divadd(int[] a, int[] result, int offset) {
long carry = 0;
for (int j=a.length-1; j >= 0; j--) {
long sum = (a[j] & LONG_MASK) +
(result[j+offset] & LONG_MASK) + carry;
result[j+offset] = (int)sum;
carry = sum >>> 32;
}
return (int)carry;
}
D7
qhatを求めた商の桁に設定し、次の桁の計算にうつります。
という感じで、この先実装では最終桁に対する処理が繰り返されますが、本質的にはD3-D6と同様の処理になるので省略します。
正直除算の実装はわからない点もいくつかあって歯がゆい感じもしますが、処理の流れは掴めたと思うのでとりあえずはここまでで除算の処理を追うのはやめておこうと思います。わからないところは春休みの宿題です (ないけど)。
Summary
- 多倍長整数の四則演算の標準的な実装は筆算をイメージするとわかりやすい
- 除算に関しては計算できる値から商を予想し本来の結果へ調整するのが一つの方針