962

記事の改善

記事の情報

編集履歴編集リクエスト一覧 LGTMしたユーザー一覧 Markdown で本文を見る

問題がある記事を報告する

この記事は最終更新日から1年以上が経過しています。

@drken (株式会社NTTデータ数理システム)

投稿日 2018年08月01日

更新日 2021年01月31日

Organization

「1000000007 で割ったあまり」の求め方を総特集！〜逆元から離散対数まで〜

アルゴリズム数学競技プログラミング整数論数学やり直し

プログラミングコンテストの問題に挑んでいると、

「ある量」を求めたい。
ただし答えがとても大きくなる可能性があるので、その量を $1000000007$ で割ったあまりを求めよ。

というタイプの出題をよく見かけます。今回はこのタイプの問題に対する対策をまとめます！！！

下表に示すように、四則演算・累乗・二項係数といった話題を集大成します。なお、表中でライブラリ化の必要性について書いていますが、そもそも「1000000007 で割ったあまり」を直感的に扱える構造体 modint を用意すると便利です。modint についても、最後の方で簡単に触れます。

演算	やり方	ライブラリ化の必要性
足し算 $a + b$	計算途中でもあまりとってよい
掛け算 $a \times b$	計算途中でもあまりとってよい
引き算 $a - b$	計算途中でもあまりとってよい、ただし最終結果が負になるなら $1000000007$ を足す
割り算 $a \div b$	$b$ の逆元をかける	「逆元」があると便利
累乗 $a^{n}$	二分累乗法	あると便利
二項係数 $_{n} C_{r}$	$a!$ , $(a!)^{- 1}$ のテーブルを作る	あると便利
離散対数 $\log_{a} b$	Baby-Step Giant-Step 法など	あると時々使える
平方剰余 $\sqrt{a}$	Tonelli-Shanks のアルゴリズムなど	あると時々使える

1. なぜ 1000000007 で割るのか？

最初はこのような設問を見るとぎょっとしてしまいますが、実はとても自然な問題設定です。 $1000000007$ で割らないと、答えの桁数がとてつもなく大きくなってしまうことがあります。このとき以下のような問題が生じます:

多倍長整数がサポートされている言語とされていない言語とで有利不利が生じる
10000 桁にも及ぶような巨大な整数を扱うとなると計算時間が膨大にかかってしまう

1 番目の事情はプログラミングコンテストに特有のものと思えなくもないですが、2 番目の事情は切実です。整数の足し算や掛け算などを実施するとき、桁数があまりにも大きくなると桁数に応じた計算時間がかかってしまいます。実用的にもそのような巨大な整数を扱うときは、いくつかの素数で割ったあまりを計算しておいて、最後に中国剰余定理を適用して復元することも多いです。

なぜ 1000000007 なのか？

これについては、それほど深い理由はなさそうですが

$10^{9}$ より大きい最小の素数で、手頃な素数である (素数であることは「割り算」などで重要)
$1000000007$ 未満の数同士を足しても 32 ビット整数におさまる
$1000000007$ 未満の数同士を掛けても 64 ビット整数におさまる

というのはあります。 $1000000007$ 以外にも $1000000009$ や $998244353$ で割ったあまりを求めさせるケースも多いです。

2. 足し算、引き算、掛け算

まずは簡単な足し算、引き算、掛け算について見て行きます。基本的なノウハウは「計算の途中過程で積極的にあまりをとってよい」ということになります。たとえば

$a = 111111111$
$b = 123456789$
$c = 987654321$

のときに $a \times b \times c$ を $1000000007$ で割ったあまりを求めたいとします。このとき $M O D = 1000000007$ として、

$(a \times b \times c)$ % $M O D$ を計算する
$((a \times b)$ % $M O D \times c)$ % $M O D$ を計算する

の両者は原理的には同じ答えになります。後者は計算の途中過程であまりをとっています。Python で実験してみます。

MOD = 1000000007
a = 111111111
b = 123456789
c = 987654321
print a * b * c % MOD
print a * b % MOD * c % MOD

結果は以下のようになります:

769682799
769682799

オーバーフローに注意！！！

C++ で足し算・引き算・掛け算を行うときには「計算の途中過程であまりをとってよい」という控えめな言い方ではなく、むしろ

掛け算を 1 回やる度に、毎回、あまりをとった方が無難

という感じです。C++ で上の実験を試してみるとおかしなことになります:

#include <iostream>
using namespace std;

const int MOD = 1000000007;

int main() {
  long long a, b, c;
  a = 111111111;
  b = 123456789;
  c = 987654321;
  cout << a * b * c % MOD << endl;       // a*b*c がオーバーフローするので間違った答えに
  cout << a * b % MOD * c % MOD << endl; // a*b で一旦 MOD をとる
}

結果は以下のようになりました。下の 769682799 の方が正しい答えです。

874753149
769682799

下の「a * b * c % MOD」の方の答えが違ってしまう理由は、a * b * c を計算した時点で 64 ビット整数におさまらずにオーバーフローしてしまっているからです。C++ で $1000000007$ で割ったあまりを計算するとき、とくに掛け算を扱うときは

64 ビット整数を使う (足し算のみなら 32 ビット整数でも OK)
掛け算する度に $1000000007$ で割っておく

という風にするのがよいです。
上記と似た事情は多倍長整数が扱える Python などでも言えます。 $1000000007$ で割るのはそれ自体もコストなので「掛け算する度に毎回割る」というのではかえって計算時間が伸びてしまう恐れもあるのですが、定期的に $1000000007$ で割るようにしないと桁数が大きくなって計算時間が増大してしまいます。

2-2. 引き算

引き算も基本的には一緒で「計算の途中過程で積極的にあまりをとってよい」という感じなのですが、「負の数のあまり」について少しだけ注意が必要です。

$- 17$ を $5$ で割ったあまり

は正しくは 3 になります。そもそも「 $a$ を $b$ で割ったあまり」というのは「 $a$ に $b$ の倍数を足し引きして得られる数のうち、 $0$ 以上 $b$ 未満のもの」を表しています。 $- 17$ に $5$ の倍数である $20$ を足すと $3$ になります。しかし C++ ではあまりの定義が少しだけ変なことになります:

#include <iostream>
using namespace std;

int main() {
  cout << (-17 % 5) << endl; // 3 になってほしいが -2 になる
}

これを実行すると $- 2$ が出力されるはずです。それは

$17$ を $5$ で割ったあまりは $2$
$- 17$ を $5$ で割ったあまりはマイナスをつけて $- 2$

というロジックです。完全に間違っているわけではなくて、 $- 2$ にさらに $5$ を足すと $3$ になって正しい答えになります。したがって「あまりを求めた結果が負になったら法 $1000000007$ を足す」という風にすれば良いです。

このようなケースが発生する例として

$a = 2000000020$
$b = 20$

として $a - b$ を $1000000007$ で割ったあまりを計算したいが、 $a$ を先に $1000000007$ で割っておこうとする場合などがあります。

#include <iostream>
using namespace std;

const int MOD = 1000000007;

// 負の数にも対応した % 演算
long long mod(long long val, long long m) {
  long long res = val % m;
  if (res < 0) res += m;
  return res;
}

int main() {
  int a = 2000000020;
  int b = 20;

  cout << "普通に計算して余りを求める: " << (a - b) % MOD << endl;
  cout << "余り求めてから計算して余りを求める: " << ((a%MOD) - (b%MOD)) % MOD << endl;
  cout << "余り求めてから計算して余りを求める (対策済): " << mod((a%MOD) - (b%MOD), MOD) << endl;
}

結果は下のようになります:

普通に計算して余りを求める: 999999993
余り求めてから計算して余りを求める: -14
余り求めてから計算して余りを求める (対策済): 999999993

なお、Python を用いる場合はずっと楽で、負の数に対してもきちんとあまりを計算することができます:

>>> -17 % 5
3

3. 割り算 a ÷ b

3-1. mod p の世界における割り算とは

掛け算では「掛け算する度に $1000000007$ で割っておく」としてよかったです。しかし割り算では少し頭を悩ませることになります。例えば、

$12345678900000 \div 100000 = 123456789$

という計算において、一度 $12345678900000$ を $1000000007$ で割ったあまり ( $678813585$ になります) を出してから計算しようとすると、

$678813585 \div 100000$

を計算することになります。しかしこれでは割り切れなくなってしまいます。割り算については「あまりをとってから割り算する」という方法が単純には通用しないことがわかります。割り算が登場する場面では、泣く泣く多倍長整数をそのままの形で扱うしかないのでしょうか？

実はそんなことはないのです。一般に $p$ を素数として $\mod p$ の世界における「割り算」について掘り下げて考えてみましょう。具体例として $\mod 13$ について考えてみます。

$9 \div 4$ はなにか？

と問われたら、普通の整数の世界では割り切れないですが、 $\mod 13$ の世界では割り切れます:

$9 \div 4 \equiv 12 (\mod 13)$

です。なぜなら、 $4 \times 12 = 48 \equiv 9 (\mod 13)$ が成り立つからです。同様にして

$1 \div 4 \equiv 10 (\mod 13)$
$2 \div 4 \equiv 7 (\mod 13)$
$3 \div 4 \equiv 4 (\mod 13)$
$4 \div 4 \equiv 1 (\mod 13)$ 　(通常の割り算と一緒です)
$5 \div 4 \equiv 11 (\mod 13)$
$6 \div 4 \equiv 8 (\mod 13)$
$7 \div 4 \equiv 5 (\mod 13)$
$8 \div 4 \equiv 2 (\mod 13)$ 　(通常の割り算と一緒です)
$9 \div 4 \equiv 12 (\mod 13)$
$10 \div 4 \equiv 9 (\mod 13)$
$11 \div 4 \equiv 6 (\mod 13)$
$12 \div 4 \equiv 3 (\mod 13)$ 　(通常の割り算と一緒です)

となります。一般に、 $p$ を素数として $\mod p$ の世界ではこのような割り算ができます (数学的には $mod p$ の世界が「体」であるということになります)。これは以下の定理から来ています。この定理の証明は

フェルマーの小定理の証明と使い方

に書いたので、参考にしていただけたらと思います。

$p$ を素数とし、 $b$ を $p$ で割り切れない整数とする。このとき $a$ を整数として

b x \equiv a (\mod p)

を満たすような $x$ は $mod p$ において一意に存在する

以上を踏まえて元の問題に戻ると、

$678813585 \div 100000 (\mod 1000000007)$

はちゃんと一意に特定できます。答えは $123456789$ です。実際

$100000 \times 123456789 \equiv 678813585 (\mod 1000000007)$

が成立しています。

さて、ここまでは「 $mod p$ の世界で割り算ができる」という事実の紹介をしましたが、具体的な計算方法については述べていませんでした。次節以降で具体的な方法を述べます。

3-2. mod p における「逆元」

「 $a \div b (\mod p)$ 」を計算する方法を考えます。実は少し考えてみると

$1 \div b (\mod p)$

だけ計算できればいいことがわかります。なぜなら、

$a \div b \equiv a \times (1 \div b) (\mod p)$

が成り立つからです。つまり $1 \div b$ さえ求めることができれば、これを $a$ 倍すればよいです。このような $1 \div b$ を「 $mod p$ における $b$ の逆元」と呼びます。これは

$b$ をかけると $1$ になる数 ( $mod p$ の意味で)

でもあります。例えば $mod 13$ において $2$ の逆元は $7$ です。 $b$ の逆元はしばしば $b^{- 1}$ という風に表記します。まとめると、

$a \div b \equiv a \times b^{- 1} (\mod p)$

です。こう書いてみると、普通の実数における計算 ( $3 \div 7 = \frac{3}{7}$ など) とピッタリ対応していることがわかります。逆元について具体的な数値を示すと、 $mod 13$ において

$1^{- 1} \equiv 1$ 　( $1 \times 1 \equiv 1 (\mod 13)$ のため)
$2^{- 1} \equiv 7$ 　( $2 \times 7 = 14 \equiv 1 (\mod 13)$ のため)
$3^{- 1} \equiv 9$ 　( $3 \times 9 = 27 \equiv 1 (\mod 13)$ のため)
$4^{- 1} \equiv 10$ 　( $4 \times 10 = 40 \equiv 1 (\mod 13)$ のため)
$5^{- 1} \equiv 8$ 　( $5 \times 8 = 40 \equiv 1 (\mod 13)$ のため)
...

となります。逆元を求めるアルゴリズム、及び $mod m$ における割り算を実行する実装を先に示すと以下のようになります:

#include <iostream>
using namespace std;

// mod. m での a の逆元 a^{-1} を計算する
long long modinv(long long a, long long m) {
    long long b = m, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= m;
    if (u < 0) u += m;
    return u;
}

// a ÷ b mod. MOD を計算してみる
int main() {
    const int MOD = 1000000007;

    long long a = 12345678900000;
    long long b = 100000;

    // a を 10000000007 で割ってから b の逆元をかけて計算
    a %= MOD;
    cout << a * modinv(b, MOD) % MOD << endl;
}

上のソースコードでは実験として

$a = 12345678900000$
$b = 100000$

として、 $a \div b (\mod 1000000007)$ を計算しています。答えは $123456789$ になってほしいです。結果は以下のようになりました:

123456789

次節以降、逆元 $b^{- 1}$ を求めるアルゴリズム (上の実装における modinv の中身) を考えて行きます。

3-3. 逆元の求め方の概要

逆元の求め方にはメジャーな方法が 2 つあります。計算時間的にはどちらも $O (\log p)$ になります。

Fermat の小定理を用いる方法
拡張 Euclid の互除法を用いる方法

前者は「実装が単純明快」というメリットがありますが、「法 $p$ が素数でないと使えない」というデメリットがあります。後者は実装はやや複雑になりますが

法 $p$ が素数でなくても、後述する逆元存在条件を満たせば使える
Fermat の小定理を用いた方法よりも平均的に高速と言われている

というメリットがあります。kirika_comp さんの整数論テクニック集でも、逆元ライブラリを作るときには「拡張 Euclid の互除法」を用いる方を推しています。

3-4. Fermat の小定理による逆元計算

Fermat の小定理とは、

$p$ を素数、 $a$ を $p$ の倍数でない整数として

a^{p - 1} \equiv 1 (\mod p)

が成立する

というものです。ちょっと確かめてみると $p = 7$ として

$1^{6} = 1 \equiv 1 (\mod 7)$
$2^{6} = 64 \equiv 1 (\mod 7)$
$3^{6} = 729 \equiv 1 (\mod 7)$
$4^{6} \equiv (- 3)^{6} = 3^{6} \equiv 1 (\mod 7)$
$5^{6} \equiv (- 2)^{6} = 2^{6} \equiv 1 (\mod 7)$
$6^{6} \equiv (- 1)^{6} = 1^{6} \equiv 1 (\mod 7)$

となって、確かに成立しています。後半 3 つは「それはそう」という感じですが、なかなか不思議な定理です。Fermat の小定理の証明は今回は置いておくとして、これを利用して逆元を求めてみます。定理の式を変形すると

$a \times a^{p - 2} \equiv 1 (\mod p)$

となります。これはつまり、 $a^{p - 2}$ が $a$ の逆元になっていることを意味しています。したがって、 $a^{p - 2} (\mod p)$ の計算をすればよいのですが、愚直に「 $a$ を $p - 2$ 回掛け算する」という方法では $O (p)$ かかってしまいます。4 章で見る二分累乗法を用いることで $O (\log p)$ で求めることができます。

#include <iostream>
using namespace std;

// a^n mod を計算する
long long modpow(long long a, long long n, long long mod) {
    long long res = 1;
    while (n > 0) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

// a^{-1} mod を計算する
long long modinv(long long a, long long mod) {
    return modpow(a, mod - 2, mod);
}

int main() {
    // mod. 13 での逆元を求めてみる
    for (int i = 1; i < 13; ++i) {
        cout << i << " 's inv: " << modinv(i, 13) << endl;
    }
}

上のコードでは $\mod 13$ での逆元を一通り出してみました。結果は以下のようになりました:

1 's inv: 1
2 's inv: 7
3 's inv: 9
4 's inv: 10
5 's inv: 8
6 's inv: 11
7 's inv: 2
8 's inv: 5
9 's inv: 3
10 's inv: 4
11 's inv: 6
12 's inv: 12

例えば $4$ などを見ると逆元は $10$ で、確かに $4 \times 10 = 40 \equiv 1 (\mod 13)$ になっていることがわかります。mod. 13 スピードといったトランプ競技では、これらの値を瞬時に言えるようになることが重要になります。

3-5. 拡張 Euclid の互除法による逆元計算

$x$ が $mod p$ における $a$ の逆元であるとは、

$a x \equiv 1 (\mod p)$

が成立するということでもあります。これは $a x - 1$ が $p$ で割り切れることを意味しています。すなわち

$a x + p y = 1$ を満たす整数 $y$ が存在すること

を意味しています。そして上記を満たす $x$ はまさに拡張 Euclid の互除法によって求めることができます。拡張 Euclid の互除法とは以下のことができるアルゴリズムです:

互いに素な 2 整数 $a, b$ が与えられたときに

$a x + b y = 1$

を満たす $(x, y)$ を 1 つ求める

拡張ユークリッドの互除法の詳細や計算時間の解析など、詳しくは

などを読んでいただけたらと思います。ここでは非再帰拡張 Euclid の互除法による逆元計算の実装を示します。

#include <iostream>
using namespace std;

long long modinv(long long a, long long m) {
    long long b = m, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= m; 
    if (u < 0) u += m;
    return u;
}

int main() {
    // mod. 13 での逆元を求めてみる
    for (int i = 1; i < 13; ++i) {
        cout << i << " 's inv: " << modinv(i, 13) << endl;
    }
}

結果は上と同じく以下のようになります。高速かつ $m$ が素数でなくてもよいということで、かなりおススメのライブラリです。

1 's inv: 1
2 's inv: 7
3 's inv: 9
4 's inv: 10
5 's inv: 8
6 's inv: 11
7 's inv: 2
8 's inv: 5
9 's inv: 3
10 's inv: 4
11 's inv: 6
12 's inv: 12

3-6. 逆元が存在する条件

上記の方法の注意点を最後に述べます。
どんな定理・アルゴリズムを使用するときにも、その前提条件をきちんと把握しておくことは重要です。「逆元」がちゃんと存在する条件を整理します。大前提として、逆元が常に存在するわけではないです。それは行列 $A$ の逆行列が常に存在するとは限らない (行列式が $0$ でないときに限り存在) のと同様です。

$mod p$ での $a$ の逆元が存在する条件は、 $p$ と $a$ とが互いに素であること

です。よくある誤解として以下のものがあります:

$p$ が素数でないときは、必ず逆元は存在しない (誤解 1)
$p$ が素数のときは、必ず逆元が存在する (誤解 2)

まず前者についてですが、 $p$ と $a$ とが互いに素でありさえすればよいです。例えば $p = 12$ のときは $a \equiv 1, 5, 7, 11$ がその条件を満たします。

そして後者についてですが、 $a = 0, p, 2 p, 3 p, \dots$ のとき、すなわち $a$ が $p$ の倍数のときには逆元は存在しないです。これは通常の整数における「ゼロ割りはできない」というのにピッタリ対応しています。特に $a = p$ のときに逆元が存在しないことは見落としがちです。 $p = 1000000007$ とかだと、そんな巨大な $a$ で割り算する機会はほとんどないのですが、 $p = 10007$ くらいだと、 $_{n} C_{r} (\mod p)$ の計算をするときに罠にハマったりします。

4. 累乗 a^n

愚直に $a$ を $n$ 回掛け算すると $O (n)$ かかってしまいます。これを $O (\log n)$ にする二分累乗法と呼ばれるテクニックがよく知られています。例えば $3^{16}$ を計算したいとします。このとき

$3$ を二乗すると $3^{2}$ が求まる
それを二乗すると $3^{4}$ が求まる
それを二乗すると $3^{8}$ が求まる
それを二乗すると $3^{16}$ が求まる

という風にすれば僅か 4 回の掛け算で求めることができます。今は冪が「16」で 2 の累乗だったから簡単だったのですが、冪が 2 の累乗でなくても $O (\log n)$ 回の掛け算で求めることができます。

例えば $3^{45}$ を計算したい場合には次のようにします。まず $45$ を二進法展開すると

$45 = 2^{0} + 2^{2} + 2^{3} + 2^{5}$

になります。従って、

$3^{45} = 3^{2^{0} + 2^{2} + 2^{3} + 2^{5}} = 3^{2^{0}} \times 3^{2^{2}} \times 3^{2^{3}} \times 3^{2^{5}}$

になるので、先程と同様にして $3^{2^{0}}, 3^{2^{1}}, 3^{2^{2}}, 3^{2^{3}}, 3^{2^{4}}, 3^{2^{5}}$ を求めておけば、 $3^{45}$ も計算することができます。以上の着想に基づいて $a^{n} (\mod m)$ を求める実装を示します:

#include <iostream>
using namespace std;

// a^n mod を計算する
long long modpow(long long a, long long n, long long mod) {
    long long res = 1;
    while (n > 0) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

int main() {
    // 3^45 mod. 1000000007 を計算してみる
    cout << modpow(3, 45, 1000000007) << endl;
}

なお、こうした二分累乗法テクニックは、なにも $\mod 1000000007$ という設定だけでなく

行列累乗
ツリー上の LCA を求めるダブリング　(蟻本「4-3. グラフマスターへの道」の LCA など)
ダブリングを用いた DP　(蟻本「4-4. 厳選！頻出テクニック(2)」のダブリングなど)
kirika_comp さんの資料の「3.1 モノイド的構造を見つけて二分累乗する (Lv. 2)」

など様々な場面で活用することができます。

5. 二項係数 nCr

「 $1000000007$ で割ったあまりを求めよ」という問題において、二項係数 $_{n} C_{r} (\mod 1000000007)$ の計算が必要になる場面は非常に多いです。典型的な求め方を

よくやる二項係数 (nCk mod. p)、逆元 (a^-1 mod. p) の求め方

に書きました。求め方の詳細な原理は上記記事に譲り、以下に実装例を示します。考え方としては

$1!, 2!, 3!, \dots$
$(1!)^{- 1}, (2!)^{- 1}, (3!)^{- 1}, \dots$

を前処理で計算しておいて、

$_{n} C_{r} = \frac{n!}{r! (n - r)!} = (n!) (r!)^{- 1} ((n - r)!)^{- 1}$

を毎回計算するというものです。前処理は $O (n)$ 、毎回の二項係数計算は $O (1)$ でできます。

#include <iostream>
using namespace std;

const int MAX = 510000;
const int MOD = 1000000007;

long long fac[MAX], finv[MAX], inv[MAX];

// テーブルを作る前処理
void COMinit() {
    fac[0] = fac[1] = 1;
    finv[0] = finv[1] = 1;
    inv[1] = 1;
    for (int i = 2; i < MAX; i++){
        fac[i] = fac[i - 1] * i % MOD;
        inv[i] = MOD - inv[MOD%i] * (MOD / i) % MOD;
        finv[i] = finv[i - 1] * inv[i] % MOD;
    }
}

// 二項係数計算
long long COM(int n, int k){
    if (n < k) return 0;
    if (n < 0 || k < 0) return 0;
    return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD;
}

int main() {
    // 前処理
    COMinit();

    // 計算例
    cout << COM(100000, 50000) << endl;
}

6. 離散対数 log_a b

「足し算」「引き算」「掛け算」「割り算」「累乗」とやって来ると「対数」も考えてみたくなります。
蟻本の「発展的トピック」のところでリストアップされているテーマの 1 つでもあります。「累乗」までとは違って数え上げ問題などで必要になることはなく、純粋に整数論的問題において必要なテーマなので、今までと気色の違う話ですが一応書いてみたいと思います。また、今までは「答えを $m$ で割ったあまり」を求める話でした。今回も $\mod m$ の世界について考察する話ではありつつも、答えを $m$ で割るという話ではなくなります。

$a^{x} \equiv b (\mod m)$

となる最小の正の整数 $x$ を求める問題を考えます。そのような $x$ が存在しないときにはその旨を報告するようにします。 $m$ は素数でなくてもよいですが、簡単のために $a$ と $m$ は互いに素であるとしておきます。

一瞬手の付けようがない絶望感を感じる問題ですが、実は $a^{x} (\mod m)$ は周期性を持ちます。例えば $m$ が素数のときは Fermat の小定理から

$a^{m - 1} \equiv 1 (\mod m)$

が成立するので、 $a^{0}, a^{1}, a^{2}, \dots (\mod m)$ は、 $a^{0}, a^{1}, a^{2}, \dots, a^{m - 2} (\mod m)$ をずっと繰り返していく感じになります。 $m$ が合成数の場合も同様で、 $a$ と $m$ が互いに素なとき、 $ϕ (m)$ を Euler 関数として、

$a^{ϕ (m)} \equiv 1 (\mod m)$

が成立します (Euler の定理)。 $ϕ (m) \leq m$ ですので、最悪でも $a^{1}, a^{2}, \dots, a^{m}$ $(m o d . m)$ まで調べれば元の問題を解くことができることがわかりました。しかしこれでは $O (m)$ の計算時間がかかってしまいます。これを平方分割的な探索方法を用いて $O (\sqrt{m})$ に高速化します。まず $a^{x} \equiv b (\mod m)$ に解があるとしたら $1 \leq x \leq m$ の範囲にあることはわかっているので、

$x = p \sqrt{m} + r$

として、 $p, r$ $(0 \leq p \leq \sqrt{m}, 0 \leq r < \sqrt{m})$ を求める方針で考えてみます。このとき

$a^{x} \equiv b (\mod m)$
⇔ $a^{p \sqrt{m} + r} \equiv b (\mod m)$
⇔ $b (a^{p \sqrt{m}})^{- 1} \equiv a^{r} (\mod m)$
⇔ $b (a^{- \sqrt{m}})^{p} \equiv a^{r} (\mod m)$
⇔ $b A^{p} \equiv a^{r} (\mod m)$ ( $A \equiv a^{- \sqrt{m}}$ とおく)

と変形できます。右辺は { $a^{0}, a^{1}, a^{2}, \dots, a^{\sqrt{m} - 1}$ } という $O (\sqrt{m})$ 個の値を取りえます。これを例えば std::unordered_map といったハッシュに蓄えておきます。この作業は平均的に $O (\sqrt{m})$ の計算時間で実行できます (ハッシュへの要素挿入にかかる時間は平均的に $O (1)$ )。この前処理を終えた後は、

$b A^{0}$ がハッシュに入っているか調べる、入っていたらそれを $a^{k}$ として、 $p = 0, r = k$ が答え
$b A^{1}$ がハッシュに入っているか調べる、入っていたらそれを $a^{k}$ として、 $p = 1, r = k$ が答え
$\dots$
$b A^{\sqrt{m}}$ がハッシュに入っているか調べる、入っていたらそれを $a^{k}$ として、 $p = \sqrt{m}, r = k$ が答え
ここまで調べて入っていなかったら解なし

という風にすればよいです。 $O (\sqrt{m})$ 回のループを回していて、各ループの判定はハッシュを用いて平均的に $O (1)$ でできるので全体として平均的に $O (\sqrt{m})$ の計算量になります。まとめると

ハッシュを作る前処理: 平均的に $O (\sqrt{m})$
ハッシュを用いた探索: 平均的に $O (\sqrt{m})$

ということになり、全体を通しても平均的に $O (\sqrt{m})$ の計算時間でできます。このように、周期の上限がわかっている問題に対して、

$a^{r}$ に関する探索: Baby-Step
$b A^{p}$ に関する探索: Giant-Step

とにうまく平方分割して探索するアルゴリズムを Baby-Step Giant-Step 法と呼びます。離散対数問題に限らず様々な場面で活用することができ、具体的な問題例としては

などがあります。離散対数を求める実装例を以下に示します。ここでは、{ $a^{0}, a^{1}, a^{2}, \dots, a^{\sqrt{m} - 1}$ } の管理については std::map を用いた実装を示します。この場合、計算量は $O (\sqrt{m} \log m)$ になります。

#include <iostream>
#include <map>
using namespace std;

// a^b
long long modpow(long long a, long long n, long long mod) {
    long long res = 1;
    while (n > 0) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

// a^-1
long long modinv(long long a, long long m) {
    long long b = m, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= m;
    if (u < 0) u += m;
    return u;
}

// a^x ≡ b (mod. m) となる最小の正の整数 x を求める
long long modlog(long long a, long long b, int m) {
    a %= m, b %= m;

    // calc sqrt{M}
    long long lo = -1, hi = m;
    while (hi - lo > 1) {
        long long mid = (lo + hi) / 2;
        if (mid * mid >= m) hi = mid;
        else lo = mid;
    }
    long long sqrtM = hi;

    // {a^0, a^1, a^2, ..., a^sqrt(m)} 
    map<long long, long long> apow;
    long long amari = a;
    for (long long r = 1; r < sqrtM; ++r) {
        if (!apow.count(amari)) apow[amari] = r;
        (amari *= a) %= m;
    }

    // check each A^p
    long long A = modpow(modinv(a, m), sqrtM, m);
    amari = b;
    for (long long q = 0; q < sqrtM; ++q) {
        if (amari == 1 && q > 0) return q * sqrtM;
        else if (apow.count(amari)) return q * sqrtM + apow[amari];
        (amari *= A) %= m;
    }

    // no solutions
    return -1;
}

int main() {
    // 使い方の一例: a <= 10, b <= 10 に対して計算してみる
    const int MOD = 1000000007;
    for (long long a = 2; a <= 10; ++a) {
        for (long long b = 1; b <= 10; ++b) {
            cout << "log_" << a << "(" << b << ") = " << modlog(a, b, MOD) << endl;
        }
    }
}

7. 平方剰余 √a

$p$ を素数として $\sqrt{a} (\mod p)$ を求める問題です。正確には

$x^{2} \equiv a (\mod p)$

の解を求める問題です。この解は常に存在するとは限らないです。解が存在するとき $a$ は平方剰余であると言います。 $p$ が奇素数のときは、 $a = 1, 2, \dots, p - 1$ が平方剰余になっているかどうかはちょうど半々になっています。

平方剰余かどうかは「平方剰余の相互法則」によって求められる
平方剰余ならば具体的な $x$ は、Tonelli-Shanks のアルゴリズムなどによって求められる

という感じです。詳しくは kirika_comp さんの

整数論テクニック集

を読んでいただければと思います。

8. modint

ここまで $m o d . 1000000007$ の扱い方をまとめました。これらの知見を集大成した構造体 modint を紹介します。たとえば

$a = 423343$
$b = 74324$
$c = 13231$
$d = 8432455$

として、 $m o d . 1000000007$ で $(a * b + c) / d$ を計算したいとします。modint なしだと以下のような実装になります。

modintなし

// modinv は省略
const int MOD = 1000000007;

int main() {
    long long a = 423343;
    long long b = 74324;
    long long c = 13231;
    long long d = 8432455;

    cout << (a * b % MOD + c) % MOD * modinv(d, MOD) % MOD << endl;
}

しかし modint があると、こんな風に直感的な演算で実装することができます。出力結果は、どちらも 79639022 になります。

modintの使い方

const int MOD = 1000000007;
using mint = Fp<MOD>;

int main() {
    mint a = 423343;
    mint b = 74324;
    mint c = 13231;
    mint d = 8432455;

    cout << (a * b + c) / d << endl;
}

modint 構造体の作り方は C++er の方にとっては、こだわりポイントが多々ある部分だと思います。ここでは一例を示します。高速化のための工夫として

constexpr を用いる
noexcept を書くことで例外処理判定のオーバーヘッドをなくす

といったものを取り入れています。その他の様々な工夫については、noshi さんの以下の記事がとても参考になります。

modint 構造体を使ってみませんか？ (C++)

modintの実装例

#include <iostream>
using namespace std;

// modint: mod 計算を int を扱うように扱える構造体
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b; swap(a, b);
            u -= t * v; swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept {
        if (n == 0) return 1;
        auto t = modpow(a, n / 2);
        t = t * t;
        if (n & 1) t = t * a;
        return t;
    }
};


const int MOD = 1000000007;
using mint = Fp<MOD>;

int main() {
    mint a = 423343;
    mint b = 74324;
    mint c = 13231;
    mint d = 8432455;

    cout << (a * b + c) / d << endl;
}

9. おわりに: さらに先へ

ここまで述べて来たことを踏まえれば、「 $1000000007$ で割ったあまりを求めよ」というタイプの問題にひるむことはなくなると思います。Have a happy contest life!

最後に、ここまでの話でカバーできない少しマニアックな話を書き並べて行きます。数え上げ問題の途中で NTT (Number-Theoretic Transform) を用いて高速化する想定の問題では、 $\mod 1000000007$ で考えるのが大変なので「 $998244353$ で割ったあまりを求めよ」と言われたりします。

ちょっと面倒な mod として以下のケースもあります。いずれも、無駄に面倒になるという理由でコンテスト参加者から嫌われがちです。

mod が素数でない場合
mod が int 型から溢れる程度に大きい

前者は「足し算・引き算・掛け算」や「累乗」をする場合は問題ないのですが、「割り算」や「二項係数」を扱うときには注意が必要になります。対策として中国剰余定理や、マニアックな二項係数の求め方 (uwi さん著) が必要になったりします。

後者は足し算・引き算では問題ないのですが、掛け算をすると long long 型を用いていてもオーバーフローしてしまいます。解決策としては多倍長整数を扱う方法もありますが、「mod の値が大きいときの mod 演算」に示した方法も有力です。

ユーザー登録して、Qiitaをもっと便利に使ってみませんか。

あなたにマッチした記事をお届けします
ユーザーやタグをフォローすることで、あなたが興味を持つ技術分野の情報をまとめてキャッチアップできます
便利な情報をあとで効率的に読み返せます
気に入った記事を「ストック」することで、あとからすぐに検索できます

より詳しく

けんちょん (Otsuki)

@drken

NTTデータ数理システムでリサーチャーをしている大槻です。機械学習やアルゴリズムに関して面白いと思ったことを記事にしていきたいと思います。記事へのリンク等についてはお気軽にしていただいて大丈夫です。よろしくお願いします。

株式会社NTTデータ数理システム

数理科学とコンピュータサイエンスの融合！！

http://www.msi.co.jp/

この記事は以下の記事からリンクされています

アルゴ式の「シグマ計算とパイ計算」をすべて発展解法で解くからリンク2 months ago

あのアルゴリズムはどこ？[Scheme(Gauche)版]からリンク3 months ago

精進したAtCoderの問題を分類する(rating~1200)からリンク4 months ago

【AtCoder】ABC209をPython3で解説からリンク8 months ago

あのアルゴリズムはどこ？　Pythonを使用してAtCoderの緑色や水色を目指す方に、30以上のアルゴリズムスニペットと100問以上の問題（ACコード付き）を紹介！からリンク10 months ago

過去の45件を表示する

@ganyariya

2018-08-01 14:52

リンクをコピー

このコメントを報告

いつも参考にさせていただいております！
割り算 a ÷ b ～ b の逆元をかける～でよくわからない点あるため質問させていただきたいです。

逆元について、インターネットで調べて、「ある元があったとき、二項演算子を成り立たせる単位元を求めるための元。(4の掛け算の元は $1 / 4$ )」ということをこの記事を通して新たに調べて学びました。
しかし、 $b x \equiv a (m o d . p)$ の合同式を解くために、 $b^{- 1}$ を求めて掛けてあげる理由がよくわかりませんでした。掛け算をしているのは、単位元を求めて、その単位元に $a$ を求めることでオリジナルの $A$ と $B$ を求めているのでしょうか。

お時間があるときに、よろしければアドバイスいただけるとありがたいです！

(Qiitaに上げる前は写真が縦なのですが、Qiitaに上げると横向きになってしまいます。申し訳ありません。)

@drken

2018-08-01 17:26

リンクをコピー

このコメントを報告

@ganariya いつもありがとうございます！
根本的なところで誤解させてしまっていたので、全体的に書き直してみました。

a % b = 0 となる必要はないです
bx ≡ a (mod. p) を解くために b^{-1} を求める理由についてもう少し詳しく書いてみました

という感じです。まだ不明点や分かりにくい部分があれば気軽にコメントしていただければと思います。よろしくお願いします。

@SSW-SCIENTIFIC

2018-08-03 05:14

（編集済み）

リンクをコピー

このコメントを報告

TeX記法の話で恐縮ですが, modを表記するには\modや\pmod, \bmodがあります(恐らくそれぞれwith parenthesesとas binary operator)。

$a \equiv b \mod m$ => $a \equiv b \mathchoice \mod m$
$a \equiv b \pmod m$ => $a \equiv b \mathchoice (\mod m)$
$a \equiv b \bmod m$ => $a \equiv b mod m$

これらはスペーシング等を自動で調整してくれるため, \mathrmや\rmを用いるより優れています。
同様に3本線の同値記号は\equivで記号両側のスペーシングがいい感じになります。

@drken

2018-08-04 09:06

リンクをコピー

このコメントを報告

@SSW-SCIENTIFIC なるほどです！！！
ありがとうございます！！！

@drken

2018-08-07 08:36

リンクをコピー

このコメントを報告

@ganariya ガナリヤさん、いつもありがとうございます。
いつもわかりにくいところを質問していただいているおかげで、記事での説明を見直すことができています。
「割り算」について、「合同方程式」という言葉を使わずに全体的に書き直してみました。恐らく少しわかりやすくなったと思うので、読んでいただけたら嬉しいです。

@nai__

2018-08-12 13:06

リンクをコピー

このコメントを報告

大変勉強になる記事ありがとうございます．
一点，「7. 平方剰余 √a」で√aをx**2≡a (mod p)なるaの意味で用いていらっしゃいますが，「6. 離散対数 log_a b」では普通に算術的な意味で√mを用いており，かなり混乱してしまいました．mは平方数とは限らないので，例えばint(√m)を他の文字で置き換えるなどするとよりわかりやすくなるかと思いましたが，いかがでしょうか．

@drken

2018-08-27 05:00

リンクをコピー

このコメントを報告

@nai__ こちらこそコメントありがとうございます！
確かに平方剰余の √a は標語的な記述なので混乱の元かもしれないですね。
どうしたらよいか少し考えてみます。

@xuzijian629

2018-10-01 13:09

（編集済み）

リンクをコピー

このコメントを報告

こんにちは、よくhatenablog読ませてもらっています！

離散対数のところで

map<long long, long long> apow;
long long amari = 1;
for (long long r = 0; r < sqrtM; ++r) {
    apow[amari] = r;
    (amari *= a) %= m;
}

となっている箇所がありますが、正しくは

map<long long, long long> apow;
long long amari = 1;
for (long long r = 0; r < sqrtM; ++r) {
    if (!apow.count(amari)) {
        apow[amari] = r;
    }
    (amari *= a) %= m;
}

です。今のままだと、たとえば、b^x == y (mod p)の(b == 10, y == 10, p == 11)の場合で誤った答えを出力します（明らかにx == 1が解ですが、もとのアルゴリズムだと3を出力します）。

あと

// check each A^p
long long A = modpow(modinv(a, m), sqrtM, m);
amari = b;
for (long long q = 0; q <= sqrtM; ++q) {

のループはsqrtM未満まででいいと思います。

@drken

2018-10-03 08:00

リンクをコピー

このコメントを報告

@xuzijian629 ありがとうございます！！！
コメント大変ありがたいです！！！
指摘いただいた通り、離散対数のコードに間違いがありましたので修正しました！
hatenablog の件もとてもありがたいです。
今後ともよろしくお願いいたします m(_ _)m

@hydrohue2037

2020-01-22 03:02

リンクをコピー

このコメントを報告

考えてみれば当たり前のことなのですが、私はよく忘れるので

負の数も考慮した a % b を計算したい場合は、 (a % b + b) % b と計算すればよいです、ということをどっかに書いておくと親切かな、と思いました（もしどっかに書いてあったらごめんなさい）

@drken

2020-01-22 03:11

（編集済み）

リンクをコピー

このコメントを報告

@hydrohue2037 コメントありがとうございます！
確かに「負の数を考慮した a % b の計算」はそのように実現することができるのですが、二度も % 演算を行っているので処理が遅くなってしまう側面もあります。

実際は、負の数を考慮した % 演算は、「途中過程が負である状態を許しつつ計算して、最後の最後に値が負である場合には MOD を足す」とするだけでよいはずです。その方が速度的に有利なので、その方法を推したい気持ちでおります (modint では常に 0 以上になるように毎回 MOD を足しています)。

@hydrohue2037

2020-01-22 03:19

リンクをコピー

このコメントを報告

@drken なるほど、たしかにそうですね。パフォーマンスの点からは、必ずしもワンライナーが全てではないのですね。

ありがとうございました。勉強になります。

@kkato233

2020-01-23 11:22

リンクをコピー

このコメントを報告

最近 atcoder 等の競技プログラミングで TLE(計算時間オーバー) が出るので modint を探していました。
独自に作成していたのですが、どうしても割り算が高速化できなくてここにたどり着きました。

割り算のアルゴリズム採用させていただきます。ありがとうございます。

私はアルゴリズムだけ採用するので問題ないのですが、直接コードを流用したい人もいるかと思います。
今日から始めるOSSライセンス講座
にもあるように

「OSSライセンスが設定されていないソースコードには著作権法が適用され、著作権者がすべての権利を留保します。」
つまり、このソフトウェアを使用（実行・ソースコードの閲読・コンパイル）することはできますが、利用（複製・再配布・二次著作物の作成）することができないのです。

ここのコードを MIT ライセンスか GPL ライセンスにできないでしょうか？

@SSW-SCIENTIFIC

2020-01-23 15:28

リンクをコピー

このコメントを報告

@kkato233 さん

何はともあれQiitaの利用規約の第9条をお読みください。

@caprest

2020-03-16 08:53

リンクをコピー

このコメントを報告

参考になりました。
ところで記事内のリンク「pow_mod, inverse の速度」のリンクが切れてしまっているようです。

@drken

2020-03-16 08:55

リンクをコピー

このコメントを報告

@caprest
コメントありがとうございます！
確かに、リンク切れしていますね。これは、はてなグループのサービス終了に伴うものですね。
該当記事が失われてしまっていることが大変残念です...

@packer_jp

2020-10-12 05:56

リンクをコピー

このコメントを報告

離散対数のところ参考にさせていただきました。とてもわかりやすくて助かりました。その上で気づいたことがあるのでご報告をと思います。

提示されている実装例では最適解が求まらないケースがあると思います。具体的には、 $b = 1$ で $p = 0, a^{r} \equiv 1 (r > 0)$ という形の解がある場合、map には $r = 0$ のメモしかないのでこれが正しく取得できていません。 $(a, b, m) = (4, 1, 5)$ などがこれに該当します。 $x = 2$ が最適解ですが、実装例では $4$ が出力されます。

@drken

2020-11-03 19:29

リンクをコピー

このコメントを報告

@packer_jp
指摘していただいてありがとうございます！！！
おかげさまでバグを認識し、修正することができました。

@anono

2021-11-20 15:32

リンクをコピー

このコメントを報告

コメント失礼します。
4. 累乗 a^nのサンプルコードの部分なんですが、aの値が大きいときオーバーフローを起こしてしまうので、while文の前にa%=mod;を入れたほうが良いと思いました。

あなたもコメントしてみませんか :)

ユーザー登録

すでにアカウントを持っている方はログイン

「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜

1. なぜ 1000000007 で割るのか？

なぜ 1000000007 なのか？

2. 足し算、引き算、掛け算

オーバーフローに注意！！！

2-2. 引き算

3. 割り算 a ÷ b

3-1. mod p の世界における割り算とは

3-2. mod p における「逆元」

3-3. 逆元の求め方の概要

3-4. Fermat の小定理による逆元計算

3-5. 拡張 Euclid の互除法による逆元計算

3-6. 逆元が存在する条件

4. 累乗 a^n

5. 二項係数 nCr

6. 離散対数 log_a b

7. 平方剰余 √a

8. modint

9. おわりに: さらに先へ

「1000000007 で割ったあまり」の求め方を総特集！〜逆元から離散対数まで〜