mpq_tの乗算 - ぺんぎんさんのおうち

GMPの有理数mpq_tの乗算を，実装を読みながら理解する．

op1が$\frac{a}{b}$, op2が$\frac{c}{d}$で表され，既約分数であるとする．
これらの積$op1 \cdot op2 = \frac{ac}{bd}$が既約分数であるとは限らないので，計算結果をいい感じに既約分数で表したい．約分せずに放置すると，無駄なメモリ領域を使用するおそれがある．

愚直な実装法（読み飛ばしてOK）

$\frac{ac}{bd}$を既約分数で表すには，分子分母を$GCD(ac, bd)$で割ればいい．

例) $\frac{3}{8} \cdot \frac{4}{9} = \frac{12}{72} = \frac{12}{72} \cdot \frac{ \frac{1}{12} }{ \frac{1}{12}} = \frac{1}{6} $

この方法の問題点は，乗算で桁数が増えてから除算を行うことである． GMPでは，先に約分してから乗算するように実装されている．分数の乗算を手計算でやるときと同じかな．

GMPの実装

実装の重要な部分のみを示す．残りはop1=op2の処理(2乗)やメモリ確保である．もしかしたら今後書くかも．

NUM(op)はopの分子要素を，DEN(op)はopの分母要素を取り出すマクロ定義である．

mpq/mul.c

  mpz_gcd (gcd1, NUM(op1), DEN(op2));
  mpz_gcd (gcd2, NUM(op2), DEN(op1));

  mpz_divexact_gcd (tmp1, NUM(op1), gcd1);
  mpz_divexact_gcd (tmp2, NUM(op2), gcd2);

  mpz_mul (NUM(prod), tmp1, tmp2);

  mpz_divexact_gcd (tmp1, DEN(op2), gcd1);
  mpz_divexact_gcd (tmp2, DEN(op1), gcd2);

  mpz_mul (DEN(prod), tmp1, tmp2);

$ t := GCD(a, d), \, s := GCD(b, c)$ を計算しておく．
$tmp1 := \frac{a}{t}, \, tmp2 := \frac{c}{s} $と除算してから， $NUM(prod) := tmp1 * tmp2$．

分母も同様に， $t, s$を使って，
$tmp1 := \frac{d}{t}, \, tmp2 := \frac{b}{s} $と除算してから， $DEN(prod) := tmp1 * tmp2$．

正当性

$$ \frac{NUM(prod)}{DEN(prod)} = \frac{ \frac{a}{t} \frac{c}{s} }{ \frac{d}{t} \frac{b}{s} } = \frac{ac}{bd} \cdot \frac{\frac{1}{ts}}{\frac{1}{ts}} = \frac{ac}{bd} $$

OK．

除算で桁数を落としてから乗算するので，メモリを無駄遣いせずに済むし，計算量も少なく済む．

ついでに2乗も

実装は次の通り．

  if (op1 == op2)
    {
      /* No need for any GCDs when squaring. */
      mpz_mul (mpq_numref (prod), mpq_numref (op1), mpq_numref (op1));
      mpz_mul (mpq_denref (prod), mpq_denref (op1), mpq_denref (op1));
      return;
    }

op1が既約分数であることが保証されているので，2乗しても既約分数である．
これ，「値が同じ」じゃなくて「ポインタが同じ」のときに，処理を2乗に切り替えてるな．「ポインタは違うけど値が同じな有理数なんてそうそうないやろｗ」という思想なんだろうか． mpz_mulを確認したら，こっちもポインタが同じなときにsqr関数呼び出してた．なるほどなあ． 2乗したいときはプログラマがmpz_sqrを呼べということか．