けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

有理数 mod 998244353 の徹底解説！

テーマ解説記事 AtCoder 有理数mod998244353 ACL

AtCoder を攻略していく上で、 998244353 に慣れることは大きな課題です。特に、有理数 \frac{p}{q} \pmod{998244353} は、初見では何を言っているか分からないと感じる方も多いかもしれないですね。下図は ABC 323 E - Playlist の補足を表示したものです1。

本記事では、有理数 \pmod{998244353} を徹底解説します！　なお、この記事は競プロアドベントカレンダー 5 日目の記事として執筆しています。

adventar.org

　

1：前提知識

まず、そもそも「 998244353 で割った余りを求めよ」という形式そのものに不慣れで、意味がわからないと感じる方は、先に次の記事を読んでもらえたらと思います。

qiita.com

最低限必要な知識を超特急でおさらいすると......

  • a + b を 998244353 で割った余りを求めたいときは、先に a, b を 998244353 で割った余りに置き換えて、それらを足して、改めて 998244353 で割った余りを求めてもよい
  • a \times b を 998244353 で割った余りを求めたいときは、先に a, b を 998244353 で割った余りに置き換えて、それらをかけて、改めて 998244353 で割った余りを求めてもよい

ということを確認しておきましょう。このように、「計算したい各値について、先に 998244353 で割った余りを求めてから、計算しても結果は変わらない」という事実はとても重要です。

また、本記事では合同式の知識を前提としています。合同式について知りたい方は、ぜひ次の記事を読んでみてください。

manabitimes.jp

　

2： \frac{y}{x} \pmod{998244353} とは何か

それでは、いよいよ \frac{y}{x} \pmod{998244353} とは何かを紐解いていきましょう！！

2-1：そもそも分数とは何か

そもそも分数って何でしょうか。分数は「割り算」です。たとえば、mod のことは一旦忘れて、

\displaystyle \frac{5}{3} = 5 \div 3

が成り立ちます。それでは「割り算」とはなんでしょうか。割り算は「かけ算の逆」です。 5 \div 3 とは、 3 \times □ = 5 となるような □ のことですね。mod の世界でも一緒です。まとめると、

\displaystyle \frac{5}{3} \pmod{998244353} とは、

3x \equiv 5 \pmod{998244353}

を満たすような x のこと

です！！！(ここでは □ を文字 x で表しました)

有理数 \pmod{998244353} とは、要するに「一次方程式の解」なのですね。それでは、一次方程式はどうやって解いたらよいでしょうか。

再び mod の世界ではなく、普通の方程式をどのように解いたかを思い出しましょう。たとえば、一次方程式 3x = 5 を解くためには、 3 の逆数である \frac{1}{3} を両辺にかけることで解けます。他にも、行列の方程式 AX = B を解くためには、 A の逆行列である A^{-1} を両辺にかけることで解けます。

　

2-2：逆元

mod の世界でも同様です。mod の世界でも、普通の数の世界における「逆数」に対応する概念として、「逆元」という概念があります。 \mod 998244353 における a の逆元とは、「 a にかけて 1 になる数」のことです。数式で書けば、

ab \equiv 1 \pmod{998244353}

を満たす b のことです。逆元は a^{-1} と表記することが多いです。つまり、

aa^{-1} \equiv 1 \pmod{998244353}

となります。いくつか具体例を見てみましょう。

\mod 998244353 において、

  • 1^{-1} \equiv 1
  • 2^{-1} \equiv 499122177
  • 3^{-1} \equiv 332748118
  • 4^{-1} \equiv 748683265
  • 5^{-1} \equiv 598946612
  • 6^{-1} \equiv 166374059
  • 7^{-1} \equiv 855638017
  • 8^{-1} \equiv 873463809
  • 9^{-1} \equiv 443664157
  • 10^{-1} \equiv 299473306

となります。1 個目の「 1 の逆元」は明らかですね。 1 \times 1 \equiv 1 ですので、 1^{-1} \equiv 1 です。2 個目の 2^{-1} \equiv 499122177 を確かめてみましょう。 2 と 499122177 をかけると、

2 \times 499122177 = 998244354 \equiv 1 \pmod{998244353}

となります。確かに 2 の逆元は 2^{-1} \equiv 499122177 であることが見て取れます。同様に、 3, 4, \dots, 10 の逆元についても、ぜひ確かめてみてください！

なお、実際のプログラムで逆元を求める方法を知りたい方は、次の記事を読んでみてください。

qiita.com

そして、AtCoder では、公式ライブラリ AtCoder Library (通称 ACL) が提供されていて、その中に modint と呼ばれる構造体があります。modint は、mod を指定して (たとえば 998244353)、四則演算や逆元計算などを直感的に操作できるようにしたものです。たとえば、 \mod 998244353 における、1 から 10 までの数の逆元は、次のように関数 inv() を用いて簡単に計算できます！！

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;

// mod 998244353 での演算を行う構造体
using mint = atcoder::modint998244353;

int main() {
    for (int a = 1; a <= 10; ++a) {
        // a の逆元を求める
        mint ra = mint(a).inv();
        
        // a の逆元を出力する
        cout << a << "^-1 = " << ra.val() << endl;
    }
}

実行結果は次のようになります。

1^-1 = 1
2^-1 = 499122177
3^-1 = 332748118
4^-1 = 748683265
5^-1 = 598946612
6^-1 = 166374059
7^-1 = 855638017
8^-1 = 873463809
9^-1 = 443664157
10^-1 = 299473306

　

2-3：分数 \pmod{98244353} を解き明かす！

それではいよいよ、 \frac{y}{x} \pmod{998244353} を実際に求める作業に入ります。例として、 \frac{5}{3} を求めてみましょう。これは、先ほどの話を踏まえると、

一次方程式 3x \equiv 5 \pmod{998244353} の解

なのでした。そして、一次方程式を解くためには、両辺に逆元をかければよいのです。今回は 3^{-1} (\equiv 332748118) を両辺にかけます。そうしますと、 3 \times 3^{-1} \equiv 1 なのですから、

x \equiv 5 \times 3^{-1}
\equiv 5 \times 332748118
\equiv 665496237 \pmod{998244353}

と求められます！　こうして、私たちは、 \frac{b}{a} のような分数 (確率や期待値など) を \mod 998244353 の世界で、整数値として表現する方法を得たのでした！！！

なお、ここまで分数を「一次方程式」を経由して解説したのですが、もっと簡単に

\displaystyle \frac{b}{a} \equiv b \times a^{-1} \pmod{998244353}

と表すこともできます。この計算は、普通の数の計算と一緒ですので、とても直感的で分かりやすいと思います。例の ACL の modint においても、割り算の演算子「/=」は次のように実装されています2。まさに、上の計算式そのものを実装していることが分かります！

mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }

つまり、まとめると「割り算とは逆元をかけることである」と一言で言えるわけですね。ACL を使って \frac{5}{3} \pmod{998244353} の値を求めてみましょう。次のコードのように書けます。

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;

// mod 998244353 での演算を行う構造体
using mint = atcoder::modint998244353;

int main() {
    mint b = 5;
    mint a = 3;
    mint answer = b / a;
    
    cout << answer.val() << endl;
}

実行結果は次のようになります。

665496237

不安な方は、同じように b * a.inv() の値も出力してみてください。きっと同じ値が返ってくるはずです！

　

補足 (1)：逆元の存在と一意性

逆元について、こんなことが疑問に浮かんだ方もいるかもしれません。

  • \mod 998244353 における a の逆元は常に存在するのだろうか？
  • \mod 998244353 における a の逆元は複数個存在することもあり得るのだろうか？

しかし、実は a が 998244353 の倍数でなければ、 a の逆元はただ 1 つだけ存在するのです。ちゃんと書くと、次の定理が成り立ちます！

p を素数として、 a を p の倍数ではない整数とする。このとき、

ab \equiv 1 \pmod{p}

を満たす整数 b ( 0 \le b \lt p) がただ 1 つ存在する。

とても美しい定理ですね！！！！　 998244353 は素数ですから、 1, 2, \dots, 998244352 の逆元はそれぞれ 1 個ずつ存在するのです。

逆元の性質の証明

証明

まず、 a, 2a, \dots, (p-1)a を p で割った余りが重複しないことを証明します。そのために、

ax \equiv ay \pmod{p}

が成り立つような x, y ( 1 \le x \lt y \le p-1) が存在すると仮定して矛盾を導きます。仮定の式を変形すると、

a(y - x) \equiv 0 \pmod{p}

となります。つまり、 a(y - x) は p の倍数となります。ここで、 a は p の倍数ではないため、 y - x は p の倍数であると言えます。しかし、 1 \le x \lt y \le p-1 ですから、 0 \lt y - x \lt p であり、 y - x は p の倍数とはなり得ません。よって、 ax \equiv ay \pmod{p} を満たす x, y ( 1 \le x \lt y \le p-1) が存在するという仮定が誤りであることが分かりました。

以上より、 a, 2a, \dots, (p-1)a を p で割った余りは互いに相異なります。特に、

ab \equiv 1 \pmod{p}

を満たす整数 b ( 0 \le b \lt p) は存在したとしてもただ 1 つしかないことが言えました。さらに、もし仮に、

ab \equiv 1 \pmod{p}

を満たす b ( 0 \le b \lt p) が存在しないと仮定すると、 a, 2a, \dots, (p-1)a を p で割った余り ( p-1 個あります) は 2, 3, \dots, p-1 のいずれか ( p-2 種類です) となります。鳩の巣原理より、 a, 2a, \dots, (p-1)a を p で割った余りの中に重複があることになります。これは矛盾です。

以上より、

ab \equiv 1 \pmod{p}

を満たす整数 b ( 0 \le b \lt p) はただ一つ存在することが証明されました。

(証明終わり)

　

補足 (2)：Fermat の小定理

せっかく美しい定理が示せたので、さらに一歩補足してみます！

Fermat の小定理について
上の証明を突き詰めると、さらに次のことも言えます。

p を素数として、 a を p の倍数ではない整数とする。このとき、 a, 2a, \dots, (p-1)a を p で割った余りを並び替えると

1, 2, \dots, p-1

となる。

実際、 p = 7 として、1 から 6 までの数同士をかけざんして九九もどきを作ると下図のようになります。

2 の列に着目すると、 2, 4, 6, 1, 3, 5 となっています。これは、 a = 2 として a, 2a, 3a, 4a, 5a, 6a を 7 で割った余りを並び替えると 1, 2, 3, 4, 5, 6 となることを表しています。

この性質を利用すると、かの有名な定理が証明できてしまいます！！　 a, 2a, \dots, (p-1)a を p で割った余りを並び替えると 1, 2, \dots, p-1 になるということは、

a \times 2a \times \dots \times (p-1)a \equiv 1 \times 2 \times \dots \times (p-1) \pmod{p}

が成り立つことを意味します。整理すると

a^{p-1} (p-1)! \equiv (p-1)! \pmod{p}

となります。 (p-1)! と p は互いに素なので、両辺を (p-1)! で割ってよくて

a^{p-1} \equiv 1 \pmod{p}

となります。これは Fermat の小定理です！！

なお、Fermat の小定理は、 a の逆元を計算する方法をも示唆しています。Fermat の小定理の式を変形すると

a \times a^{p-2} \equiv 1 \pmod{p}

となります。先ほど、まさに逆元の一意性は示したので、次のことが言えます！

a^{-1} \equiv a^{p-2} \pmod{p}

つまり、 a^{p-2} \pmod{p} を計算することで逆元が求められます。実際、ACL の modint においても、逆元を求める関数 inv() は次のように実装されています。このコードで pow(umod() -2) の部分がそれを表しています。

mint inv() const {
    if (prime) {
        assert(_v);
        return pow(umod() - 2);
    } else {
        auto eg = internal::inv_gcd(_v, m);
        assert(eg.first == 1);
        return eg.second;
    }
}

ただし、逆元を Fermat の小定理で求める方法は、mod が素数の場合にしか通用しません。合成数である場合 (else 文の中身) は、別の方法で逆元を求めています。

　

補足 (3)：なぜ mod をとるのか

そもそも根本的な疑問として、確率や期待値の値を double 型で計算するのではダメなのだろうか......とも思われるかもしれません。実際、そのような形式で出題される問題もあります。

しかし、double 型で確率や期待値を求めると、誤差の問題が常につきまといます。一方、 \mod 998244353 で求めた整数値を答えさせると、誤差の問題はまったくありません。出題者も安心して出題できますし、参加者も安心して回答できますね。

　

3：分数の足し算、引き算、かけ算、割り算

さて、次のステップに行きましょう！！！

3-1：分数の足し算

早速ですが、次の問題を解いてみてください。いずれ、ABC-D などで出題されるかもしれません。

例題 (1)

正の整数 a, b, c, d が与えられるので、

\displaystyle \frac{b}{a} + \frac{d}{c}

を計算してください。ただし、計算結果を \displaystyle \frac{Q}{P} と表したとき、その値を \mod 998244353 で出力してください。

制約

  • 1 \le a, b, c, d \lt 998244353

解法 1

この問題を読んで、次のように考えた方もいるかもしれないですね。通分すると、

\displaystyle \frac{b}{a} + \frac{d}{c} = \frac{bc + ad}{ac}

となりますので、ACL の modint を用いて (b * c + a * d) / (a * c) の値を出力すれば解けると。つまり、次のように実装できます。

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;

// mod 998244353 での演算を行う構造体
using mint = atcoder::modint998244353;

int main() {
    int ia, ib, ic, id;
    cin >> ia >> ib >> ic >> id;
    mint a = ia, b = ib, c = ic, d = id;
    
    mint bunshi = b * c + a * d;
    mint bunbo = a * c;
    mint answer = bunshi / bunbo;
    cout << answer.val() << endl;
}

解法 2

もちろん、上の解法 1 は正解です！！　

......が、 \frac{b}{a} + \frac{d}{d} の式変形の計算を手元でやるのはミスのリスクがあります。実は、もっと簡単に次のコードでよいのです！

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;

// mod 998244353 での演算を行う構造体
using mint = atcoder::modint998244353;

int main() {
    int ia, ib, ic, id;
    cin >> ia >> ib >> ic >> id;
    mint a = ia, b = ib, c = ic, d = id;
    
    mint answer = b / a + d / c;
    cout << answer.val() << endl;
}

どちらのコードを用いても、たとえば a = 4, b = 5, c = 6, d = 7 とすると、答えは 915057326 となります。ぜひ確かめてみてください！

補足： \frac{b}{a} + \frac{d}{c} の計算について

上の解法では、

  • \displaystyle \frac{b}{a} \equiv p \pmod{998244353}
  • \displaystyle \frac{d}{c} \equiv q \pmod{998244353}

としたときに、

\displaystyle \frac{b}{a} + \frac{d}{c} \equiv p + q \pmod{998244353}

が成り立つということを証明なしで使いました。一応これを証明しておきましょう。一見すると、整数 a, b, x, y に対して a \equiv x かつ b \equiv y ならば a + b \equiv x + y が成り立つことから自明だと思われるかもしれません。しかし、 \frac{b}{a} や \frac{d}{c} は整数ではないので、証明が必要ですね。

証明

証明

\displaystyle \frac{b}{a} + \frac{d}{c} = \frac{bc + ad}{ac}

であるから、

\displaystyle \frac{b}{a} + \frac{d}{c} \equiv p + q \pmod{998244353}

を示すためには、

ac(p + q) \equiv bc + ad \pmod{998244353}

であることを示せばよいです。ここで、 \frac{b}{a} \equiv p, \frac{d}{c} \equiv q より、

  • ap \equiv b \pmod{998244353}
  • cq \equiv d \pmod{998244353}

が成り立つことに着目して、

ac(p + q) = (ap)c + a(cq) \equiv bc + ad \pmod{998244353}

が成り立ちます。以上より、示せました。

(証明終わり)

　

3-2：分数のかけ算

足し算に続いて、かけ算もやってみましょう！

例題 (2)

正の整数 a, b, c, d が与えられるので、

\displaystyle \frac{b}{a} \times \frac{d}{c}

を計算してください。ただし、計算結果を \displaystyle \frac{Q}{P} と表したとき、その値を \mod 998244353 で出力してください。

制約

  • 1 \le a, b, c, d \lt 998244353

解説

足し算の場合と同様に、次のコードのように、式をそのまま実装すればよいと期待できます。実際、これで正解です！！

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;

// mod 998244353 での演算を行う構造体
using mint = atcoder::modint998244353;

int main() {
    int ia, ib, ic, id;
    cin >> ia >> ib >> ic >> id;
    mint a = ia, b = ib, c = ic, d = id;
    
    mint answer = (b / a) * (d / c);
    cout << answer.val() << endl;
}

補足： \frac{b}{a} \times \frac{d}{c} の計算について

上の解法では、足し算の場合と同様に、

  • \displaystyle \frac{b}{a} \equiv p \pmod{998244353}
  • \displaystyle \frac{d}{c} \equiv q \pmod{998244353}

としたときに、

\displaystyle \frac{b}{a} \times \frac{d}{c} \equiv pq \pmod{998244353}

が成り立つことを使っています。一応、これを証明しておきましょう。

証明

証明

\displaystyle \frac{b}{a} \times \frac{d}{c} = \frac{bd}{ac}

であるから、

\displaystyle \frac{b}{a} \times \frac{d}{c} \equiv pq \pmod{998244353}

を示すためには、

ac(pq) \equiv bd \pmod{998244353}

であることを示せばよいです。ここで、 \frac{b}{a} \equiv p, \frac{d}{c} \equiv q より、

  • ap \equiv b \pmod{998244353}
  • cq \equiv d \pmod{998244353}

が成り立つことに着目して、

ac(pq) = (ap)(cq) \equiv bd \pmod{998244353}

が成り立ちます。以上より、示せました。

(証明終わり)

　

3-3：分数の四則演算のまとめ

ここまで、分数の足し算と引き算について見てきました。同様のことは引き算や割り算でも成り立ちます。まとめると、次のことが言えます。

a, b, c, d を 998244353 の倍数ではない整数として、

  • \displaystyle \frac{b}{a} \equiv p \pmod{998244353}
  • \displaystyle \frac{d}{c} \equiv q \pmod{998244353}

とする。このとき、

  • \displaystyle \frac{b}{a} + \frac{d}{c} \equiv p + q \pmod{998244353}
  • \displaystyle \frac{b}{a} - \frac{d}{c} \equiv p - q \pmod{998244353}
  • \displaystyle \frac{b}{a} \times \frac{d}{c} \equiv pq \pmod{998244353}
  • \displaystyle \frac{b}{a} \div \frac{d}{c} \equiv \frac{p}{q} \pmod{998244353}

が成り立つ。

こうしてみると、 \mod 998244353 における分数計算は、本当に直感的にできることが分かりますね！

　

4：実際の確率計算の例

ここまでに得た知見をもとに、実際に簡単な問題を解いてみましょう。

例題 (3)

1 から N までの目が等確率で出るサイコロがあります。このサイコロを 3 回振るとき、目の最大値が A である確率を求めてください。

ただし、確率を \displaystyle \frac{Q}{P} と表したとき、その値を \mod 998244353 で出力してください。

制約

  • 2 \le N \le 100
  • 1 \le A \le N

解法

目の最大値が A である確率は、次のように求められます。

(3 個のサイコロの目がすべて A 以下である確率)
- (3 個のサイコロの目がすべて A-1 以下である確率)

最大値が A であるということは、すべての目が A 以下であることが必要です。しかし、 A の目が存在しなければなりません。そこで、 A の目が存在しない場合......つまり、すべての目が A-1 以下であるような場合の確率を引くのです。そして、

  • 3 個のサイコロの目がすべて A 以下である確率： \displaystyle (\frac{A}{N})^{3}
  • 3 個のサイコロの目がすべて A-1 以下である確率： \displaystyle (\frac{A-1}{N})^{3}

であるから、求める確率は次のようになります。

\displaystyle (\frac{A}{N})^{3} - (\frac{A-1}{N})^{3}

あとは、これを実装すればよいです。modint を使えば、次のように、上の式をそのまま直感的に実装できます！

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;

// mod 998244353 での演算を行う構造体
using mint = atcoder::modint998244353;

int main() {
    int N, A;
    cin >> N >> A;
    
    mint p = mint(A) / N;
    mint q = mint(A-1) / N;
    mint res = p * p * p - q * q * q;
    cout << res.val() << endl;
}

　

5：おわりに

ここまで、有理数 \pmod{998244353} について解説してきました。確率や期待値のように「分数で表したいもの」を求めるときに、時々そのような形式で出力することを要求されることがあります。

しかし modint さえあれば、結局まったく難しいことはなく、数式をそのまま直感的に実装すればよいのです！

  1. このような問題文中の注釈は、分からない人が分かるようにするためというよりは、「コンテストの問題文として成立するようにするため」につけているものだと思われます。有理数 \pmod{998244353} という概念はそれほど一般的なものではないため、このような注釈が必要となっているのだと思われます。↩
  2. ACL は実装も公開されていますので、中身を読むと勉強になります。↩