ﾕﾕﾕﾕﾕ

webエンジニアです

合同式におけるモジュラ逆数 (mod_inv) の求め方

C++

　合同不定式 ax\equiv 1{\pmod {m}} を x について解く。

　 ax = 1 という不定方程式であれば、両辺に a の逆数 a^{-1} = \dfrac{1}{a} をかけて x = \dfrac{1}{a} としてあげればよい。

　同じように合同式 ax\equiv 1{\pmod {m}} においても両辺に a の逆数をかければ x が求められる。ただしこの文脈では、単に「 a の逆数」と言うよりも「 m を法とする a の逆数」と呼ぶのが正確である。単に「逆数」や「逆元」とだけ名指すことも可能とはいえ、合同式における逆数は拡張された概念であるため、単に分数のような計算を想像すると誤る。 wikipedia では「モジュラ逆数」という項目にまとまっているように、特殊な手続きを必要とする逆数であると認識しなければならない。とはいえ、例えば atcoder/ac-library においても modinv.hpp に inv() が定義され、この計算がまとめられていることから、多くのユースケースを伴う汎用的な概念であるのだろうと思っている。

　僕自身はそもそも合同式にすらさほど慣れ親しんでおらず、このあたりの概念についてはまだわからないことだらけである。だからこそ、まずはこうして学んだことを記録するところを第一歩とし、将来より高度な考察に至れるようにと願いを込めながら書く。

　以下はいくつかのソースから得た情報を僕なりに解釈したものとして読んでもらえればと思う。基本的に知識を持たないところにバラバラな資料を与えて理解を構築しているから、誤った論理もあると思う。誤謬は修正できるように未来の自分自身に期待するが、読者におかれては以上のような前提をもった上で読んでほしいと思う。

問題

　 ax\equiv 1{\pmod {m}} を x について解く。

問題文の読み替え

　問題文は次のように、 x についての合同式に変換できる。

ax\equiv 1{\pmod {m}} \Longleftrightarrow x \equiv a^{-1}\times 1{\pmod {m}}

　ここで現れる a のモジュラ逆数を x’ とおいて、さらに次のように式変換する。なお、 a と m は互いに素であるとする。

\begin{align} x’\equiv a^{-1} {\pmod {m}} &\Longleftrightarrow ax’ \equiv 1 {\pmod {m}}\\\ &\Longleftrightarrow ax’ \times my' \equiv 1 {\pmod {m}} \end{align}

　つまり一次不定式 ax’ + my’ = \gcd(a, m) を満たす任意の x’ を求めれば、それがそのまま a のモジュラ逆数 x’ のひとつの元となることにある。「ひとつの」と書いたのは、法を m とする合同式を満たす x’ は無限に存在し、ここで求めることのできる任意の x’ はその合同類のうちのひとつでしかないためである。とはいえ、解答する上では x’ は任意の一元で構わず、これは問題にならない。

一次不定式の特殊解を求める

　 ax’ + my’ = 1 を満たす任意の x’ を求めるには、適当な特殊解を見出せばよいのであって、一般解を導く必要はない。特殊解を求めるにあたって、例えば次の不定式を見てみよう。

4x + 5y = 1

これを見ると、 (x, y) = (-1, 1) と直感的にひとつの解を見出すことができるだろう。そしてこのときの x = -1 こそが、 4x \equiv 1 \pmod {5} における 4 のモジュラ逆数なのである。実際、明らかに 4 \times (-1) \equiv 1 \pmod{5} である。

　一次不定式の解法に戻ると、上のように係数が小さなケースでは単純な論理や推測で x を求めることができた。しかし、大きな二つの係数の最大公約数を解とする、次のような不定式はどうであろうか。

2021x + 1763y = \gcd(2021, 1763) = 43

これについて、直感ないし暗算で (x, y) のペアを見出すのは困難であろう。ではどうするか？

　実はここで、二つの係数 2021 と 1763 についてユークリッドの互除法を適用すると、ひとつの x, y の組み合わせが求められる。手続き的な計算手順としてはこうなる。

　まずは係数とその剰余を、剰余が発生しなくなるまで切り下げていく。これは一般的なユークリッドの互除法である。

\begin{align} 2021 &= 1763 \times 1 + 258\\\ 1763 &= 258 \times 6 + 215\\\ 258 &= 215 \times 1 + 43\\\ 215 &= 43 \times 5 \end{align}

　続いて今度は、上記のそれぞれの計算過程における剰余を根拠に、 2021x + 1763y = 43 の形を復元していく。つまり

\begin{align} 258 &= 2021 - 1763\\\ 215 &= 1763 - 258 \times 6\\\ & = 1763 - (2021 - 1763) \times 6\\\ & = -6 \times 2021 + 7 \times 1763\\\ 43 &= 258 - 215\\\ &= (2021 - 1763) - (-6 \times 2021 + 7 \times 1763)\\\ &= 7 \times 2021 - 8 \times 1763 \end{align}

　以上より 7 \times 2021 - 8 \times 1763 = 43 が求められる。つまり

\begin{align} 2021x &\equiv 1 &\pmod{1763}\\\ x &\equiv 1 \times 2021^{-1} &\pmod{1763}\\\ x &\equiv 7 &\pmod{1763} \end{align}

このように解くことができる。

アルゴリズム

　こうして確立することのできた手続きを、今度はアルゴリズムに実装していきたい。すでに材料は出そろっているから、あとはコードに置き換えていくだけ…と言いたいところだが、再帰を制御するやり方にまだしっくりきておらず、いまは写経したもので間に合わせておく。

　なお、モジュラ逆数が存在する必要十分条件は「 a と m が互いに素であること」であるが、ここでは a と m をそれぞれ \gcd(a, m) で除算してから演算を行うことで、 a と m は必ず互いに素 \Longleftrightarrow \gcd(a, m) = 1 であることが保証される。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

ll ext_gcd(ll a, ll b, ll &x, ll &y)
{
  if (b == 0)
  {
    x = 1;
    y = 0;
    return a;
  }

  ll d = ext_gcd(b, a % b, y, x);
  y -= a / b * x;
  return d;
}

ll mod(ll a, ll m)
{
  return (a % m + m) % m;
}

ll mod_inv(ll a, ll m)
{
  ll x, y;
  ll g = ext_gcd(a, m, x, y);
  return mod(x, m / g);
}

int main() {
  // 3x ≡ 1 mod 7
  cout << mod_inv(3, 7) << endl; // 5

  // 3x ≡ gcd(3, 11) = 1 mod 11
  cout << mod_inv(3, 11) << endl; // 4

  // 1071x ≡ gcd(1071, 1029) = 21 mod 1029
  cout << mod_inv(1071, 1029) << endl; // 25

  // 2021x = gcd(2021, 1763) = 43 mod 1763
  cout << mod_inv(2021, 1763) << endl; // 7

  return 0;
}