$O(n)$ 时间复杂度的解法

前言

在逛 OI WiKi 的时候发现了一个 $O(n)$ 时间复杂度的求逆元方法.
于是想到这题 Y 总给出的解法需要循环 $n$ 次每次需要用 $O(\log{n})$ 的时间运行费马小定理求逆元. 所以思考了一下能不能使用这个线性求逆元的算法预处理一下逆元, 然后就能够将整体的时间复杂度优化为 $O(n)$ 了.

运行时间对比

可以看到使用线性求逆元确实是快了.

尽管 Oi WiKi 已经有了一些证明, 但是我觉得还是说得不够清楚. 所以我尝试证明了一下这个算法. 下面来介绍它.

线性求逆元算法

要计算 $1$ 到 $p$ 所有数模 $p$ (其中, $p$ 没有模 $p$ 意义下的逆元) 的逆元时有一种线性的做法.

首先我们所以可以构造出除数 $i$, 余数 $j$, 商 $k$
$$ p = ki + j $$
则有
$$ p = ki + j \equiv 0 (\mod p) $$
将两边同时乘上逆元 $i^{-1}$ 和 $j^{-1}$ 可得
$$ kj^{-1} + i^{-1} \equiv 0 (\mod p) $$
所以有
$$ i^{-1} \equiv -kj^{-1} (\mod p) $$
又因为 $k = \lfloor \frac{p}{i} \rfloor$, $j = p \mod i$, 所以代入可得
$$ i^{-1} = (- \lfloor \frac{p}{i} \rfloor)(p \mod i)^{-1} \mod p $$
由于 $j$ 小于 $i$, 所以 $j^{-1}$ 依然先于 $i^{-1}$ 计算. 所以可以得出递推公式 :
$$ inv[i] = -(p / i) \cdot inv[p \ \% \ i] \ \% \ p $$

如果要防止计算出负值, 则可以给同余方程两边同时乘上 $1 - i$, 可得 :
$$ i^{-1} (1 - i) \equiv (1 - i) (- \lfloor \frac{p}{i} \rfloor)(p \mod i)^{-1} (\mod p) $$

展开可得

$$ i^{-1} - i^{-1}i\equiv (p - \lfloor \frac{p}{i} \rfloor)(p \mod i)^{-1} (\mod p) $$

$$ i^{-1} - 0\equiv (p - \lfloor \frac{p}{i} \rfloor)(p \mod i)^{-1} (\mod p) $$

$$ i^{-1} \equiv (p - \lfloor \frac{p}{i} \rfloor)(p \mod i)^{-1} (\mod p) $$
因此, 我们的递推公式可以写成

$$ inv[i] = (p - p / i)*inv[p \% i] \% p $$

使用限制

• 只能求 $1$ 到 $p$ 内的数在模 $p$ 意义下的逆元

线性求逆元实现代码

typedef long long LL;

// 获得 1 到 n 区间内的所有乘法逆元 (n 要小于 p, 且 p 为质数)
vector<int> getRangeModularMultInv(int n, int p) {
    vector<int> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; ++i)
        inv[i] = (LL)(p - p / i) * inv[p % i] % p;
    return inv;
}

QA

Q : 为什么只能计算 $1$ 到 $p$

A : 因为 $inv[i]$ 由 $inv[p \% i]$ 推导, 如果 $i$ 大于 $p$, 那么 $inv[i] = inv[p]$, 而 $p$ 不存在模 $p$ 意义下的逆元, 所以只能计算 $1$ 到 $p - 1$

题解

注意到题目数据 $1≤b≤a≤10^5 < MOD = 10^9 + 7$, 所以可以使用线性求逆元算法

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

typedef long long LL;

vector<int> getRangeModularMultInv(int n, int p) {
    vector<int> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; ++i)
        inv[i] = (LL)(p - p / i) * inv[p % i] % p;
    return inv;
}

const int N = 100010, MOD = 1e9 + 7;
int fac[N], invFac[N];
void calcFac(int n, int m) {
    auto invI = getRangeModularMultInv(n, m);
    fac[0] = 1; invFac[0] = 1;
    for (int i = 1; i < N; i++) {
        fac[i] = (LL)fac[i - 1] * i % m;
        invFac[i] = (LL)invFac[i - 1] * invI[i] % m;
    }
}

int getC(int a, int b, int m) {
    return (LL)fac[a] * invFac[a - b] % m * invFac[b] % m;
}

int main() {
    int q; cin >> q;
    calcFac(N, MOD);
    while (q--) {
        int a, b; cin >> a >> b;
        cout << getC(a, b, MOD) << endl;
    }
    return 0;
}

结语

证明真快乐.
没学过同余这些东西, 就自己硬证的. 不知道过程有没有错误, 欢迎大佬们指正 :)

参考文献

OI WiKi