前置知识

历史资料

一些以前的资料，扫描备份（然后本体就可以扔了，芜湖～）。

目录
整除 1
最大公约数和最小公倍数 5
素数及唯一分解定理 11
同余 15
母函数 19
递推数列 22

逆元

若 $ax \equiv 1 \pmod{b}$ ，则 $x$ 为 $a$ $\mathrm{mod}$ $b$ 的逆元，记为 $a^{-1}$ .

逆元概念的理解：

逆元约等于模 b 意义下的倒数，这个概念和离散数学中的差不多。
在一些算法竞赛中，由于各种缘故，并不是很想处理浮点数，而是用整数替代。遇到两数相除的情况，就不会算出小数，而是通过“逆元”的概念，得到一个整数。

欧几里得算法

这个算法用来求两个数的最大公约数，它基于以下定理：

$\mathrm{gcd}(a,b) = \mathrm{gcd}(b, a \ \mathrm{mod} \ b)$

这个定理可以通过几何直观理解。

由此写出代码：

int gcd(int a, int b)
{
    return b ? gcd(b, a % b) : a;
}
// 或者使用函数 __gcd();

线性同余方程

形如 $ax \equiv b \pmod{n}$ 的方程称为线性同余方程（Linear Congruence Equation），其中 a, b, n 为给定整数，$x$ 为未知数，需要从 $[0, n-1]$ 中求解 $x$，当解不唯一时需要求出全体解。

用逆元求解

ref: https://oi-wiki.org/math/number-theory/linear-equation/

首先考虑简单的情况，当 $a$ 和 $n$ 互素（coprime 或 relatively prime）时，即 $\gcd(a, n) = 1$。

此时可以计算 $a$ 的逆元，并将方程的两边乘以 $a$ 的逆元，可以得到唯一解。

$x\equiv ba ^ {- 1} \pmod n$

接下来考虑 $a$ 和 $n$ 不互素（not coprime），即 $\gcd(a, n) \ne 1$ 的情况。此时不一定有解。例如，$2x\equiv 1\pmod 4$ 没有解。

设 $g = \gcd(a, n)$，即 $a$ 和 $n$ 的最大公约数，其中 $a$ 和 $n$ 在本例中大于 1。

当 $b$ 不能被 $g$ 整除时无解。此时，对于任意的 $x$，方程 $ax\equiv b\pmod n$ 的左侧始终可被 $g$ 整除，而右侧不可被 $g$ 整除，因此无解。

如果 $g$ 整除 $b$，则通过将方程两边 $a$、$b$ 和 $n$ 除以 $g$，得到一个新的方程：

$a^{'}x\equiv b^{'} \pmod{n^{'}}$

其中 $a^{‘}$ 和 $n^{‘}$ 已经互素，这种情形已经解决，于是得到 $x^{‘}$ 作为 $x$ 的解。

很明显，$x^{‘}$ 也将是原始方程的解。这不是唯一的解。可以看出，原始方程有如下 $g$ 个解：

$x_i\equiv (x^{'} + i\cdot n^{'}) \pmod n \quad \text{for } i = 0 \ldots g-1$

总之，线性同余方程的 解的数量 等于 $g = \gcd(a, n)$ 或等于 $0$。

一点补充说明：

同余式可以三个数同时除掉 gcd. 即若 $d \ | \ a,b,m$ , 且 $a \equiv b \pmod{m}$，则：

$\frac{a}{d} \equiv \frac{b}{d} \pmod{\frac{m}{d}}$

证明比较简单，略。

扩展欧几里得（exgcd）

ref: https://blog.csdn.net/weixin_43872728/article/details/107289833

一、线性同余方程 $ax\equiv b \pmod n$ 可以改写为如下线性不定方程：

$ax + nk = b$

其中 $x$ 和 $k$ 是未知数。这两个方程是等价的，有整数解的充要条件为 $\gcd(a,n) \mid b$。

二、由裴属定理，对于整数 a、b，存在 a、b 的线性组合等于 $\mathrm{gcd} (a,b)$ .

换句话说，如果 ax+by=m 有解，那么 m 一定是 gcd(a,b) 的若干倍（可以来判断一个这样的式子有没有解）。

进一步提出疑问：在有解的情况下，这个解是多少，也就是 x 和 y 的值是多少？

扩展欧几里得模板：

// 求x, y，使得ax + by = gcd(a, b)
int exgcd(int a, int b, int &x, int &y)
{
    if (!b)
    {
        x = 1; y = 0;
        return a;   //到达递归边界开始向上一层返回
    }
    int d = exgcd(b, a % b, x, y);
    int temp = y;    //推出这一层的x，y
    y = x - (a/b)*y;
    x = temp;
    return d;
}

接下来解释一下上面代码。

编写递归代码的核心是将还没有发生的事当作已经发生，然后严格从定义出发描述行为，规定边界情形。

我们考虑相邻两层之间的状态关系，由回溯过来的 x,y 得到新的 x,y .

假设栈中当前层得到的解是 x1，y1；如何求下一层的状态？代入方程得：

$b x_1 + (a \ \mathrm{mod} \ b)y_1 = \mathrm{gcd}(a,b)$ $bx_1 + (a-(a/b)b)y_1 = \mathrm{gcd}(a,b)$ $ay_1 + b(x_1-(a/b)y_1) = \mathrm{gcd}(a,b)$

于是，由上一层的解 x1,y1，可以推出下一层 x,y 的状态：

$x = y_1 \qquad y=x_1-(a/b)y_1$

接下来考虑边界条件。

当递归到达边界时，余数 b==0 ，除数 a 就是最大公约数，此时可得出方程的一组解：

$x = 1, \ y = 0$ $a \cdot 1 + b \cdot 0 = \mathrm{gcd}(a,b)$

由上述模板可以求出方程 $ax+by=\mathrm{gcd}(a,b)$ 的任意一组解 $x_0$，$y_0$ ，那么方程 $ax+by=k$ 的通解的形式是什么样的呢？注意， $ax+by=k$ 中的 k 一定是 gcd(a,b) 的倍数，方程才会有解。

因为由 ax+by=gcd(a,b) 变成 ax+by=k 方程扩大了 k/gcd(a,b)倍，所以通过 exgcd 求出来的 x0，y0 也扩大了相同的倍数。

$x_0 = x_0 \cdot \frac{k}{\mathrm{gcd}(a,b)} \qquad y_0 = y_0 \cdot \frac{k}{\mathrm{gcd}(a,b)}$

通解形式为：

$x = x_0 + \frac{b}{\mathrm{gcd}(a,b)} \cdot n$ $y = y_0 + \frac{a}{\mathrm{gcd}(a,b)} \cdot n$

相当于 x 每次可以增减 b/gcd 的整数倍，y 每次可以增减 a/gcd 的整数倍。注意：x 求出来后，y 通常由 x 代入方程求得。

最小正整数解：

$x = (x + \frac{b}{\mathrm{gcd}(a,b)} \cdot n) \ \mathrm{mod} \ \frac{b}{\mathrm{gcd}(a,b)} = x \ \mathrm{mod} \ \frac{b}{\mathrm{gcd}(a,b)}$

其中 b/gcd(a,b) 应取正。若 x<=0，则 x += b/gcd .

我们发现，exgcd 虽然是用来求线性同余方程的，但只需稍作改动，将其中一个数改为 1，就满足逆元的定义。因此 exgcd 也可以用来求逆元。

快速幂法求逆元

刚才介绍了扩展欧几里得法求逆元，其实也可以快速幂法求逆元。

逆元定义式（这里设 b 为素数）：

$ax \equiv 1 \pmod{b}$

这里显然 a,b 是互素的。

由费马小定理得：

$a^b \equiv a \pmod{b}$ $a^{b-1} \equiv 1 \pmod{b}$

所以有：

$a^bx \equiv ax \equiv 1 \equiv a^{b-1} \pmod{b}$

即：

$ax \equiv a^{b-1} \pmod{b}$

所以：

$x \equiv a^{b-2} \pmod{b}$

然后我们就可以用快速幂来求了。

拉格朗日插值法

算法的原理部分，参考工程数学相关部分。

本文只探讨对于任意点直接求解的最简单情形。

题目链接： https://www.luogu.com.cn/problem/P4781

ac 代码（下标从 1 开始计数）：

// 参考：https://www.cnblogs.com/purinliang/p/13670493.html
#include <cstdio>

using ll = long long;


namespace LagrangeInterpolation {

static const int MOD = 998244353;

ll qpow (ll x, ll n) {
    ll res = 1;
    while (n) {
        if (n & 1) {
            res = res * x % MOD;
        }
        x = x * x % MOD;
        n >>= 1;
    }
    return res;
}

ll inv (ll x) {
    return qpow (x, MOD - 2);
}

/**
 * 用 n 个点 (x[1], y[1]), (x[2], y[2]), ..., (x[n], y[n])
 * 插值得到 n - 1 次多项式 L，直接求解 L(x0) 的值
 *
 * 时间复杂度: O(n^2)
 * 空间复杂度：O(1)
*/
ll lagrange_interpolation (ll x0, ll *x, ll *y, int n) {
    for (int i = 1; i <= n; ++i) {
        if (x0 == x[i]) {
            return y[i];
        }
    }
    ll Lx0 = 0;
    ll P = 1;
    for (int i = 1; i <= n; ++i) {
        P = P * (x0 - x[i] + MOD) % MOD;
    }
    for (int i = 1; i <= n; ++i) {
        ll q = 1;
        for (int j = 1; j <= n; ++j) {
            if (j == i) {
                continue;
            }
            q = q * (x[i] - x[j] + MOD) % MOD;
        }
        ll p = P * inv (x0 - x[i] + MOD) % MOD;
        ll lix0 = p * inv (q) % MOD;
        Lx0 = (Lx0 + lix0 * y[i] % MOD) % MOD;
    }
    return Lx0;
}

};

const ll MAXN = 2e3+5;
ll xx[MAXN], yy[MAXN];

int main()
{
    ll n, k;
    scanf("%lld%lld", &n, &k);
    for(int i = 1; i <= n; i++){
        ll x, y;
        scanf("%lld%lld", &x, &y);
        xx[i] = x;
        yy[i] = y;
    }

    printf("%lld", LagrangeInterpolation::lagrange_interpolation(k, xx, yy, n));
    return 0;
}