数论

数论

逆元

若 \(b, m\) 互质，并且 \(b | a\)，则存在一个整数 \(x\) 使得

\[\dfrac{a}{b} \equiv a \times x \pmod m \]

则称 \(x\) 为 \(b\) 的模 \(m\) 乘法逆元，记作 \(b^{-1}\pmod m\)。

例如 \(b=31, m=5, a=93\)，存在 \(x=1\) 满足 \(\dfrac{93}{31} \equiv 93 \times 1 \ (\bmod 5)\)。

因为乘法满足模运算的分配律，即 \((a \times b) \bmod m = \left( \left(a \bmod m\right) \times \left(b \bmod m\right) \right ) \bmod m\)。

而除法不满足，故乘法逆元就可以很好地解决这个问题，即

\[\dfrac{a}{b} \bmod m=((a \bmod m) \times (b^{-1} \bmod m)) \bmod m \]

这里介绍三种求逆元的方法：

1.

因为

\[\dfrac{a}{b}\equiv a \times b^{-1} \equiv \dfrac{a}{b} \times b \times b^{-1} \pmod m \]

所以

\[b \times b^{-1} \equiv 1 \pmod m \]

其实也就是要找到一个 \(x\)，满足 \(b \times x \equiv 1 \pmod m\)。

当 \(m\) 为质数时，根据费马小定理 \(b^m \equiv b \pmod m\) 可得：\(b \times b^{m-2} \equiv 1 \pmod m\)。

因此，当模数 \(m\) 为质数时，\(b^{m-2}\) 就是 \(b\) 的模 \(m\) 乘法逆元。

2.

第一种方法说过，

其实也就是要找到一个 \(x\)，满足 \(b \times x \equiv 1 \pmod m\)。

也就是解 \(my + bx = 1\) 这个方程，故我们可以使用扩展欧几里得来求解 \(x\)。

3.

这里介绍一个线性时间求 \(1 \sim n\) 在模 \(m\) 下的逆元。

我们设 \(m = q \times d + r ~ (r < d)\)。

可以得到；

1. \(q \times d + r \equiv 0 \pmod m\)，
2. \(q=\left\lfloor\dfrac{m}{d}\right\rfloor\)
3. \(r = m \bmod d\)

对于式子 1. 两边同乘 \(d^{-1} \times r^{-1}\)，得到

\[q \times r^{-1} + d^{-1} \equiv 0 \pmod m \\ d^{-1} \equiv -q \times r^{-1} \pmod m \\ d^{-1} \equiv -\left\lfloor\dfrac{m}{d}\right\rfloor \times (m \bmod d) ^ {-1} \pmod m \]

故可以根据前面的 \((m \bmod d) ^ {-1}\) 推出后面的 \(d^{-1}\)。

#include <bits/stdc++.h>
const int N = 3e6 + 10;
using namespace std;

int n, p;
long long ans[N];

int main()
{
    cin >> n >> p;
    ans[1] = 1;
    printf("%lld\n", 1);

    for (int i = 2; i <= n; ++ i)
        ans[i] = ((- (p / i) * ans[p % i]) % p + p) % p,
        printf("%lld\n", ans[i]);

    return 0;
}

扩展欧几里得定理

问题引入

我们要解决一个这样的问题：

\[ax+by=\gcd(a, b) \]

分析

在欧几里得算法的最后一步，即 \(b=0\) 时，显然有一对整数 \(x=1, y=0\) 满足 \(a\times 1+b \times 0=\gcd(a, 0)\)。

我们假设有 \(x_2, y_2\) 满足 \(bx_2+(a \bmod b)y_2=\gcd(b, a\bmod b)\)。

由于 \(\gcd(a, b)=\gcd(b, a\bmod b)\)，故 \(ax+by=bx_2+(a \bmod b)y_2\)

又因为 \((a \bmod b) = a - \left\lfloor\dfrac{a}{b}\right\rfloor \times b\)。

故

\[ax+by=bx_2+(a - \left\lfloor\dfrac{a}{b}\right\rfloor b)y_2 \\ ax+by=bx_2+ay_2-\left\lfloor\dfrac{a}{b}\right\rfloor by_2 \\ ax+by=ay_2+b(x_2-\left\lfloor\dfrac{a}{b}\right\rfloor y_2) \]

由于 \(a=a, b=b\)，故 \(x=y_2, y=x_2-\left\lfloor\dfrac{a}{b}\right\rfloor y_2\)，这样从 \(x=1, y=0\) 一步一步就可以推到 \(ax+by=\gcd(a, b)\) 的解。

对于一般的方程 \(ax+by=c\)，当且仅当 \(\gcd(a,b) | c\) 时有解。

我们可以先求出 \(ax+by=\gcd(a, b)\) 的特殊解 \(x_0, y_0\)，然后令

\[x = x_0 \times \dfrac{c}{\gcd(a, b)}, y=y_0\times \dfrac{c}{\gcd(a, b)} \]

此时 \(ax+by=c\)，\(x, y\) 可能为负数，有的题目会要求 \(x, y\) 为正整数，我们可以把原式进行修改，变为

\[ax+by+k\dfrac{ab}{\gcd(a, b)}-k\dfrac{ab}{\gcd(a, b)}=c \\ a(x+k\dfrac{b}{\gcd(a, b)})+b(y-k\dfrac{a}{\gcd(a, b)})=c \]

即

\[x=x+k\dfrac{b}{\gcd(a, b)} \\ y=y-k\dfrac{a}{\gcd(a, b)} \]

其中 \(k \in \Z\)。

类欧几里得算法

问题引入

设

\[f(a, b, c, n) = \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor \]

其中 \(a, b, c, n\) 是常数，需要 \(\mathcal O(\log n)\) 的做法。

若 \(a \geq c\) 或 \(b \geq c\)，我们可以将 \(a, b\) 对 \(c\) 取模以简化问题。

考虑到

\[x = \left \lfloor \frac{x}{c}\right \rfloor c + x \bmod c \]

故

\[\begin{split} f(a, b, c, n) &= \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor \\ &= \sum_{i=0}^n \left\lfloor\frac{(\left \lfloor \frac{a}{c}\right \rfloor c + a \bmod c)i + (\left \lfloor \frac{b}{c}\right \rfloor c + b \bmod c)}{c}\right\rfloor \\ &= \frac{n(n + 1)}{2} \left \lfloor \frac{a}{c}\right \rfloor + (n+1) \left \lfloor \frac{b}{c}\right \rfloor + f(a \bmod c , b\bmod c, c, n) \end{split} \]

此时一定有 \(a < c\) 且 \(b < c\)。

令

\[S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \]

再进行转化

\[\begin{split} \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor &= \sum_{i=0}^n \sum_{j=0}^{S(i)} 1 \\ &= \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ j \leq S(i) \right] \end{split} \]

考虑到

\[\begin{split} j \leq S(i) &\iff j + 1 \leq {\left\lfloor\frac{ai + b}{c}\right\rfloor} \\ j + 1 \leq {\left\lfloor\frac{ai + b}{c}\right\rfloor} &\iff j + 1 \leq \frac{ai + b}{c} \\ j + 1 \leq \frac{ai + b}{c} &\iff jc + c \leq ai + b \\ jc + c \leq ai + b &\iff jc + c - b \leq ai \\ jc + c - b \leq ai &\iff jc + c - b - 1 < ai \\ jc + c - b - 1 < ai &\iff \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor < i \end{split} \]

故

\[\begin {split} \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ j \leq S(i) \right] &= \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ i > \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor\right] \\ &= \sum_{j=0}^{S(n)} \left(n - \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \right) \\ &= (S(n) + 1)n - \sum_{j=0}^{S(n)} \left \lfloor \frac{c j + (c - b - 1)}{a} \right \rfloor \\ &= (S(n) + 1)n - f(c, c - b - 1, a, S(n)) \end {split} \]

故

\[f(a, b, c, n) = (S(n) + 1)n - f(c, c - b - 1, a, S(n)) \]

可以发现，上述式子是一个递归式，我们不断重复上述过程，先取模，后递归，其实就是辗转相除的过程，时间复杂度 \(\mathcal O(\log n)\)。

拓展

我们再来推导两个变种求和式

\[g(a, b, c, n) = \sum_{i=0}^n i\left\lfloor\frac{ai + b}{c}\right\rfloor \\ h(a, b, c, n) = \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor^2 \]

推导 \(g\)

引理 1.

\[\sum_{i=0}^n i^2 = \frac{n (n+1)(2n+1)}{6} \]

证明如下：

考虑到

\[(n+1)^3=n^3+3n^2+3n+1 \]

故

\[\begin {split} (n+1)^3-n^3 &= 3n^2+3n+1 \\ n^3-(n-1)^3 &=3(n-1)^2+3(n-1)+1 \\ &\cdots \\ 2^3 - 1^3&=3 \times (2-1)^2+3 \times(2-1) + 1 \\ \end {split} \]

将这 \(n\) 个等式左右两边相加，得到

\[(n+1)^3 - 1 = 3(1^2 + 2^2 + \cdots +n^2) + 3(1 + 2 + \cdots +n) + n \]

即

\[n^3+3n^2+3n = 3(1^2 + 2^2 + \cdots +n^2) + 3\frac{n(1 +n)}{2} + n \]

整理后得

\[\sum_{i=0}^n i^2 = \frac{n (n+1)(2n+1)}{6} \]

首先和 \(f\) 一样，对其取模（根据引理 1. 可将其展开）。

\[\begin{split} g(a, b, c, n) &= \sum_{i=0}^n i\left\lfloor\frac{ai + b}{c}\right\rfloor \\ &= \sum_{i=0}^n i\left\lfloor\frac{(\left \lfloor \frac{a}{c}\right \rfloor c + a \bmod c)i + (\left \lfloor \frac{b}{c}\right \rfloor c + b \bmod c)}{c}\right\rfloor \\ &= \frac{n(n + 1)(2n+1)}{6} \left \lfloor \frac{a}{c}\right \rfloor + \frac{n(n + 1)}{2} \left \lfloor \frac{b}{c}\right \rfloor + g(a \bmod c , b\bmod c, c, n) \end{split} \]

其他部分推导与 \(f\) 类似。

令

\[S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \]

有

\[\begin {split} g(a, b, c, n) &= \sum_{j=0}^{S(n)} \sum_{i=0}^n i\left[ j \leq S(i) \right] \\ &= \sum_{j=0}^{S(n)} \sum_{i=0}^n i\left[ i > \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor\right] \\ \end {split} \]

令

\[t(j) =\left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \]

则有

\[\begin{split} g(a, b, c, n) &= \sum_{j=0}^{S(n)} \left((t(j) + 1) + (t(j) + 2) + \cdots + n\right) \\ &= \sum_{j=0}^{S(n)} \left(\frac{(t(j)+1+n)\times(n-t(j))}{2}\right) \\ &= \frac{1}{2}\left[ (S(n) + 1) n(n+1) - \sum_{j=0}^{S(n)}(t(j))^2 - \sum_{j=0}^{S(n)}t(j)\right] \\ &= \frac{1}{2} \left[(S(n) + 1) n(n+1) - h(c, c - b - 1, a, S(n)) - f(c, c - b - 1, a, S(n))\right] \end {split} \]

推导 \(h\)

同样套路，先取模

\[\begin{split} h(a, b, c, n) &= \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor^2 \\ &= \sum_{i=0}^n \left\lfloor\frac{(\left \lfloor \frac{a}{c}\right \rfloor c + a \bmod c)i + (\left \lfloor \frac{b}{c}\right \rfloor c + b \bmod c)}{c}\right\rfloor^2 \\ &= \sum_{i=0}^n \left(i \left \lfloor \frac{a}{c}\right \rfloor + \left \lfloor \frac{b}{c}\right \rfloor + \frac{(a \bmod c)i + (b\bmod c)}{c} \right)^2 \\ \end{split} \]

令

\[Q(i)=\frac{(a \bmod c)i + (b\bmod c)}{c} \]

拆开则有

\[\begin {split} h(a, b, c, n) &= \sum_{i=0}^n \left(i^2\left \lfloor \frac{a}{c}\right \rfloor^2 + \left \lfloor \frac{b}{c}\right \rfloor^2 + (Q(i))^2 + 2 i\left \lfloor \frac{a}{c}\right \rfloor \left \lfloor \frac{b}{c}\right \rfloor + 2i\left \lfloor \frac{a}{c}\right \rfloor Q(i) + 2\left \lfloor \frac{b}{c}\right \rfloor Q(i)\right) \\ &= \frac{n(n + 1)(2n + 1)}{6} \left \lfloor \frac{a}{c}\right \rfloor^2 + (n + 1)\left \lfloor \frac{b}{c}\right \rfloor^2 + h(a \bmod c, b \bmod c, c, n) + \\ & n(n+1)\left \lfloor \frac{a}{c}\right \rfloor \left \lfloor \frac{b}{c}\right \rfloor + 2\left \lfloor \frac{a}{c}\right \rfloor g(a \bmod c, b \bmod c, c, n) + 2\left \lfloor \frac{b}{c}\right \rfloor f(a \bmod c, b \bmod c, c, n) \end {split} \]

照样令

\[S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \\ t(j) =\left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \]

发现平方不好处理，将平方转化为加法

\[\begin{split} n^2 &= 2\frac{n(n + 1)}{2} - n \\ &= \left(2\sum_{i=0}^n i \right) - n \end{split} \]

则

\[\begin{split} h(a, b, c, n) &= \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor^2 \\ &= \sum_{i=0}^n \left[\left(2\sum_{j=0}^{S(i) + 1} j \right) - \left\lfloor\frac{ai + b}{c}\right\rfloor\right] \\ &= \left(2 \sum_{i=0}^n \sum_{j=0} ^{S(i) + 1} j\right) - f(a, b, c, n) \end{split} \]

现在只要解决前面的式子即可

\[\begin{split} \sum_{i=0}^n \sum_{j=0} ^{S(i) + 1} j &= \sum_{i=0}^n \sum_{j=0} ^{S(i)} (j+1) \\ &= \sum_{j=0}^{S(n)} (j+1)\sum_{i=0}^n \left[j < \left\lfloor\frac{ai + b}{c}\right\rfloor \right] \\ &= \sum_{j=0}^{S(n)} (j+1) (n - t(j)) \\ &= (S(n) + 1)n + \frac{S(n)(S(n) + 1)}{2} n - g(c, c - b - 1, a, S(n)) - f(c, c - b - 1, a, S(n)) \end{split} \]

综上所述

\[h(a, b, c, n) =(S(n)+1)(S(n) + 2) n - 2g(c, c - b - 1, a, S(n)) - 2f(c, c - b - 1, a, S(n)) - f(a, b, c, n) \]

当 \(a=0\) 时，上述三个式子都可以 \(\mathcal O(1)\) 计算。

代码实现时，因为三个函数交错递归，考虑三个一起整体递归，同步计算，时间复杂度为 \(\mathcal O(\log n)\)。

Luogu5170

#include <bits/stdc++.h>
const int mod = 998244353;
int InvTwo, InvSix;
inline int read()
{
    int cnt = 0; char ch = getchar(); bool op = 1;
    for (; ! isdigit(ch); ch = getchar())
        if (ch == '-') op = 0;
    for (; isdigit(ch); ch = getchar())
        cnt = cnt * 10 + ch - 48;
    return op ? cnt : - cnt;
}

inline int quick_pow(int a, int b)
{
    int ret = 1;
    for (; b; b >>= 1)
    {
        if (b & 1) ret = 1ll * ret * a % mod;
        a = 1ll * a * a % mod;
    }
    return ret % mod;
}

struct eu
{
    int f, g, h;
};

inline int S(int a, int b, int c, int i)
{
    return ((1ll * a * i + b) / c) - 1;
}

inline eu solve(int a, int b, int c, int n)
{
    if (a == 0)
    {
        eu now = {0, 0, 0};
        now.f = 1ll * (b / c) * (n + 1) % mod;
        now.h = 1ll * (b / c) * (b / c) % mod * (n + 1) % mod;
        now.g = 1ll * n * (n + 1) % mod * InvTwo % mod * (b / c) % mod;
        return now;
    }
    if (a >= c || b >= c)
    {
        eu now = {0, 0, 0};
        now.f = (now.f + 1ll * n * (n + 1) % mod * InvTwo % mod * (a / c) % mod) % mod;
        now.f = (now.f + 1ll * (n + 1) * (b / c) % mod) % mod;
        now.g = (now.g + 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * InvSix % mod * (a / c) % mod) % mod;
        now.g = (now.g + 1ll * n * (n + 1) % mod * InvTwo % mod * (b / c) % mod) % mod;
        now.h = (now.h + 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * InvSix % mod * (a / c) % mod * (a / c) % mod) % mod;
        now.h = (now.h + 1ll * (n + 1) * (b / c) % mod * (b / c) % mod) % mod;
        now.h = (now.h + 1ll * n * (n + 1) % mod * (a / c) % mod * (b / c) % mod) % mod;
        eu nxt = solve(a % c, b % c, c, n);
        now.f = (now.f + nxt.f) % mod;
        now.g = (now.g + nxt.g) % mod;
        now.h = (now.h + nxt.h) % mod;
        now.h = (now.h + 1ll * 2 * (a / c) % mod * nxt.g % mod) % mod;
        now.h = (now.h + 1ll * 2 * (b / c) % mod * nxt.f % mod) % mod;
        return now;
    }
    else
    {
        eu now = {0, 0, 0};
        eu nxt = solve(c, c - b - 1, a, S(a, b, c, n));
        now.f = 1ll * (S(a, b, c, n) + 1) * n % mod;
        now.f = (now.f - nxt.f) % mod;
        now.f = (now.f + mod) % mod;
        now.g = 1ll * InvTwo * ((((1ll * (S(a, b, c, n) + 1) * n % mod * (n + 1) % mod - nxt.h) % mod) - nxt.f) % mod) % mod;
        now.g = (now.g + mod) % mod;
        now.h = (1ll * (S(a, b, c, n) + 1) * (S(a, b, c, n) + 2) % mod * n % mod) % mod;
        now.h = (now.h - 1ll * 2 * nxt.g) % mod;
        now.h = (now.h - 1ll * 2 * nxt.f) % mod;
        now.h = (now.h - now.f) % mod;
        now.h = (now.h + mod) % mod;
        return now;
    }
}

int main()
{
    int t = read();
    InvTwo = quick_pow(2, mod - 2);
    InvSix = quick_pow(6, mod - 2);
    while (t --)
    {
        int n, a, b, c;
        n = read(), a = read(), b = read(), c = read();
        eu ans = solve(a, b, c, n);
        printf("%d %d %d\n", ans.f, ans.h, ans.g);
    }
    return 0;
}

Reference

oi-wiki

中国剩余定理（CRT）

问题引入

我们需要解决满足下列方程组的一个解 \(x\)，并且保证所有的 \(b\) 都互相互质。

\[\begin{cases}x \equiv a_1 \ (\bmod \ b_1) \\ x \equiv a_2 \ (\bmod \ b_2) \\ ~~~~~~~~~~~~~\vdots \\ x \equiv a_n \ (\bmod \ b_n)\end{cases} \]

分析

我们考虑设 \(M=\prod\limits_{i=1}^{n} b_i\)，\(m_i=\dfrac{M}{b_i}\)，\(t_i\) 为 \(m_ix\equiv 1\pmod {b_i}\) 的一个解，则解为：

\[x=\sum\limits_{i=1}^{n}a_im_it_i \]

我们现在来证明一下 \(x\) 是一个满足条件的解。

证明：

对于任意的 \(k \in [1, n]\)，\(a_km_kt_k \equiv a_k \pmod {b_k}\)，

而 \(a_km_kt_k \equiv 0 \pmod {b_i}(i \in [1, k - 1] \ \bigcup \ [k + 1, n])\)，故 \(x\) 是一个满足条件的解。

我们知道，\(M \equiv 0 \pmod {b_i} (i \in [1, n])\)，故 \(x + kM \ (k \in \Z)\) 都是满足答案的解，若题目要求最小正整数解，则：

\[x = (x \bmod M + M) \bmod M \]

#include <bits/stdc++.h>
const int N = 10 + 10;
using namespace std;

int n;
long long M = 1, m[N], t[N], ans;
int a[N], b[N];

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

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++ i)
        cin >> a[i] >> b[i], M = M * a[i];

    for (int i = 1; i <= n; ++ i)
    {
        long long y;
        m[i] = M / a[i];
        exgcd(m[i], a[i], t[i], y);
        ans = (ans + m[i] * t[i] * b[i]) % M;
    }

    cout << (ans + M) % M << endl;

    return 0;
}

拓展中国剩余定理 （EXCRT）

问题引入

我们要解决的还是这么上文那个问题，其中不保证所有的 \(b\) 都互相互质。

分析

我们假设已经有一个解 \(x\) 满足前 \(k - 1\) 个同余方程，此时我们设 \(d=\operatorname{lcm}\{b_1, b_2, \cdots, b_{k-1}\}\)，显然，\(x + jd \ (j \in \Z)\) 都是满足条件的解。

那么对于第 \(k\) 个同余方程，即要满足 \(x+jd\equiv a_k \pmod {b_k}\)，可转换为 \(jd \equiv a_k - x \ (\bmod \ b_{k})\)，显然，可以用扩展欧几里得来求解 \(j\)，若 \(\gcd(d, b_k)\nmid a_k-x\) 那么就无解，否则 \(x = x+jd, d=\operatorname{lcm}\{b_1, b_2, \cdots, b_k\}\)，此时的 \(x\) 是满足前 \(k\) 个方程的，这样一直做到最后就可以求出答案。

#include <bits/stdc++.h>
const int N = 1e5 + 10;
using namespace std;

int n; long long X, m;
long long a[N], b[N];

inline long long exgcd(long long a, long long b, long long &x, long long &y)
{
    if (b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    long long d = exgcd(b, a % b, x, y);
    long long t = x; x = y; y = t - (a / b) * y;
    return d;
}

inline long long quick_mul(long long a, long long b, long long p)
{
    a %= p, b %= p;
    long long res = 0, op = 1;
    if (b < 0) b = -b, op = -1;
    for (; b; b >>= 1)
    {
        if (b & 1) res = (res + a) % p;
        a = (a + a) % p;
    }
    return res * op;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++i)
        cin >> a[i] >> b[i];
    X = b[1], m = a[1];

    for (int i = 2; i <= n; ++i)
    {
        long long d, x, y, c = (b[i] - X);
        d = exgcd(m, a[i], x, y);
        x = quick_mul(x, (c / d), (a[i] / d));
        x = (x % (a[i] / d) + (a[i] / d)) % (a[i] / d);
        X = (X + x * m);
        m = m / d * a[i];
        X = (X % m + m) % m;
    }

    cout << X << endl;

    return 0;
}

BSGS

问题引入

我们需要解决 \(a^x \equiv b \pmod M\) 的一个解 \(x\)，其中 \(a, M\) 互质。

分析

首先，如果有解，那么解一定在 \([0, M]\) 这个区间内，根据抽屉原理，\(a^x \bmod M\) 只会有 \(M\) 个取值，故 \([0, M]\) 中至少会存在两个相同余数的数，即一定存在一对 \(i, j\) 满足 \(a^i\equiv a^j \pmod M\)，那么在 \(i\) 之后，\(a^i \bmod \ M, a^{i + 1} \bmod \ M, \cdots a^{j} \bmod \ M\) 这些数会循环出现，故如果有解，那么解一定在 \([0, M]\) 这个区间内。

我们考虑把 \(x\) 进行拆分，\(x=qm-p\)，其中 \(m=\left\lfloor\sqrt{x}\right\rfloor+1\)，\(q \in [1, m]\)，\(p \in [0, m]\)。

显然，\(qm-p\) 可以表示出 \(0 \sim m\) 中的所有数。

我们把原式修改一下得到 \(a^{qm-p} \equiv b \pmod M\)，\(a^{qm} \equiv b \times a^{p} \pmod M\)（该步后往前推导需要有乘法逆元）

我们发现 \(q\) 与 \(p\) 没有关系了，考虑先 \(\mathcal O(\sqrt{x})\) 枚举 \(p\)，把 \((b \times a^p) \bmod M\) 放入哈希表（map 也行），然后再 \(\mathcal O(\sqrt{x})\) 枚举 \(a^{qm}\)，找到最小的 \(q\) 满足 \(a^{qm} \bmod M\) 在哈希表中出现过，那么 \(qm+p\) 就是最小的答案。

#include <bits/stdc++.h>

int p, b, n, m;
std::map <int, int> mp;

inline int quick_pow(int x, long long y)
{
    int res = 1; x %= p;
    for (; y; y >>= 1)
    {
        if (y & 1) res = 1ll * res * x % p;
        x = 1ll * x * x % p;
    }
    return res;
}

int main()
{
    std::cin >> p >> b >> n;

    m = sqrt(p) + 1;
    for (int i = 0; i <= m; ++ i)
        mp[1ll * n * quick_pow(b, i) % p] = i;
    for (int i = 1; i <= m; ++ i)
    {
        if (mp[quick_pow(b, 1ll * i * m)])
        {
            std::cout << 1ll * i * m - mp[quick_pow(b, i * m)] << std::endl;
            return 0;
        }
    }
    std::cout << "no solution" << std::endl;

    return 0;
}

扩展 BSGS

问题引入

我们要解决的还是 \(a^x\equiv b \pmod p\) 这个问题，只是不满足 \(a\) 与 \(p\) 互质。

分析

若 \(a\) 与 \(p\) 互质，那么就会存在逆元，可以将 \(x\) 分解为 \(qm - p\) 来求解，所以我们要想办法让 \(a, p\) 变得互质。

我们考虑把这个式子转化为不定方程，\(a^x -yp = b\)，首先，我们求出 \(d_1=\gcd(a, p)\)，然后整个式子都除以 \(d_1\)，若 \(b\) 与 \(d_1\) 互质，那么就无解，否则我们继续求出 \(d_2=\gcd(a, p)\)，一直进行下去，直到 \(a\) 与 \(p\) 互质。

假设进行了 \(k\) 次，最后得到式子就是：

\[\dfrac{a^k}{\prod\limits_{i=1}^{k}d_i}\times a^{x-k}-\dfrac{yp}{\prod\limits_{i=1}^{k}d_i}=\dfrac{b}{\prod\limits_{i=1}^{k}d_i} \]

其中 \(a\) 与 \(p\) 互质。

我们再化为同余方程的形式：

\[\frac{a^k}{\prod\limits_{i=1}^{k}d_i}\times a^{x-k}\equiv \frac{b}{\prod\limits_{i=1}^{k}d_i} \left( \bmod {\frac{p}{\prod\limits_{i=1}^{k}d_i}}\right) \]

显然，求出

\[\dfrac{a^k}{\prod\limits_{i=1}^{k}d_i} \]

的逆元，然后直接移到右边，直接用 BSGS 来解决，然后最后把答案加上 \(k\) 即可。

注意，不保证答案小于 \(k\) 的情况，故在消因子前做一下 \(O(k)\) 枚举，直接验证 \(a^i \equiv b \pmod p\)，就能避免这种情况。

code

高斯消元

问题引入

给出 \(n\) 组方程，每组方程形如 \(\sum\limits_{i=1}^n a_ix_i = b_i\)，要求求出 \(x_1, x_2 \cdots x_n\) 或告知无解。

分析

我们把这些方程转化为矩阵

\[\begin{vmatrix} & a_{1,1} & a_{1,2} & \cdots & a_{1, n} & b_1 \\ & a_{2,1} & a_{1,2} & \cdots & a_{2, n} & b_2 \\ & \vdots & \vdots & \ddots & \vdots & \vdots & \\ & a_{n,1} & a_{n,2} & \cdots & a_{n, n} & b_n \\ \end{vmatrix} \]

我们希望最后的矩阵形如

\[\begin{vmatrix} & 1 & 0 & 0 & 0 & \cdots & 0 & c_1 \\ & 0 & 1 & 0 & 0 & \cdots & 0 & c_2 \\ & 0 & 0 & 1 & 0 & \cdots & 0 & c_3 \\ & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \\ & 0 & 0 & 0 & \cdots & 1 & 0 & c_{n-1} \\ & 0 & 0 & 0 & \cdots & 0 & 1 & c_n \\ \end{vmatrix} \]

对于第 \(i\) 行就能求出 \(x_i = c_i\)，即可求出解。

流程如下：

1. 我们从上往下做，先将当前系数绝对值最大的移到当前这一行，减小误差。
2. 让当前系数化为 \(1\)，即让当前行的每一个数都除以当前系数。
3. 再用当前第 \(i\) 行的第 \(i\) 个系数（即为 \(1\)）消去 $j \ (j \not= i) $ 的第 \(i\) 个系数，把它消为 \(0\)。

重复做以上的操作，直到做到第 \(n\) 行，此时最后的矩阵形如上述。

我们可以在 \(\mathcal O(n^3)\) 的时间复杂度完成。

判断解情况：

最后说一下如何判无解和无限解：

1. 无解：当前行系数全为 \(0\)，但值不为 \(0\)。

2. 无限解：当前行系数全为 \(0\)，且值为 \(0\)。

一、积性函数

1. 积性函数：

对于任意互质的整数 \(a\) 和 \(b \ (\gcd(a,b)=1)\) 有性质 \(f(ab)=f(a)\times f(b)\)。

2. 完全积性函数：

对于任意的整数 \(a\) 和 \(b\) 有性质 \(f(ab) = f(a) \times f(b)\)。

3. 证明 \(\varphi(n)\) 是积性函数。

1. 定理：

当 \(\gcd(x, y)=1\) 时，\(\varphi(x \times y)=\varphi(x) \times \varphi(y)\)。

2. 证明：

我们将 \(1 \sim x \times y\) 排列成如下方阵：

\[ \begin{matrix} 1 & 2 &\cdots & i & \cdots & y \\ 1 + y & 2 + y & \cdots & i + y & \cdots & 2 \times y \\ 1 + 2 \times y & 2 + 2 \times y & \cdots & i + 2 \times y & \cdots & 3 \times y \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 1 + j \times y & 2 + j \times y & \cdots & i+j \times y & \cdots & (j + 1) \times y \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 1 + (x-1) \times y & 2 + (x-1) \times y & \cdots & i + (x-1) \times y & \cdots & x \times y \\ \end{matrix} \]

由于 \(\gcd(x, y)=1\)，每一列 \(i + j \times y\) 都能构成一个 \(\bmod y\) 的完全剩余系 \(\{\overline0, \overline1, \cdots, \overline{y-2}， \overline{y-1}\}\)，则每一列与 \(y\) 互质的数就有 \(\varphi(y)\)。同理，对于每一行，都能构成一个 \(\bmod x\) 的完全剩余系 \(\{\overline0, \overline1, \cdots, \overline{x-2}， \overline{x-1}\}\)，并且它们都是对齐的，则每一行与 \(x\) 互质的数就有 \(\varphi(x)\)。

若 \(a\) 与 \(x \times y\) 互质，则 \(a\) 与 \(x\) 互质且 \(a\) 与 \(y\) 互质，即 \(\varphi(x \times y)=\varphi(x) \times(y)\)。

二、欧拉函数

1. 欧拉函数的定义及求法

\(\varphi(x)\) 表示有多少个 \(i \ (1 \leq i \leq x)\) 满足 \(\gcd(i, x) = 1\)（即有多少个 \(\leq x\) 的 \(i\) 与 \(x\) 互质）。

显然当 \(x\) 为质数时，\(\varphi(x) = x - 1\)。

把 \(x\) 唯一分解：

\[x=\prod\limits_{i=1}^{m}p_i^{c_i} \]

有

\[\varphi(x)=x \times \prod\limits_{i=1}^{m}\dfrac{p_i-1}{p_i} \]

我们来证明一下这个定理：

引理 1：设 $x=p^k \ (p $ 为质数\()\)，那么 \(\varphi(x)=p^{k-1}\times(p-1)\)。

证明 1：此时对于一个 \(y\) 与 \(x\) 互质，当且仅当 \(y \not\equiv 0 \pmod p\)，我们考虑将 \(x\) 分成许多段，长度为 \(p\)，有 \(\frac{x}{p}=\frac{p^k}{p}=p^{k-1}\) 段。其中每一段都有 \(p-1\) 个数与 \(p\) 互质，有 \(p^{k-1}\) 段，故 \(\varphi(x)=p^{k-1}\times(p-1)\)。

根据欧拉函数的积性性质，我们可以得到：

\[\varphi(x)=\prod\limits_{i=1}^{m}\varphi(p_i^{c_i})=\prod\limits_{i=1}^mp^{c_i-1}\times(p_i-1)=\prod\limits_{i=1}^m\dfrac{p^{c_i}}{p}\times(p_i-1)=\prod\limits_{i=1}^mp^{c_i}\times\dfrac{p_i-1}{p_i}=x\times\prod\limits_{i=1}^m\dfrac{p_i-1}{p_i} \]

故我们可以使用试除法在 \(O(\sqrt{x})\) 的时间内求出 \(\varphi(x)\)。

2. 欧拉函数的一些性质

定理 1：\(\forall \ n > 1\)，\(1 \sim n\) 之间与 \(n\) 互质的数的和为：\(\dfrac{n \times \varphi(n)}{2}\)。

证明：因为 \(\gcd(n, x)=\gcd(n, n-x)\)，故与 \(n\) 不互质的数是成对出现的，且平均值为 \(\dfrac{n}{2}\)，又考虑到共有 \(\varphi(n)\) 个与其互质的数，故 \(\forall \ n > 1\)，\(1 \sim n\) 之间与 \(n\) 互质的数的和为 \(\dfrac{n \times \varphi(n)}{2}\)。

定理 2：当 \(x\) 为奇数时，\(\varphi(2 \times x)=\varphi(x)\)。

证明 2：把 \(x\) 唯一分解：

（其中 \(p_i \not=2\)），\(2\times x=2\times\prod\limits_{i=1}^mp_i^{c_i}\)，根据上面的欧拉函数的求法可得：

\[\varphi(2\times x)=2 \times x \times \dfrac{2-1}{2} \times \prod\limits_{i=1}^m\dfrac{p_i - 1}{p_i}=x\times \prod\limits_{i=1}^m\dfrac{p_i-1}{p_i} \\ \varphi(x)=x\times \prod\limits_{i=1}^m\dfrac{p_i-1}{p_i} \]

所以当 \(x\) 为奇数时，\(\varphi(2\times x)=\varphi(x)\)。

定理 3：若 \(x > 2\) ，那么 \(\varphi(x)\) 是偶数。

证明 3：我们考虑把 \(x\) 唯一分解：

\[x=\prod\limits_{i=1}^{m}p_{i}^{c_i} \\ \varphi(x)=x\times\prod\limits_{i=1}^m\dfrac{p_i-1}{p_i} \]

显然，除开 \(2\) 这个质数以外，其他质数都是奇数，减去 \(1\) 之后就都是偶数了，故不包含 \(2\) 这个质因子的数的欧拉函数一定是偶数，而 \(2-1=1\)，对结果没有影响，并且如果这个数有 \(2\) 这个质因子，那么它一定是偶数，所以若 \(x > 2\) ，那么 \(\varphi(x)\) 是偶数。

定理 4：设 \(p\) 为质数，若 \(p | n\) 且 \(p^2 | n\) 则 \(\varphi(n)=\varphi(\dfrac{n}{p})\times p\)。

证明 4：由于 \(p^2 | n\) 所以 \(\dfrac{n}{p}\) 也有 \(p\) 这个因子，故 \(\varphi(n)\) 和 \(\varphi(\dfrac{n}{p})\) 的区别只有 \(n\) 和 \(\dfrac{n}{p}\) 这个乘积的差别，故可得 \(\varphi(n)=\varphi(\dfrac{n}{p})\times p\)。

定理 5：设 \(p\) 为质数，若 \(p | n\) 且 \(p^2 \nmid n\) 则 \(\varphi(n)=\varphi(\dfrac{n}{p})\times(p-1)\)。

证明 5：由于 \(p^2 \nmid n\) 所以 \(\dfrac{n}{p}\) 与 \(p\) 互质，故 \(\varphi(n)\) 和 \(\varphi(\dfrac{n}{p})\) 的比，多了 \(\dfrac{p-1}{p}\) 和 \(p\) 这两个乘积，故可得：\(\varphi(n)=\varphi(\dfrac{n}{p})\times (p-1)\)。

定理 6：

\[x=\sum\limits_{d|x}\varphi(d) \]

证明：设 \(f(n)=\sum\limits_{d|n}\varphi(d)\)， 根据乘法分配律以及 \(\varphi\) 的积性性质，我们可以得到：若 \(n, m\) 互质，则

\[f(nm) = \sum\limits_{d | nm} \varphi(d)=\left(\sum\limits_{d|n}\varphi(d)\right)\times \left(\sum\limits_{d|m}\varphi(d) \right)=f(n)\times f(m) \]

故 \(f(n)\) 是积性函数。对于单个质因子 ：

\[f(p_i^{c_i})=\sum\limits_{d|p_i^{c_i}}\varphi(d)=\varphi(1)+\varphi(p)+\cdots+\varphi(p_i^{c_i}) \]

我们发现这是一个等比数列求和再加 \(1\)，结果为 \(p^{c_i}\)。所以

\[f(n)=\prod\limits_{i=1}^mf(p_i^{c_i})=\prod\limits_{i=1}^mp_i^{c_i}=n \]

拉格朗日插值

多项式定义

\(n\) 次多项式定义为：

\[f(x) = a_nx^n + a_{n-1}x^{n-1} + a_{n-2}x^{n-2} + \cdots + a_0 \]

问题引入 1

给出 \(n + 1\) 个点 \((x_i, y_i)\)，表示 \(f(x_i)=y_i\)，要求求出 \(f(t)\)。

做法 1：

高斯消元。

我们把这 \(n + 1\) 个点都写出来：

\[\left\{ \begin{aligned} a_nx_1^n + a_{n-1}x_1^{n-1} + a_{n-2}x_1^{n-2}+\cdots+a_0 = y_1 ~~~~ \\ a_nx_2^n + a_{n-1}x_2^{n-1} + a_{n-2}x_2^{n-2}+\cdots+a_0 = y_2 ~~~~ \\ \vdots ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \\ a_nx_{n+1}^n + a_{n-1}x_{n+1}^{n-1} + a_{n-2}x_{n+1}^{n-2}+\cdots+a_0 = y_{n+1} \\ \end{aligned} \right. \]

把 \(a_n, a_{n-1}, \cdots, a_0\) 看作未知量，把 \(x,y\) 看作已知量，这样就可以高斯消元进行求解，解出 \(a_n,a_{n-1},\cdots,a_0\)。

时间复杂度：\(O(n^3)\)。

做法 2：

拉格朗日插值。

我们从二次多次项入手，我们给出 \(3\) 个点，\((x_1, y_1), (x_2, y_2), (x_3,y_3)\)。

可以列出方程组：

\[\left\{ \begin{aligned} a_2x_1^2 + a_1x_1^1+a_0 = y_1 \\ a_2x_2^2 + a_1x_2^1+a_0 = y_2 \\ a_2x_3^2 + a_1x_3^1+a_0 = y_3 \\ \end{aligned} \right. \]

这里我们换一种思路，我们构造三个函数，然后相加得到我们要的函数。

构造函数：

\[f_1(x)=\left\{ \begin{aligned} 1 && (x=x_1) \\ 0 && (x\not=x_1) \\ \end{aligned} \right. , f_2(x)=\left\{ \begin{aligned} 1 && (x=x_2) \\ 0 && (x\not=x_2) \\ \end{aligned} \right. ,f_3(x)=\left\{ \begin{aligned} 1 && (x=x_3) \\ 0 && (x\not=x_3) \\ \end{aligned} \right. \]

那么当 \(f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)\)，此时就满足 \(x\) 的取值为 \(x_1,x_2,x_3\) 时，\(f(x)\) 的取值就为 \(y_1,y_2,y_3\)。

现在考虑 \(n\) 次多项式，给出 \(n + 1\) 个点：\((x_1, y_1), (x_2, y_2),\cdots,(x_{n + 1}, y_{n + 1})\)。

我们考虑构造 \(f_1 \sim f_{n+1}\)。

设

\[f_i(k)=\prod\limits_{j\not=i}\dfrac{k-x_j}{x_i-x_j} \]

则 \(f_i(k)\) 满足在 \(k=x_i\) 时函数值为 \(1\)，在 \(k=x_j(j \not= i)\) 时函数值为 \(0\)。

那么多项式

\[f(k)=\sum\limits_{i=1}^{n+1}y_if_i(k)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\not=i}\dfrac{k-x_j}{x_i-x_j} \]

就满足当 \(k\) 为 \(x_1, x_2, \cdots, x_{n+1}\) 时 \(f(k)\) 的取值就为 \(y_1, y_2, \cdots, y_{n + 1}\)。

如果我们只要求解特定函数值，我们可以使用逆元来代替除法。

时间复杂度：\(\mathcal O(n^2)\)。

优化 1

当 \(x_a\) 与 \(a\) 有关系时，例如 \(x_a=a\) 时，可以使用前缀乘来优化，时间复杂度：\(\mathcal O(n)\)。

优化 2

重心拉格朗日插值法

当我们要解决两个操作：

1. 往点集里加一个点 \((x, y)\)
2. 给出一个 \(k\)，求出当前点集对应的多项式 \(f(k)\) 的值。

考虑把原式子拿去转换。

\[f(k)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\not=i}\dfrac{k-x_j}{x_i-x_j} = \sum\limits_{i=1}^{n + 1} \dfrac{y_i \times \prod\limits_{j\not=i} (k-x_j)}{\prod\limits_{j\not=i} (x_i-x_j)}=\sum\limits_{i=1}^{n+1}\dfrac{y_i \prod\limits_{j=1}^{n + 1}(k-x_j)}{(k-x_i) \prod\limits_{j \not=i}(x_i-x_j)} \]

记

\[g(k) = \prod\limits_{i=1}^{n + 1}(k-x_i) \]

则

\[f(k) = g(k)\times \sum\limits_{i=1}^{n+1}\dfrac{y_i }{(k-x_i) \prod\limits_{j \not=i}(x_i-x_j)} \]

再记

\[t(i) = \dfrac{y_i}{\prod\limits_{j \not= i} (x_i - x_j)} \]

则

\[f(k)=g(k) \times \sum\limits_{i=1}^{n+1}\dfrac{t(i)}{k-x_i} \]

新加点时 \(\mathcal O(n \log n)\) 求出 \(t(i)\)，故每加一个点就能做到 \(\mathcal O(n \log n)\) 解决。

一些常识

1. \(\sum\limits_{i=1}^n i^k\) 是一个 \(k + 1\) 次的多项式，可以插值。