円周率.jp (http://xn--w6q13e505b.jp/method/add.html)

多倍長加減算

多倍長の加減算については,基本的に筆算と同様のことをすればよい。 本ページではプログラムを実現する上での注意点をいくつか書き留める. 多倍長整数多倍長浮動小数とがあり, 処理を分ける必要がある.

多倍長整数

多倍長整数での加減算はそれこそ筆算通りの計算でできる。 気持ち程度の選択肢として、計算途中に正規化するものと最後に正規化するものとがある。 また、符号を別に持つ形で処理を行う場合は両者の符号関係や絶対値の大きさを比較して処理する必要がある一方、 符号を持たない、非負整数を扱う場合には減算結果が負になったらどうするか決めておくか、 必ず負にならないように演算を適用させる必要がある.

\[ A = a_{n_A-1}R^{n_A-1} + a_{n_A-2}R^{n_A-2} + \cdots + a_1R + a_0 \] \[ B = b_{n_B-1}R^{n_B-1} + b_{n_B-2}R^{n_B-2} + \cdots + b_1R + b_0 \] \[ \Rightarrow C = A \pm B = c_{n_C-1}R^{n_C-1} + c_{n_C-2}R^{n_C-2} + \cdots + c_1R + c_0 \]

上記のような計算モデルで計算を行う.

(1-A) 後で正規化を行う

3141592653
589793
 
3141117123146 正規化処理
3142182446

まずは後で正規化を行うバージョンから. 右の計算イメージのように,とりあえず要素ごとの和を求め, 後で正規化を行う. この方式では複数の正規化ルールを取り換える場合に加算ルーチンが 1 つで済む. また,本来やりたいことだけ記述するのでコードの可読性が高い.

2 つの多倍長数 A と B の和を C に代入するプログラムを書く.

Natural UnsignedAdd(const Natural& a, const Natural& b) {
  Natural c(max(a.size(), b.size()));  // C の桁数は A と B の大きい方に合わせる
  for (int64 i = 0; i < c.size(); ++i) {
    Digit dig_a = (i < a.size()) ? a[i] : 0;  // 桁数を超えると 0 にする
    Digit dig_b = (i < b.size()) ? b[i] : 0;
    c[i] = dig_a + dig_b;
  }
  // C を正規化
  Digit carry = 0;
  for (int64 i = 0; i < c.size(); ++i) {
    c[i] = c[i] + carry;
    carry = c[i] / R;
    c[i] = c[i] % R;
  }
  if (carry)
    c.push_back(carry);  // 最後に繰り上がりがあれば桁数を増やす
  return c;
}

このプログラムでは $R$ を少し超えた数でもオーバーフローしない前提で書いているが、例えば 64 ビット計算機で $R=2^{64}$ とするなど、オーバーフローすることが前提となる場合は正規化の部分が少し変わって

Natural UnsignedAdd(const Natural& a, const Natural& b) {
  ... (足し算の部分は同じ) ...
  // C を正規化
  Digit carry = 0;
  for (int64 i = 0; i < c.size(); ++i) {
    Digit next_carry = (c[i] < a[i]) ? 1 : 0;
    c[i] = c[i] + carry;
    if (c[i] < carry)
      next_carry = 1;
    carry = next_carry;
  }
  if (carry)
    c.push_back(carry);
  return c;
}

この実装では繰り上がりが 1 以下になる前提を利用している。 最初の実装でも同様の前提を利用して構わないが、分岐が無いので除算する形にしている。

(1-B) 計算中に正規化を行う

続いて,計算途中で正規化を行うバージョン. こちらはキャッシュ系マシンでキャッシュ外のメモリを舐める回数が少なくなるため, 若干ではあるが高速になる(と思う).

3141592653
589793
 
(1)46
3141592653
589793
 
(1)2446
3141592653
589793
 
(1)182446
3141592653
589793
 
3142182446

この場合のプログラムは以下の通り.

Natural UnsignedAdd(const Natural& a, const Natural& b) {
  Natural c(max(a.size(), b.size()));  // C の桁数は A と B の大きい方に合わせる
  Digit carry = 0;
  for (int64 i = 0; i < c.size(); ++i) {
    Digit dig_a = (i < a.size()) ? a[i] : 0;  // 桁数を超えると 0 にする
    Digit dig_b = (i < b.size()) ? b[i] : 0;
    c[i] = dig_a + dig_b + carry;
    carry = c[i] / R;
    c[i] = c[i] % R;
  }
  if (carrry)
    c.push_back(carry);
  return c;
}
Natural UnsignedAdd(const Natural& a, const Natural& b) {
  Natural c(max(a.size(), b.size()));  // C の桁数は A と B の大きい方に合わせる
  Digit carry = 0;
  for (int64 i = 0; i < c.size(); ++i) {
    Digit dig_a = (i < a.size()) ? a[i] : 0;  // 桁数を超えると 0 にする
    Digit dig_b = (i < b.size()) ? b[i] : 0;
    Digit dig_c = dig_a + carry;

    carry = 0;
    if (dig_c < dig_a)
      carry = 1;
    dig_c = dig_c + dig_b;
    if (dig_c < dig_b)
      carry = 1;
  }
  if (carry)
    c.push_back(carry);
  return c;
}

(2) 符号を考慮した処理

符号のデータを別に持っている場合,加算と減算の処理を符号だけではなく分ける必要がある。

Integer SignedAdd(const Integer& a, const Integer& b) {
  if (a.sign() != b.sign()) {  // 符号が違う場合は符号を合わせて減算
    return SignedSubtract(a, -b);
  } else {  // 同符号なので演算
    Integer c(UnsignedAdd(a, b));
    c.sign() = a.sign();
    return c;
  }
}

Integer SignedSubtract(const Integer& a, const Integer& b) {
  if (a.sign() != b.sign()) {  // 符号が違う場合は符号を合わせて加算
    return SignedAdd(a, -b);
  } else if (a.abs() >= b.abs()) {    // |A|>=|B| かつ同符号という前提になるので絶対値レベルでの減算
    Integer c(UnsignedSubtract(a, b));
    c.sign() = a.sign();
    return c;
  } else {  // |B| > |A| では C = A - B = -(B - A) に変形して演算
    return -SignedSubtract(-B, -A);
  }
}

多倍長浮動小数

浮動小数では仮数部の有効桁数と指数部との関係で足し合わせるインデックスにズレが出てくるが, 基本的には整数と違いが無い. 指数が異なる場合には片方の数値の末尾が端数となるが,好きな方法で丸めると良い。