前置知识：Farey 序列以及其相关理论，根号平衡。
这里稍微提一下什么是 Farey 序列 以及其性质：
- \(F_n\) 是第 \(n\) 阶的 Farey 序列，序列中存储所有不可约的分母 \(\le n\) 的值在 \([0,1]\) 之间的分数（这里认为 \(\frac01,\frac11\) 也是不可约的分数）按照其值从小到大排序得到的有序分数序列。
- 性质 \(1\)：对于 \(F_{m-1}\) 序列上任意两个相邻分数 \(\frac ab,\frac cd\)，在 \(F_m\) 序列上一定按照顺序连续的存在分数 \(\frac ab,\frac{a+c}{b+d},\frac bd\)。
- 推论 \(1\)：对于 \(F_m\) 中任意两个相邻分数 \(\frac ab,\frac cd\)，必然有 \(bc-ad=1\)。
- 性质 \(2\)：\(|F_n|=1+\sum\limits_{i=1}^n\varphi(i)\)。
- 性质 \(3\)：对任意整数 \(n\ge 2\) 和任意实数 \(v\in[0,1]\)，一定可以在 \(F_{n-1}\) 序列中找到至少一个分数 \(\frac xy\)，满足 \(|v-\frac xy|\le \frac1{yn}\)。
- 推论 \(3\)：分数 \(\frac xy\) 一定是数 \(v\) 在序列 \(F_{n-1}\) 中向左查或者向右查能找到的第一个分数。
- 性质 \(4\)：对于 \(F_n\) 序列内的每个分数 \(\frac xy\)，\(\lfloor\frac{xn^2}{2y}\rfloor\) 的值互不相同。
回到这个题目。考虑固定 \(n\)，记 \(v=\frac ap\)。根据上面得到的性质和推论，可以知道此时必然存在一个分数 \(\frac xy\in F_{n-1}\) 满足 \(|\frac ap-\frac xy|\le\frac1{yn}\)。此时把不等式的左右两侧同时乘以 \(py\)（显然有 \(py>0\)），得到：\(|ay-px|\le\lfloor\frac pn\rfloor\)。
记 \(u=ay-px\)，则此时又能得出两个结论：
- \(ay\equiv u\pmod p\)
- \(|u|\le \lfloor\frac pn\rfloor\)。
考虑把上面提到的结论 \(1\) 改写，得到：\(y\equiv u\times a^{-1}\pmod p\)。因为这里想要求的是 \(a\) 在模 \(p\) 意义下的逆元即 \(a^{-1}\)，因此想到在同余式两侧同时乘以 \(u^{-1}\)，得到：\(a^{-1}\equiv u^{-1}\times y\pmod p\)。
然后再利用结论 \(2\) 即 \(|u|\le\lfloor\frac pn\rfloor\) 也就是 \(u\) 的绝对值很小，因此对这样的 \(a\) 只需在对应 Farey 序列上找出其对应的分数 \(\frac xy\)，计算 \(u=ay-px\)，即可套用公式 \(a^{-1}\equiv u^{-1}y\pmod p\) 解决。而此时 \(u^{-1}\) 可以直接线性求逆元。
然而此时存在 \(u<0\) 的情况。对这个东西的处理是比较容易的，有 \(p^{-1}\equiv -(-p)^{-1}\pmod p\) 然后转化成上一种情况了。
问题在于怎么快速构造 \(n\) 阶 Farey 序列 \(F_n\)。利用性质 \(2\) 可知 \(|F_n|=1+\sum\limits_{i=1}^n\varphi(i)\) 这个东西是 \(O(n^2)\) 量级的，而利用性质 \(4\) 可以想到开桶做到 \(O(V)=O(n^2)\) 排序，因此构造 \(F_n\) 的总时间复杂度为 \(O(n^2)\)。
而现在还需要找出 \(n\) 阶 Farey 序列上一个分数 \(\frac xy\) 的前驱 / 后继分数，容易想到直接在值域上维护前缀和，然后简单处理即可 \(O(n^2)-O(1)\) 查询。
因此总时间复杂度分为两个部分：
- 线性递推 \(\lfloor\frac pn\rfloor\) 以内数的逆元，时间复杂度为 \(O(\lfloor\frac pn\rfloor)\)。
- 求 \(n\) 阶 Farey 序列并在上面做一些处理，时间复杂度为 \(O(n^2)\)。
因此总时间复杂度为 \(O(n^2+\lfloor\frac pn\rfloor)\)，根号平衡一下就是 \(n^2=\lfloor\frac pn\rfloor\) 解得 \(p=n^{\frac13}\)，此时时间复杂度为 \(O(p^{\frac23}+Q)\)。
代码实现
namespace Loyalty
{int pre[N], idx, mod, n, m, Inv[N];pair<int, int> frac[N], farey[N];inline void init() { }inline void nr_init(int p){n = cbrt(p), m = n * n;for (int i = 1; i < n; ++i)for (int j = 0; j <= i; ++j){int val = j * m / i;if (!pre[val])pre[val] = 1, frac[val] = {j, i};}for (int i = 1; i <= m; ++i)pre[i] += pre[i - 1];idx = 0;for (int i = 0; i <= m; ++i)if (frac[i].first || frac[i].second)farey[idx++] = frac[i];Inv[0] = Inv[1] = 1;for (int i = 2; i <= p / n; ++i)Inv[i] = 1ll * Inv[p % i] * (p - p / i) % p;mod = p;}inline int inv(int a){int val = 1ll * a * m / mod;if (pre[val] < idx){auto [x, y] = farey[pre[val]];long long w = 1ll * a * y - 1ll * x * mod;if (abs(w) <= mod / n){int inv_t;if (w >= 0)inv_t = Inv[w];elseinv_t = (mod - Inv[-w]) % mod;return 1ll * y * inv_t % mod;}}if (pre[val] - 1 >= 0){auto [x, y] = farey[pre[val] - 1];long long w = 1ll * a * y - 1ll * x * mod;if (abs(w) <= mod / n){int inv_t;if (w >= 0)inv_t = Inv[w];elseinv_t = (mod - Inv[-w]) % mod;return 1ll * y * inv_t % mod;}}throw;return -1;}
}int p, Q, op, ol;namespace Your_Code
{int inv[20][20];inline void init(){for (int i = 2; i < 20; ++i)for (int j = 1; j < i; ++j)inv[i][j] = power(j, i - 2, i);Loyalty::nr_init(p);}inline int solve(int x){if (p < 20)return inv[p][x];return Loyalty::inv(x);}
}
参考文献
alfalfa_w 的文章 科技-在线 \(O(1)\) 逆元。