先日のABCで中国剰余定理による連立合同式が出たことでまた冷えたのでこれをきっかけに中国剰余定理を勉強しました。

勉強にはこちらのけんちょんさんの中国剰余定理の解説記事を使いましたが、いくつか理解が詰まった箇所があるので自分のメモも兼ねた補足記事を書くことにしました(本人了承済み)。けんちょんさんの記事を読みながらこれを読むことを想定しています。

書いたCRT関連のライブラリ

視聴者プレゼントです。Verifyにはけんちょんさんの記事の下にある問題たちを使いました。MITライセンスで、参考するなどお好きにしてください。

使い方は大体日本語のコメントとして残しています。基本的にACLのと使い方は同じになるように作りました。

これより下はけんちょんさんの記事の補足になります。

2 中国剰余定理のアルゴリズム

2-2 拡張 Euclid の互除法を用いた方法

該当する箇所

記事の通りに、二元のやり方さえわかれば、2つを選んで1つにまとめるを繰り返すことでn本の連立方程式でもできます。

ところで、今 $$ x _ i := 1からi番目のまでの式を連立させたときの解 (\mod m _ 1 m _ 2 \cdots m _ i) $$ とすると、$ i + 1 $番の式を連立させるとき、以下のような条件を満たすことになる。 $$ x _ {i + 1} \equiv x _ {i} (\mod m _ 1 m _ 2 \cdots m _ i) \ x _ {i + 1} \equiv b _ {i + 1} (\mod m _ {i + 1}) $$

これらから、$ k, l \in 整数, x _ {i + 1} = k (m _ 1 m _ 2 \cdots m _ i) + x _ i = l m _ {i + 1} + x _ {i + 1} $を満たすとわかる。

これを式変形すると、 $$ k (m _ 1 m _ 2 \cdots m _ i) - l m _ {i + 1} = x _ {i + 1} - x _ i $$

という形となり、これは$ x _ {i + 1} - x _ i $が$ \gcd(m _ 1 \cdots m _ i m _ {i + 1}) $で割り切れるなら、拡張ユークリッドで係数$ k, l $を求めることができます。

よって、$ x _ {i + 1} = x _ i + k (m _ 1 m _ 2 \cdots m _ i) = l m _ {i + 1} + x _ {i + 1} $を代入すればいいが、けんちょんさんの記事では$ k $に代入して計算している。なぜでしょうか。

拡張ユークリッドで$ ax + by = \gcd(a, b) $を求めるとき、$ (x, y) $は共に$ -|b| < x < |b|, -|a| < y < |a| $となることが保証されます。これを$ k, l $について考えると、

$ k $は$ (-m _ {i + 1}, m _ {i + 1}) $の間に存在。
$ l $は$ (-m _ 1 \cdots m _ i, m _ 1 \cdots m _ i) $の間に存在。

ここで、明らかに$ l $の方が存在する範囲が広いとわかります。範囲が広い場合、64bit整数でもオーバーフローすることによって実質的に計算できなくなることがあります。特に$ 10 ^ 6 $をレベルを3つ連立する場合は$ l $を計算する場合はオーバーフローします。(原理的には$ k $も$ l $も掛ければ同じぐらいのオーダーになるが、細かく見ると$ k $の方が微妙にオーバーフローしにくい。)

これを使うところ：中華風(Easy)

最後に、計算量を考えてみます。拡張ユークリッド互除法の計算量は最悪で見積もって、$ O(\log (m _ 1 \cdots m _ n)) $となります。よって、全体では、$ O(n \log (m _ 1 \cdots m _ n)) $となります。

このアルゴリズム、欠点としてはオーバーフローが挙げられる。$ n $が多くなると、$ m _ 1 \cdots m _ n $は明らかに64bit整数では収まらない場合が多く、結果的にはこのアルゴリズムは大量の連立方程式を解くのは(常に多倍長演算が必須という意味では)向いていないと言えます。 この常に多倍長整数での演算が必須という欠点を解消したのが、Garner's Algorithmである。これは、計算量におおよそ$ n $がもう1つつくが、その代わりに任意modで求めたり、多倍長整数の保存に必要なメモリを大幅に削減することができるアルゴリズムです。

Garnerのアルゴリズム

該当する箇所

けんちょんさんの記事では(勝手ながらも)かなりごちゃごちゃしていて見づらいので、自分が考えたのを載せます。やってることは同じです。(コードに起こすと少し煩雑かもしれない)

$$ x _ i := 式1からiまで統合した時のxの値 $$ とします。この時、以下のように考えることができます。拡張ユークリッドの互除法では、2つの合同式から倍数関係が分かる$ ax + by = \gcd(a, b) $の式を求めましたが、こちらでは1つの等式から倍数関係で$ x $を表し、それをもう1つの合同式に代入している形です。

なお、$ m _ 1, \cdots m _ n $はお互い互いに素であるとします。(そうでなくとも、指定のアルゴリズムでそのようにすることができます。後述)

$$ x _ {i + 1} \equiv x _ i (\mod m _ 1 \cdots m _ n), k \in 整数, x _ {i + 1} = x _ i + k m _ 1 \cdots m _ n $$

これを$ x _ {i + 1} \equiv b _ {i + 1} (\mod m _ {i + 1}) $に代入します(ここが拡張ユークリッドの互除法で$ ax + by = \gcd(a, b) $と違うところ)。

$$ x _ i + k m _ 1 \cdots m _ n \equiv b _ {i + 1} (\mod m _ {i + 1}) $$ $$ k m _ 1 \cdots m _ n \equiv b _ {i + 1} - x _ i (\mod m _ {i + 1}) $$ すると、$ m _ 1, \cdots m _ n $はお互いに互いに素なので、$ m _ 1, \cdots m _ i $はすべてmod$ m _ {i + 1} $に対して、逆元を持ちます。 よって、 $$ k \equiv (b _ {i + 1} - x _ i) (m _ 1 \cdots m _ i) ^ {-1} $$ と、$ k $を求めることができます。

よって、$ x _ {i + 1} $は$ k $を代入することで求めることができます。また、$ x _ {i + 1} $を求めるときの$ k=p _ {i + 1} $とする($ p _ {i} $は存在しない)と、最終的な$ x _ n $は以下のようになる。

$$ x _ n = b _ 1 + m _ 1 p _ 2 + m _ 1 m _ 2 p _ 3 + \cdots + (m _ 1 \cdots m _ {n - 1})p _ n $$ $$ p _ i = (b _ i - x _ {i - 1})(m _ 1 \cdots m _ {i-1}) ^ {-1} $$ ただし、2つ目の式の逆元は$ \mod m _ i $である。

この式を(時間計算量を犠牲にしても、できるだけ多倍長整数を使わずに値を保持するように)計算するのがGarner's Algorithmである。

実際の自分の実装(一般的と思われる)では、以下にように計算を行った。

求める最終的な$ x _ n $はMODで割る値であるとする。
$ i $が小さい順に、$ x _ i $を計算していく。
$ x _ i $をMODで割ったもの(コード内ではx_MOD)を保持。
$ m _ 1 \cdots m _ {i - 1} $をMODで割ったもの(コード内ではm_mul_MOD)を保持。
$ p _ i $の値は、MODで割らずに配列で保持する。 拡張ユークリッドの互除法により、$ p _ i $自体の値の範囲は$ (-m _ i, m _ i) $にあるため、多倍長を使わずに保持可能。
$ x _ {i - 1} $の値は、$ \mod m _ {i} $では保持できないので、毎回$ O(i) $で、modを取ってない$ {t _ j} $や$ (m _ 1 \cdots m _ {i - 1}) ^ {-1} $で計算し直す。

拡張ユークリッドの互除法を用いた前の方法と比べて、毎回ループの中で_余計に$ O(n) $の計算量が懸かっているので、全体では$ O(n(n + \log(\max({m _ i})) $程度になります。)

問題例

4-2 中国剰余定理を用いる応用問題たち

該当する箇所はここ。

ここでは、$ m _ 1, \cdots m _ n $がお互いに互いに素になるようにするアルゴリズム(written by snuke)の解説や動作例を挙げてみる。

けんちょんさんの記事にあるように、$ m _ i, m _ j $が互いに素ではないとき、共通する因数をいい感じに割り振ることで、$ lcm(m _ i, m _ j) $をそのままに、$ \gcd(m _ i, m _ j)=1 $とすることができます。 割り振り方はけんちょんさんのをご覧ください。

割り振りの一例をあげてみます。$ 120, 84 $の場合は、

$$ 120 = 2 ^ 3 \times 3 \times 5 $$ $$ 252 = 2 ^ 2 \times 3 ^ 2 \times 7 $$ であるので、

$$ 40 = 2 ^ 3 \times 5 $$ $$ 63 = 3 ^ 2 \times 7 $$

とすることで、互いに素を達成できます。

この時、わざわざ素因数分解して計算するのも問題はないのですが、難解ですがsnukeさんのアルゴリズムを採用することで、gcdのみを使用して素因数の割り振りが可能となります。(そしてこっちの方が早い)

基本的には、各素因数に関して、肩の指数がより多い方に割り振ります(先ほどの例では、120は$ 2 ^ 3 $で、252は$ 2 ^ 2 $であるので、120のほうに素因数2を割り振り、252の方の$ 2 ^ 2 $を消します。)

では、snukeさんのアルゴリズムを解説します。アルゴリズムは以下です。

互いに素ではない$ m _ i, m _ j $を見つけます(全部で計$ O(n ^ 2) $の計算量をかけてみることになります) $ g = \gcd(m _ i, m _ j) $とします。
そもそも解なしの場合を判定してはじきます
$ (m _ i)' \leftarrow m _ i / g, (m _ j)' \leftarrow m _ j / g $と代入します。$ g $で割ることによって、共通してる因数を削る(最後に適宜かけて正しい$ m _ i, m _ j $に戻す)
$ gi \leftarrow \gcd((m _ i)', g) $この操作は、$ m _ i $から$ m _ j $との共通因数を取り除いたものと、$ g $のGCDを求めてることになる。これが1ではないとき、$ g=\gcd(m _ i, m _ j) $に含まれてる因数で、まだ何かが$ (mi)' $に残ってるということである。その残る因数たちで次の$ gj $にあるものを、6~7によって、$ gi $にかき集める。
$ gj \leftarrow g / gi $。$ gi $に含まれてない$ g $の因数を$ gj $に置く。
$ \gcd(gi, gj) $が1でない限り、ずっと7. を行う。
$ h = \gcd(gi, gj), gi \leftarrow gi \times h, gj \leftarrow gi / h $ と更新。これは、$ gj $にある$ gi $にあるべき因数を$ gi $を移す作業に当たる。
最後に、$ m _ i \leftarrow (m _ i)' \times mi, m _ j \leftarrow (m _ j)' \times mj $とする。

実際のコードはけんちょんさんの記事にありますが、以下にも載せておきます。

bool preGarner(vector<ll>& val, vector<ll>& mod) {
    //計算量O(n^2 * log(max{mod[i]}))
    int n = val.size();

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < i; j++) {
            ll g = gcd(mod[i], mod[j]);
            if ((val[j] - val[i]) % g != 0)
                return false;

            //ここの部分はブログで詳しく例を挙げて紹介している。
            mod[i] /= g, mod[j] /= g;
            ll gi = gcd(mod[i], g);
            ll gj = g / gi;

            //ちなみにこのループは一度では終わらない。
            do{
                ll g = gcd(gi, gj);
                gi *= g, gj /= g;
            } while (gcd(gi, gj) != 1);

            mod[i] *= gi, mod[j] *= gj;

            //val[i], val[j]は小さくなったmod[i], mod[j]に合わせて再計算する。
            val[i] %= mod[i], val[j] %= mod[j];

        }
    }
    return true;
}

これは実際どのように動くのかをシミュレーションしてみる。

シミュレーション1

$$ m _ i = 2 ^ 5 \times 3 ^ 2 \times 5 $$ $$ m _ j = 2 ^ 2 \times 3 ^ 3 \times 5 $$

この時、$ g = 2 ^ 2 \times 3 ^ 2 \times 5 $となるので、

$$ (m _ i)' = m _ i / g = 2 ^ 3 $$ $$ (m _ j)' = m _ j / g = 3 $$

$$ mi = \gcd( (m _ i)', g) = \gcd(2 ^ 3, g) = 2 ^ 2 $$ $$ mj = g / 2 ^ 2 = 3 ^ 2 \times 5 $$

これは、$ \gcd(mi, mj) = 1 $であるので特に特にステップ6~7を実行することがなく(すでに約数の振り分けが済んでいるパターン)、

$$ m _ i = (m _ i)' \times mi = 2 ^ 3 \times 2 ^ 2 = 2 ^ 5 $$ $$ m _ j = (m _ j)' \times mj = 3 \times (3 ^ 2 \times 5) = 3 ^ 3 \times 5 $$

シミュレーション2

今度は$ m _ i, m _ j $を入れ替えてみます。

$$ m _ j = 2 ^ 2 \times 3 ^ 3 \times 5 $$ $$ m _ j = 2 ^ 5 \times 3 ^ 2 \times 5 $$

この時、やはり$ g = 2 ^ 2 \times 3 ^ 2 \times 5 $となるので、

$$ (m _ i)' = m _ j / g = 3 $$ $$ (m _ j)' = m _ i / g = 2 ^ 3 $$ $$ mi = \gcd( (m _ i)', g) = \gcd(3, g) = 3 $$ $$ mj = g / 3 = 2 ^ 2 \times 3 \times 5 $$

これは、$ \gcd(mi, mj) = 3 $であるので、$ m _ i $側に振り分けるべき因数の3がまだ$ mj $にあることから、$ mi \leftarrow mi \times 3, mj \leftarrow mj / 3 $とします。その結果以下のようになります。

$$ mi = 3 ^ 2 $$ $$ mj = 2 ^ 2 \times 5 $$

$$ m _ i = (m _ i)' \times mi = 3 \times 3 ^ 2 = 3 ^ 3 $$ $$ m _ j = (m _ j)' \times mj = 2 ^ 3 \times (2 ^ 2 \times 5) = 3 ^ 3 \times 5 $$

snukeさんのアルゴリズムでGCDのみを使用して、片方に存在している因数を巧妙に振り分けられているとわかります。

これらを受けて、Garner’s Algorithmは時間計算量$ O(n ^ 2 \log)(\max{m _ i})) $となります。

最後に

以上で行間を埋める解説は終了します。ここまで一読ありがとうございました。いずれこれを作らなければならなかったのでABCで冷えてよかったと思います。(小並)

Special Thanks

(敬称略)

この記事を作るにあたって

37zigen　Garner's Algorithmの理解に大いに手伝っていただきました。
けんちょん　記事のsnukeさんのアルゴリズム部分の理解に大いに手伝っていただきました。
ikefumy　拡張ユークリッド互除法の実を使う方法とGarner's Algorithmの違いの検討をしてくれました。

他にも多くの方に助けられています。ありがとうございました！谢谢诸位无私帮助！

Senの競技プログラミング備忘録

こけた問題を自分用の解説で載せる。けんちょんさんのブログを目指したい。質的にも量的にも。こけた問題だけに限定するけど

けんちょんさんの中国剰余定理の補足