乘法逆元

本文介紹模意義下乘法運算的逆元（Modular Multiplicative Inverse），並介紹如何使用擴展歐幾里德算法（Extended Euclidean algorithm）求解乘法逆元。

定義

如果一個線性同餘方程 \(ax \equiv 1 \pmod b\)，則 \(x\) 稱為 \(a \bmod b\) 的逆元，記作 \(a^{-1}\)。

如何求逆元

擴展歐幾里得法

實現

C++Python

void exgcd(int a, int b, int& x, int& y) {
  if (b == 0) {
    x = 1, y = 0;
    return;
  }
  exgcd(b, a % b, y, x);
  y -= a / b * x;
}

def exgcd(a, b):
    if b == 0:
        x = 1
        y = 0
        return x, y
    x1, y1 = exgcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return x, y

擴展歐幾里得法和求解線性同餘方程是一個原理，在這裏不展開解釋。

快速冪法

證明

因為 \(ax \equiv 1 \pmod b\)；

所以 \(ax \equiv a^{b-1} \pmod b\)（根據費馬小定理）；

所以 \(x \equiv a^{b-2} \pmod b\)。

然後我們就可以用快速冪來求了。

實現

C++Python

int qpow(long long a, int b) {
  int ans = 1;
  a = (a % p + p) % p;
  for (; b; b >>= 1) {
    if (b & 1) ans = (a * ans) % p;
    a = (a * a) % p;
  }
  return ans;
}

def qpow(a, b):
  ans = 1
  a = (a % p + p) % p
  while b:
      if b & 1:
          ans = (a * ans) % p
      a = (a * a) % p
      b >>= 1
  return ans

注意：快速冪法使用了費馬小定理，要求 \(b\) 是一個素數；而擴展歐幾里得法只要求 \(\gcd(a, b) = 1\)。

線性求逆元

求出 \(1,2,\dots,n\) 中每個數關於 \(p\) 的逆元。

如果對於每個數進行單次求解，以上兩種方法就顯得慢了，很有可能超時，所以下面來講一下如何線性（\(O(n)\)）求逆元。

首先，很顯然的 \(1^{-1} \equiv 1 \pmod p\)；

證明

對於 \(\forall p \in \mathbf{Z}\)，有 \(1 \times 1 \equiv 1 \pmod p\) 恆成立，故在 \(p\) 下 \(1\) 的逆元是 \(1\)，而這是推算出其他情況的基礎。

其次對於遞歸情況 \(i^{-1}\)，我們令 \(k = \lfloor \frac{p}{i} \rfloor\)，\(j = p \bmod i\)，有 \(p = ki + j\)。再放到 \(\mod p\) 意義下就會得到：\(ki+j \equiv 0 \pmod p\)；

兩邊同時乘 \(i^{-1} \times j^{-1}\)：

\(kj^{-1}+i^{-1} \equiv 0 \pmod p\)

\(i^{-1} \equiv -kj^{-1} \pmod p\)

再帶入 \(j = p \bmod i\)，有 \(p = ki + j\)，有：

\(i^{-1} \equiv -\lfloor\frac{p}{i}\rfloor (p \bmod i)^{-1} \pmod p\)

我們注意到 \(p \bmod i < i\)，而在迭代中我們完全可以假設我們已經知道了所有的模 \(p\) 下的逆元 \(j^{-1}, j < i\)。

故我們就可以推出逆元，利用遞歸的形式，而使用迭代實現：

\[ i^{-1} \equiv \begin{cases} 1, & \text{if } i = 1, \\ -\lfloor\frac{p}{i}\rfloor (p \bmod i)^{-1}, & \text{otherwise}. \end{cases} \pmod p \]

實現

C++Python

inv[1] = 1;
for (int i = 2; i <= n; ++i) {
  inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}

inv[1] = 1
for i in range(2, n + 1):
    inv[i] = (p - p // i) * inv[p % i] % p

使用 \(p-\lfloor \dfrac{p}{i} \rfloor\) 來防止出現負數。

另外我們注意到我們沒有對 inv[0] 進行定義卻可能會使用它：當 \(i | p\) 成立時，我們在代碼中會訪問 inv[p % i]，也就是 inv[0]，這是因為當 \(i | p\) 時不存在 \(i\) 的逆元 \(i^{-1}\)。線性同餘方程中指出，如果 \(i\) 與 \(p\) 不互素時不存在相應的逆元（當一般而言我們會使用一個大素數，比如 \(10^9 + 7\) 來確保它有着有效的逆元）。因此需要指出的是：如果沒有相應的逆元的時候，inv[i] 的值是未定義的。

另外，根據線性求逆元方法的式子：\(i^{-1} \equiv -kj^{-1} \pmod p\)

遞歸求解 \(j^{-1}\), 直到 \(j=1\) 返回 \(1\)。

中間優化可以加入一個記憶化來避免多次遞歸導致的重複，這樣求 \(1,2,\dots,n\) 中所有數的逆元的時間複雜度仍是 \(O(n)\)。

注意：如果用以上給出的式子遞歸進行單個數的逆元求解，目前已知的時間複雜度的上界為 \(O(n^{\frac 1 3})\)，具體請看知乎討論。算法競賽中更好地求單個數的逆元的方法有擴展歐幾里得法和快速冪法。

線性求任意 n 個數的逆元

上面的方法只能求 \(1\) 到 \(n\) 的逆元，如果需要求任意給定 \(n\) 個數（\(1 \le a_i < p\)）的逆元，就需要下面的方法：

首先計算 \(n\) 個數的前綴積，記為 \(s_i\)，然後使用快速冪或擴展歐幾里得法計算 \(s_n\) 的逆元，記為 \(sv_n\)。

因為 \(sv_n\) 是 \(n\) 個數的積的逆元，所以當我們把它乘上 \(a_n\) 時，就會和 \(a_n\) 的逆元抵消，於是就得到了 \(a_1\) 到 \(a_{n-1}\) 的積逆元，記為 \(sv_{n-1}\)。

同理我們可以依次計算出所有的 \(sv_i\)，於是 \(a_i^{-1}\) 就可以用 \(s_{i-1} \times sv_i\) 求得。

所以我們就在 \(O(n + \log p)\) 的時間內計算出了 \(n\) 個數的逆元。

實現

C++Python

s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2);
// 當然這裏也可以用 exgcd 來求逆元,視個人喜好而定.
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;

s[0] = 1
for i in range(1, n + 1):
    s[i] = s[i - 1] * a[i] % p
sv[n] = qpow(s[n], p - 2)
# 當然這裏也可以用 exgcd 來求逆元,視個人喜好而定.
for i in range(n, 0, -1):
    sv[i - 1] = sv[i] * a[i] % p
for i in range(1, n + 1):
    inv[i] = sv[i] * s[i - 1] % p

逆元練習題

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用