■「中国剰余定理って知ってる？」
●「先輩いきなりですね。中国剰余定理は、連立合同式の解に関する定理、ですよね」
中国剰余定理：
$$m_1,m_2,\dots,m_n$$
が、どの二つも互いに素であれば、
$$x\equiv r_1\ \ ({\rm mod}\ m_1)$$
$$x\equiv r_2\ \ ({\rm mod}\ m_2)$$
$$\vdots$$
$$x\equiv r_n\ \ ({\rm mod}\ m_n)$$
を満たす x が存在し、法 m_1*m_2*…m_n の下で一意に定まる。
■「そうそう。この定理の主張は知ってるんだけど、具体的な解を構成する方法は知らないなと思って」
●「Garnerのアルゴリズムというのが、汎用性が高くて良いみたいですよ」
拡張Euclid互除法
■「Garnerのアルゴリズムっていうのはどういうものなの？」
●「Garnerは解の形を決めて再帰的に更新するアルゴリズムですが、まず拡張Euclid互除法をやってみましょう」
■「拡張互除法もちゃんとは知らないんだよね。確か、こんな問題を解くアルゴリズムだよね」
$$ax+by={\rm gcd}(a,b)$$
を満たす x, y を求めよ
●「はい。とりあえず a > b として、a を b で割った商を q 、余りを r とすると
$$a=bq+r$$
が成り立つので、問題が
$$(bq+r)x+by={\rm gcd}(a,b)$$
になって、bとrでまとめると
$$b(y+qx)+rx={\rm gcd}(b,r)$$
に変わります。Euclidの互除法でも利用する性質として、gcd(a,b)=gcd(b,r) というものがありますね」
■「なるほど、つまり、bx’+ry’=gcd(b,r) の解 (x’, y’) が得られたら、この対応から
$$x’=y+qx,\ y’=x$$
となるから、
$$x=y’,\ y=x’-qy’$$
という元の問題の解が得られるんだね」
●「はい。これを再帰で実装するわけですが、Euclidの互除法と同じく、最終的に gcd(a,b)x+0y=gcd(a,b) という問題を解くことになります。こうなると、x=1, y=0 としておしまいです。計算量は O(log a) ですね」
■「実装すると、こんな感じかな」
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
// ax+by=gcd(a,b) の解 (x,y) を計算 // ついでに gcd(a,b) を返す long long ext_gcd(long long a, long long b, long long &x, long long &y) { if (b == 0) { x = 1, y = 0; return a; } long long x2, y2; long long d = ext_gcd(b, a%b, x2, y2); x = y2; y = x2 - (a / b)*y2; return d; } |
●「この拡張互除法を使うと、モジュラ逆数が計算できるようになるんですよ」
法が合成数の場合のモジュラ逆数
■「モジュラ逆数って、ある数を法として見たときに、掛けて1になる数のことだよね」
$$a*x\equiv 1\ \ ({\rm mod}\ m)$$
を満たす x を、a の法 m におけるモジュラ逆数という。
●「法 m が素数なら、フェルマーの小定理を使って a^(m-2)%m を使うんですが、m が素数じゃない場合にはオイラーの定理を使うか、拡張互除法を使うことになります。オイラーの定理は面倒なので拡張互除法で解きましょう」
■「m が素数じゃない場合には、a のモジュラ逆元が存在するとは限らないよね」
●「そうですね、a と m が互いに素の場合にのみ存在します」
■「拡張互除法をどう使うの？」
●「こんな問題を解きます。
$$ax+my=1$$
この問題の解 x は、ax=1-my を満たすので、mod m で見ると、ax≡1 を満たします。」
■「なるほど、これだと、gcd(a,m)=1 のときには、拡張互除法の性質から a の逆元が1つ見つかるんだね」
|
1 2 3 4 5 6 7 8 9 10 11 12 |
// m を法としたときの a のモジュラ逆数を計算 long long mod_inv(long long a, long long m) { long long x, y; long long d = ext_gcd(a, m, x, y); if (d != 1) { // モジュラ逆数が存在しない return -1; }<br /> x %= m; if (x < 0)x += m; // 負の数の剰余を考慮 return x; } |
●「Garnerのアルゴリズムでは、合成数を法とした逆元を求める場合があるので、これを使うことになります」
Garnerのアルゴリズム
■「やっと本題だね」
●「Garnerのアルゴリズムでは、以下の条件を満たす x を計算します」
$$x\equiv r_1\ \ ({\rm mod}\ m_1)$$
$$x\equiv r_2\ \ ({\rm mod}\ m_2)$$
$$\vdots$$
$$x\equiv r_n\ \ ({\rm mod}\ m_n)$$
ただし、任意の i ≠ j について
$${\rm gcd}(m_i,m_j)=1$$
●「まず、1つ目の条件を満たす x を暫定の x とします」
■「1つ目は x≡r_1 (mod m_1) だから、x は r_1 を m_1 で割った余りにすればいいのかな」
●「そうですね。この x が、1つ目の条件を満たす、正で最小の整数になります。この x に m_1 の倍数をいくら足しても、1つ目の条件は満たしたままですよね？」
■「うん、m_1 の倍数を足しても、m_1 で割った余りは変わらないね」
●「次に、この x に適切な値を足して、1つ目の条件を満たしたまま、2つ目の条件を満たさせます」
■「1つ目の条件を満たしたまま、ということは、m_1 の倍数を足すことで2つ目の条件を達成したいということだね」
●「はい。足す数を t_1*m_1 として、
$$x+t_1m_1\equiv r_2\ \ ({\rm mod}\ m_2)$$
となるようにしたいです。式変形すると、
$$t_1m_1\equiv r_2-x\ \ ({\rm mod}\ m_2)$$
となるので、
$$t_1\equiv (r_2-x)m_{1,2}^{-1}$$
となります。ここで、m_1,2^-1 というのは、m_2 を法とした場合の m_1 のモジュラ逆数を表します。」
■「このモジュラ逆数を計算するのに、拡張互除法を使うんだね。これで得られる t_1 を使うと、新しい x は
$$x\leftarrow x+(r_2-x)m_{1,2}^{-1}m_1$$
この新しい x は、m_1 で割った余りは元の x と同じで、m_2 を法として r_2 と合同になってるね」
●「では、次は 1,2個目の条件はそのままにして、3つ目の条件 x≡r_3 (mod m_3) を満たすように変更しましょう！」
■「えっと、次は m_1 で割った余りも m_2 で割った余りもそのままにしたいから、m_1*m_2 の倍数を足さないとダメなのか。ということは、
$$x+t_2m_1m_2\equiv r_3\ \ ({\rm mod}\ m_3)$$
だから、
$$t_2\equiv (r_3-x)m_{1,3}^{-1}m_{2,3}^{-1}\ \ ({\rm mod}\ m_3)$$
となって、新しい x は
$$x\leftarrow x+(r_3-x)m_{1,3}^{-1}m_{2,3}^{-1}m_1m_2$$
とすればいいのか」
●「実際に計算するときは、m_1*m_2 のモジュラ逆数を計算するようにすれば、逆数計算が1回で済みますね」
■「実装するときは、m_1*m_2*… の値を持つようにすればいいんだね。こんな感じかな」
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
// x % m[i] = r[i] % m[i] を満たす正で最小の x を返す // i!=j に対して gcd(m[i], m[j])=1 であると仮定 // とりあえず解の存在判定は保留 long long garner(vector<long long> r, vector<long long> m) { int n = r.size(); long long m_prod = 1; // m_prod には m[i] の積を入れていく long long x = r[0] % m[0]; // 最初の条件を満たすような x の初期値 for (int i = 1; i < n; i++) { m_prod *= m[i - 1]; // m の累積積を更新 long long t = ((r[i] - x) * mod_inv(m_prod, m[i])) % m[i]; if (t < 0)t += m[i]; // 負の数の剰余の対策 x += t * m_prod; // x を更新 } return x; } |
●「これを使うと、DISCO presents ディスカバリーチャンネルコードコンテスト2019予選のD問題が解けるようになりますよ」
チップ・ストーリー ～黄金編～
■「えっと、この問題は、2～30進数で表現したときの桁和が与えられるから、それを満たす10^12以下の数を返す問題だね」
●「この問題では、x を b 進数で表現した場合の桁和は x と b-1 を法として合同である、という事実を使います」
■「10進数表現での桁和を9で割った余りが、その数を9で割った余りと一致するっていう、9の倍数判定で使う事実だね。これが10以外の基数でも成り立つんだね」
●「証明は b≡1 (mod b-1) から b^n≡1 (mod b-1) を言えば大丈夫です。これを利用すると、今回は 2～29で割った余りが与えられることになります」
■「この余りの条件を満たす数を探すのは確かに中国剰余定理っぽいけど、法が互いに素じゃないよね？」
●「法が互いに素ではないので、互いに素なものだけを使用します。例えば、27で割った余りが分かっているとき、3で割った余りは必要無いですよね？」
■「そうか、27が3の倍数だから、27を割った余りを3で割れば、3で割った余りが得られるのか。x=27p+q, q=3r+s のとき、x=3(9p+r)+s になるんだね」
●「はい。なので、29以下の素数について、その素数の累乗になっているもので最も大きいものを残せば大丈夫です。2なら16, 3なら27を残して、2,4,8 や 3,9 は無視します。このとき、6や18も無視できます」
■「ということは、使う数は 7,11,13,16,17,19,23,25,27,29 だけか」
●「こうすると、どの2数も共通の素因数を持ちませんよね」
■「そもそも1種類しか素因数を持たないんだもんね。その約数は除外してあるし」
●「この10個の数を掛けると、なんと 2329089562800 なんですよ」
■「なんとって？」
●「2329089562800 は、10^12 より大きい数です。中国剰余定理より、この余り条件を満たす数は 10^12 までに1つあるか無いかというところなんですよ」
■「なるほど、じゃあ他の数を考慮する必要が無いのか。それで 10^12 以下なんだね」
●「そういうことです。後はGarnerで計算して、出てきた解がちゃんと桁和の条件をすべて満たすかをチェックすればOKです」
■「実装はこれで大丈夫かな。C++の剰余計算 a % m は、a が負の場合に、-m+1～0 の値になる場合があるから、+m しないとダメだね」
●「C++の仕様の厄介なところですねー」
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 |
// a を b 進数表現したときの桁和 long long digit_sum(long long a, long long b) { long long ans = 0; while (a > 0) { ans += a % b; a /= b; } return ans; } int main(void) { vector<long long> a(29), r, m; set<int> use = { 7,11,13,16,17,19,23,25,27,29 }; for(int i = 0; i < 29; i++) { cin >> a[i]; if (use.find(i + 1) != use.end()) { r.push_back(a[i] % (i + 1)); // 桁和から余りを計算 m.push_back(i + 1); } } long long ans = garner(r, m); // Garnerで解の候補を計算 // 範囲チェック if (ans > 1000000000000) { cout << "invalid" << endl; return 0; } // 桁和条件のチェック for(int i = 2; i <= 30; i++) { if (digit_sum(ans, i) != a[i - 2]) { cout << "invalid" << endl; return 0; } } cout << ans << endl; } |
●「今回の問題では解をそのまま使いましたけど、問題によっては解のmodを取る場合があるので要注意ですね」
■「解を10^9+7で割った余りを出力せよ、ってやつだね」
●「その場合は、法の累積積 m_prod について、10^9+7で割った余りと、各 m[i] で割った余りと両方必要になります。このために、累積積を m[i] で割った余りをループごとに計算しないといけなくなって、O(N^2) になってしまいますね」
■「それが無ければ O(N*log(m[i])) だね。拡張互除法のモジュラ逆数計算で O(log m[i]) かかるからね」
●「そうですね。yukicoderのNo.187 中華風 (Hard)だとO(N^2)になっちゃいます」
■「複雑なことをするのかと思ったけど、意外とシンプルな考え方だね」
●「どうやら多倍長整数の実装やNTTにも使えるらしいので、そっちの応用も勉強してみたいですねー」