GMPの有理数mpq_tの加減算を，実装を読みながら理解する． 加算も減算も本質的には同じなので、ここでは単に加算と呼ぶ。
op1が\(\frac{a}{b}\), op2が\(\frac{c}{d}\)で表され, それぞれ既約分数であるとする． これらの和\(op1 + op2 = \frac{ad + bc}{bd}\)が既約分数であるとは限らないので，計算結果をいい感じに既約分数で表したい． ようは，通分して足し算した後に，必要なら約分したい．という話．
愚直な実装法（読み飛ばしてOK）
op1の分子分母にop2の分母を，op2の分子分母にop1の分母を掛けて通分し，分子の加算を行う．さいごに，分子分母をGCD(ad+bc, bd)で割る．
GMPの実装
NUM(op)はopの分子要素を，DEN(op)はopの分母要素を取り出すマクロ定義である． aorc.cにmpz_aorsという関数があり，引数としてmpz_addかmpz_subのどちらかを渡して加減算を切り替えている．aorsは多分"add or subtraction"の略．
GMPの実装でポイントになるのは，
- divexactを使う
割り切れることが保証されているときの除算は速くなる（らしい，あんまりわかってない）． Exact Division (GNU MP 6.2.1) - 乗算をなるべく後ろに持っていく
桁を落としてから掛け算を行う．
の2点だと思う
mpq/aors.c
mpq_aors (mpq_ptr rop, mpq_srcptr op1, mpq_srcptr op2,
void (*fun) (mpz_ptr, mpz_srcptr, mpz_srcptr))
{
...
}
void
mpq_add (mpq_ptr rop, mpq_srcptr op1, mpq_srcptr op2)
{
mpq_aors (rop, op1, op2, mpz_add);
}
void
mpq_sub (mpq_ptr rop, mpq_srcptr op1, mpq_srcptr op2)
{
mpq_aors (rop, op1, op2, mpz_sub);
}
以下にmpq_aorsの実装の，主要な部分を抜粋して記載する．
mpz_gcd (gcd, DEN(op1), DEN(op2));
if (! MPZ_EQUAL_1_P (gcd))
{
mpz_t t;
MPZ_TMP_INIT (t, MAX (op1_num_size + op2_den_size,
op2_num_size + op1_den_size) + 2 - SIZ(gcd));
mpz_divexact_gcd (t, DEN(op2), gcd);
mpz_divexact_gcd (tmp2, DEN(op1), gcd);
mpz_mul (tmp1, NUM(op1), t);
mpz_mul (t, NUM(op2), tmp2);
(*fun) (t, tmp1, t); // _add か _subを実行
mpz_gcd (gcd, t, gcd);
if (MPZ_EQUAL_1_P (gcd))
{
mpz_set (NUM(rop), t);
mpz_mul (DEN(rop), DEN(op2), tmp2);
}
else
{
mpz_divexact_gcd (NUM(rop), t, gcd);
mpz_divexact_gcd (tmp1, DEN(op2), gcd);
mpz_mul (DEN(rop), tmp1, tmp2);
}
}
else
{
/* The common divisor is 1. This is the case (for random input) with
probability 6/(pi**2), which is about 60.8%. */
mpz_mul (tmp1, NUM(op1), DEN(op2));
mpz_mul (tmp2, NUM(op2), DEN(op1));
(*fun) (NUM(rop), tmp1, tmp2);
mpz_mul (DEN(rop), DEN(op1), DEN(op2));
}
まず最初に，分母b, dの最大公約数を求めている．
\(gcd := GCD(b, d)\)
b, dが互いに素である
b, dが互いに素のとき， ad + bcとbdが互いに素である確率はリーマンゼータ関数で求めることができ，これは約60%である．そこそこ高い確率で既約分数となってくれるので， GCD(分子, 分母)の計算と，約分の処理を省くことができる．
b, dが互いに素でない
b, dが互いに素でなければ，\(GCD(ad + bc, bd) \neq 1\) であるため，約分が確実に必要となる．
\(t := d / gcd, tmp2 := b / gcd\)
最大公約数での除算は，確実に割り切れるということが保証されているのでmpz_divexact_gcdという関数を使うことができる．
\(tmp1 \leftarrow a \cdot t = \frac{ad}{gcd} \)
\(t \leftarrow c \cdot tmp2 = \frac{bc}{gcd} \)
\(t \leftarrow tmp1 \pm t = \frac{ad \pm bc}{gcd} \)
次に，
\(gcd \leftarrow GCD(t, gcd) = GCD(\frac{ad \pm bc}{GCD(b, d)} , GCD(b, d)) \)
を計算して，gcdが1かどうかで処理を切り替える．
\(gcd = 1\)なら，
\(NUM(rop) \leftarrow t = \frac{ad \pm bc}{GCD(b, d)}\) // 分子が確定
\(DEN(rop) \leftarrow d \cdot tmp2 = \frac{bd}{GCD(b, d)} \) // 分母が確定\(gcd \neq 1\)なら，
\(NUM(rop) \leftarrow t/gcd = \frac{ \frac{ad \pm bc}{GCD(b, d)} }{ GCD(\frac{ad \pm bc}{GCD(b, d)} , GCD(b, d)) } \) // 分子が確定
\(tmp1 \leftarrow d/gcd = \frac{d}{ GCD(\frac{ad \pm bc}{GCD(b, d)} , GCD(b, d)) }\)
\(DEN(rop) \leftarrow tmp1 \cdot tmp2 =\frac{ \frac{bd}{GCD(b, d)} }{ GCD(\frac{ad \pm bc}{GCD(b, d)} , GCD(b, d)) }\) // 分母が確定
正当性
ぜんぶで3パターンあるが，いずれも \(\frac{ad + bc}{bd}\)になっていることが確認できる．
