OI笔记 | 基础数学知识与模板
组合数学
可爱的组合数学。
组合数的计算方式
设模数为 $P$，下列方式中如果可能要求 $P$ 是素数，则带 $\to$ 号。
加法公式 $O(n^2)$：
\[\binom{n}{m} = \binom{n - 1}{m} + \binom{n - 1}{m - 1}\]拆素因子计算 $O(n + \frac{n}{\ln n}\log n)$：
时间复杂度比较粗略。具体做法是枚举所有 $n$ 以内的素数 $p_i$，设 $x!$ 中有 $\operatorname{cal}(x)$ 个素因子是 $p_i$，有公式：
\[\operatorname{cal(x)} = \sum\limits_{i=1}^{\infty}\lfloor \frac{x}{p^i} \rfloor\]则答案为
\[\sum_{i} p_i^{\operatorname{cal}(n)-\operatorname{cal}(m)-\operatorname{cal}(n-m)} \pmod P\]$\to$ 定义式 $O(n)$：
\[\binom{n}{m}=\frac{n!}{m!(n-m)!}=\frac{n^{\underline{m}}}{m!}\]$\to$ Lucas 定理 $O(p+\log m)$：
\[\binom{n}{m}\bmod P = \binom{\left\lfloor n/p \right\rfloor}{\left\lfloor m/p\right\rfloor}\cdot\binom{n\bmod p}{m\bmod p}\bmod P\]$\to$ 线性递推 $O(m)$：
\[\binom{n}{k} = \frac{n-k+1}{k}\cdot \binom{n}{k-1}, \binom{n}{0}=1\]组合恒等式
二项式定理：
\[(a+b)^n = \sum\limits_{i=0}^n \binom{n}{i} a^i b^{n-i}\]推论：
\[a=b=1\Rightarrow 2^n=\sum\limits_{i=0}^{n} \binom{n}{i}\] \[a=1,b=-1\Rightarrow 0=\sum\limits_{i=0}^{n} \binom{n}{i}(-1)^i\]列求和（考虑杨辉三角）：
\[\sum\limits_{i=n}^m \binom{i}{n}= \binom{m+1}{n+1}\]吸收恒等式：
\[\binom{n}{m}=\frac{n}{m} \binom{n-1}{m-1}\]范德蒙德卷积：
\[\sum\limits_{i=0}^k \binom{n}{i}\binom{m}{k-i} = \binom{n+m}{k}\]组合数应用
格路计数
从 $(0,0)$ 走到 $(n,m)$。
无限制：
\[\binom{n+m}{n}\]若不能碰到某条直线，例如 $y=x+1$，则考虑反射，我们假设一条线路碰到了这条直线，到达终点 $(n,m)$ 对 $y=x+1$ 的对称点 $(m-1,n+1)$，则把这些线路减掉：
\[\binom{n+m}{m}-\binom{n+m}{m - 1}\]当 $n=m$ 时上式是卡特兰数。
错排数
\[D_n = (n-1)(D_{n-2} + D_{n-1})\]Lucas定理
求 $\binom{n+m}{m}\bmod p$。
题解
Lucas 定理的内容：
\[\binom{n}{m}\bmod p = \binom{\left\lfloor n/p \right\rfloor}{\left\lfloor m/p\right\rfloor}\cdot\binom{n\bmod p}{m\bmod p}\bmod p\]其中 $n\bmod p$ 和 $m\bmod p$ 都小于 $p$，直接暴力计算即可。而 $\binom{\left\lfloor n/p \right\rfloor}{\left\lfloor m/p\right\rfloor}$ 可以递归计算。这样我们预处理只需要 $O(p)$ 处理阶乘，然后 $O(\log^2 n)$ 计算。总体可以认为接近 $O(p)$。
注意算组合数时要判 $n<m$ 时 $\binom{n}{m}=0$。
参考代码：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 1e5 + 10;
ll fac[MAX_N], p;
inline ll read() {...}
inline void write(ll x) {...}
void exgcd(ll a, ll b, ll &x, ll &y) {
if(!b) {x = 1, y = 0; return;}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
inline ll inv(ll v) {
ll x, y;
exgcd(v, p, x, y);
return (x % p + p) % p;
}
inline ll C(ll n, ll m) {
if(n < m) return 0;
return fac[n] * inv(fac[m]) % p * inv(fac[n - m]) % p;
}
ll Lucas(ll n, ll m) {
if(!m) return 1;
return C(n % p, m % p) * Lucas(n / p, m / p) % p;
}
int main() {
int T = read();
fac[0] = fac[1] = 1;
while(T--) {
ll n = read(), m = read();
p = read();
for(int i = 2; i <= p; i++) fac[i] = (fac[i - 1] * i) % p;
write(Lucas(n + m, m));
putchar('\n');
}
return 0;
}
数论
可爱的数论。
整除分块
利用整除分块，可以计算类似于这样的式子：
\[\sum\limits_{i=1}^n f(\lfloor \frac{n}{i} \rfloor)\]模板题
求 $\sum\limits_{i=1}^n \lfloor \frac{n}{i} \rfloor$。
核心就是 $\lfloor \frac{n}{i} \rfloor$ 的取值只有 $O(\sqrt n)$ 种。
我们考虑 $\lfloor \frac{n}{i} \rfloor$ 的值在某一段 $i\in [l,r]$ 内的值是相同的。例如当 $i\in [16,30], n = 30$ 时，恒有 $\lfloor \frac{n}{i} \rfloor = 2$。下面的代码可以找到跳过这样的区间：
for(int i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ans += (n / i) * (j - i + 1);
}
例题1
设 $f(x)$ 表示 $x$ 所有约数的和，求 $\sum\limits_{i=l}^r f(i)$。
考虑贡献，即 $i$ 在 $[1,n]$ 上能作为 $\lfloor \frac{n}{i} \rfloor$ 个数的约数。我们设 $g(x)=\sum\limits_{i=1}^n i\cdot \lfloor \frac{n}{i} \rfloor$，则答案即为 $g(r)-g(l-1)$。
和模板相比，其实就是和式中乘个 $i$。只要把 (j-i+1) 改为等差数列求和公式即可。
ll f(ll n) {
ll ret = 0;
for(ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ret += (n / i) * ((i + j) * (j - i + 1) / 2);
}
return ret;
}
例题2
给定 $n, k$，求 $\sum\limits_{i = 1}^n k \bmod i$。
由于 $M\bmod N = M - N \lfloor \frac{M}{N}\rfloor $，原式可以改写成：
\[\sum\limits_{i = 1}^n k \bmod i=\sum\limits_{i = 1}^n k-i \lfloor \frac{k}{i}\rfloor=nk-\sum\limits_{i=1}^n i\lfloor \frac{k}{i} \rfloor\]套用 例题1 即可。注意细节处理。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
inline ll read() {...}
inline void write(ll x) {...}
int main() {
ll n = read(), k = read();
ll ans = 0;
for(ll i = 1, j; i <= n; i = j + 1) {
j = min(n, (k / i == 0) ? n : k / (k / i));
ans += (i + j) * (j - i + 1) / 2 * (k / i);
}
write(n * k - ans);
return 0;
}
线性筛
其实一般埃氏筛法不会被卡，不过这题要写欧拉筛。
#include <bits/stdc++.h>
using namespace std;
const int MAX_N = 1e8+3;
vector<int> a;
bool vis[MAX_N];
inline int read() {...}
void write(int x) {...}
int main() {
int n = read(), q = read();
for(int i = 2; i <= n; i++) {
if(!vis[i]) a.push_back(i);
for(int k = 0; k < a.size() && a[k] * i <= n; k++) {
vis[a[k] * i] = 1;
if(i % a[k] == 0) break;
}
}
for(int i = 0; i < q; i++) {
int k = read();
write(a[k - 1]);
putchar('\n');
}
return 0;
}
可以顺带筛任何积性函数。给出一个筛欧拉函数的模板：
void pre() {
phi[1] = 1;
for(int i = 2; i < MAX_N; i++) {
if(!p[i]) primes[++cnt] = i, phi[i] = i - 1;
for(int j = 1; j <= cnt && i * primes[j] < MAX_N; j++) {
p[i * primes[j]] = 1;
if(i % primes[j] == 0) {
phi[i * primes[j]] = phi[i] * primes[j];
break;
}
else phi[i * primes[j]] = phi[i] * (primes[j] - 1);
}
}
}
中国剩余定理(CRT)
模板题：洛谷 P4777
给出下面的同余方程组：
\[\begin{cases} x &\equiv b_1 \pmod {a_1} \\ x &\equiv b_2 \pmod {a_2} \\ &\vdots \\ x &\equiv b_k \pmod {a_k} \\ \end{cases}\]若 $a_1,a_2,\dots a_n$ 两两互质，则可以用中国剩余定理求解这个方程组。流程如下：
-
求出所有模数的乘积 $prod=\prod\limits_{i=1}^k a_i。$
-
记 $m_i=\frac{prod}{a_i}$，求出它在模 $a_i$ 意义下的逆元 $m_i^{-1}$。
-
答案即为 $\sum\limits_{i=1}^k b_i m_i m_i^{-1} \pmod {prod}$。
按照上面的流程就可以写出模板题的代码。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 1e5 + 10;
ll a[MAX_N], b[MAX_N], prod = 1, n, ans;
inline ll read() {...}
inline void write(ll x) {...}
void exgcd(ll a, ll b, ll &x, ll &y) {
if(!b) {x = 1, y = 0; return;}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
ll inv(ll v, ll p) {
ll x, y;
exgcd(v, p, x, y);
return (x % p + p) % p;
}
int main() {
n = read();
for(int i = 1; i <= n; i++) {
a[i] = read(), b[i] = read();
prod *= a[i];
}
for(int i = 1; i <= n; i++) {
ll m = prod / a[i];
(ans += b[i] * m * inv(m, a[i])) %= prod;
}
write(ans);
return 0;
}
[例题] 古代猪文
给出 $n, G$，求：
\[G^{\sum_{k\mid n}\binom{n}{k}} \pmod{999911659}\]题解
考虑扩展欧拉定理：$a^{b}\equiv a^{b\bmod \varphi(m)} \pmod m$，由于 $999911659$ 为质数，则 $\varphi(999911659)=999911658$，所以原式相当于求：
\[G^{\sum_{k\mid n}\binom{n}{k}\bmod 999911658} \pmod{999911659}\]其实就是要求 $\sum_{k\mid n}\binom{n}{k}\bmod 999911658$。如果直接求解，无法预处理出范围内的所有阶乘，也无法保证 $[1,999911657]$ 中的每个数都存在逆元。
考虑我们可以用 $Lucas$ 定理计算模数较小时的组合数取模，我们将 $999911658$ 分解质因数为 $2 \times 3 \times 4679 \times 35617$，分别计算出组合数在模这四个小质数下的答案 $b_1,b_2,b_3,b_4$，然后再用中国剩余定理求解以下同余方程组：
\[\begin{cases} x \equiv b_1 \pmod 2\\ x \equiv b_2 \pmod 3\\ x \equiv b_3 \pmod {4679}\\ x \equiv b_4 \pmod {35617} \end{cases}\]则方程组的解就是答案的指数，用快速幂算出答案即可。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 999911658;
ll a[5] = {0, 2, 3, 4679, 35617}, b[5], fac[40000];
inline ll read() {...}
inline void write(ll x) {...}
ll qpow(ll a, ll b, ll p) {
ll ret = 1;
while(b) {
if(b & 1) (ret *= a) %= p;
(a *= a) %= p;
b >>= 1;
}
return ret;
}
ll inv(int v, int p) {return qpow(v, p - 2, p);}
void init_fac(ll p) {
fac[0] = 1;
for(ll i = 1; i <= p; i++)
fac[i] = (fac[i - 1] * i) % p;
}
ll C(ll n, ll m, ll p) {
if(n < m) return 0;
return fac[n] * inv(fac[n - m], p) % p * inv(fac[m], p) % p;
}
ll Lucas(ll n, ll m, ll p) {
if(!m) return 1;
return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
}
ll CRT() {
ll ans = 0;
for(int i = 1; i <= 4; i++) {
ll m = MOD / a[i];
(ans += b[i] * m % MOD * inv(m, a[i]) % MOD) %= MOD;
}
return ans;
}
int main() {
ll n = read(), G = read();
if(G == MOD + 1) {
write(0);
return 0;
}
for(ll t = 1; t <= 4; t++) {
ll p = a[t];
init_fac(p);
for(ll i = 1; i * i <= n; i++) {
if(n % i == 0) {
(b[t] += Lucas(n, i, p)) %= p;
if(i * i != n)
(b[t] += Lucas(n, n / i, p)) %= p;
}
}
}
write(qpow(G, CRT(), MOD + 1));
return 0;
}
多项式
可爱的多项式。
拉格朗日插值
请你用给定的 $n$ 个点确定一个多项式，并求出 $f(k) \bmod 998244353$ 的值。
$1 \le n \leq 2\times 10^3$
题解
可以构造拉格朗日差值的基本形式：
\[f(x)=\sum\limits_{i=1}^{n}y_i \prod\limits_{j\neq i}\frac{x-x_j}{x_i-x_j}\]容易验证这是正确的。把任意一点 $(x_i,y_i)$ 代入都满足式子，因为除了第 $i$ 项其他的分子都为 $0$。
构造的方法是构造线性同余方程组然后用 CRT 求解，这里不详细说了。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 2e3 + 10, MOD = 998244353;
ll x[MAX_N], y[MAX_N], ans, a, b, n, k;
inline ll read() {...}
inline void write(ll x) {...}
void exgcd(ll a, ll b, ll &x, ll &y) {
if(!b) {x = 1, y = 0; return;}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
ll qpow(ll x, ll y) {
ll ret = 1;
while(y) {
if(y & 1) (ret *= x) %= MOD;
(x *= x) %= MOD;
y >>= 1;
}
return ret;
}
ll inv(int v) {return qpow(v, MOD - 2);}
int main() {
n = read(), k = read();
for(int i = 1; i <= n; i++)
x[i] = read(), y[i] = read();
for(int i = 1; i <= n; i++) {
a = y[i] % MOD, b = 1;
for(int j = 1; j <= n; j++)
if(i != j) (a *= k - x[j]) %= MOD, (b *= x[i] - x[j]) %= MOD;
(ans += a * inv(b) % MOD) %= MOD;
}
write((ans % MOD + MOD) % MOD);
return 0;
}
[例题] The Sum of the k-th Powers
可以证明 $\sum_{i=1}^n i^k$ 可以表示为关于 $n$ 的 $k + 1$ 次多项式。
例如：
若 $k = 1$，则 $f(n)=\frac{n(n + 1)}{2}$；
若 $k=2$，则 $f(n)=\frac{n(n+1)(2n+1)}{2}$；
若 $k=3$，则 $f(n) = \frac{n^2(n+1)^2}{4}$；
$\dots$
那么我们取 $k+2$ 个点进行拉格朗日插值即可得出这个多项式。
由于我们可以任意取点，考虑如果我们取 $(1,f(1)),(2,f(2)),\dots (k + 2,f(k+2))$ 这样的连续整数，则可以把拉格朗日插值从 $O(n^2)$ 优化到 $O(n)$。此时基本形式转化为：
\[f(x)=\sum\limits_{i=1}^{k+2}y_i \prod\limits_{j\neq i}\frac{x-j}{i-j}\]在本题中，$y_i=\sum_{j=1}^i j^k$，要求的值就是 $f(n)$。
考虑分子，处理出 $(n-i)$ 的前后缀积即可。
考虑分母，它等于 $(i-1)! \times (k + 2 - i)! \times (-1)^{k+2-i}$，预处理出阶乘逆元即可。
时间复杂度： $O(k\log k)$。$\log$ 是快速幂带的。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_K = 1e6 + 10, MOD = 1e9 + 7;
int prodl[MAX_K], prodr[MAX_K], fac[MAX_K], inv[MAX_K];
inline ll read() {...}
inline void write(ll x) {...}
ll qpow(ll x, ll y) {
ll ret = 1;
for(; y; y >>= 1, (x *= x) %= MOD)
if(y & 1) (ret *= x) %= MOD;
return ret;
}
int main() {
int n = read(), k = read(), y = 0, ans = 0, a, b;
prodl[0] = prodr[k + 3] = fac[0] = 1;
for(int i = 1; i <= k + 3; i++) {
prodl[i] = (1ll * prodl[i - 1] * (n - i)) % MOD;
fac[i] = (1ll * fac[i - 1] * i) % MOD;
}
inv[k + 3] = qpow(fac[k + 3], MOD - 2);
for(int i = k + 2; i >= 0; i--) {
prodr[i] = (1ll * prodr[i + 1] * (n - i)) % MOD;
inv[i] = (1ll * inv[i + 1] * (i + 1)) % MOD;
}
for(int i = 1; i <= k + 2; i++) {
y = (1ll * y + qpow(i, k)) % MOD;
a = 1ll * prodl[i - 1] * prodr[i + 1] % MOD;
b = ((k + 2 - i) & 1 ? -1ll : 1ll) * inv[i - 1] * inv[k + 2 - i] % MOD;
ans = (ans + 1ll * y * a % MOD * b % MOD) % MOD;
}
write((ans % MOD + MOD) % MOD);
return 0;
}
[例题] calc
小数据可以 dp。设 $dp_{i,j}$ 为 $i$ 个数的序列，值域为 $[1,j]$ 时的权值和，则有：
\[dp_{i,j} = j \cdot dp_{i-1,j-1} + dp_{i,j-1}\]然后可以证明这个 $dp$ 的通项是 $2n+1$ 次多项式，拉格朗日插值即可。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 2000;
ll n, k, MOD, facn, ans, dp[MAX_N][MAX_N];
inline ll read() {...}
inline void write(ll x) {...}
ll qpow(ll x, ll y) {
ll ret = 1;
for(; y; y >>= 1, (x *= x) %= MOD)
if(y & 1) (ret *= x) %= MOD;
return ret;
}
int main() {
k = read(), n = read(), MOD = read();
facn = 1;
for(int i = 1; i <= n; i++)
(facn *= i) %= MOD;
for(int i = 0; i <= k && i <= (2 * n + 1); i++)
dp[0][i] = 1;
for(int i = 1; i <= n; i++)
for(int j = 1; j <= k && j <= (2 * n + 1); j++)
dp[i][j] = (j * dp[i - 1][j - 1] + dp[i][j - 1]) % MOD;
if(k <= 2 * n + 1) ans = facn * dp[n][k] % MOD;
else {
for(int i = 1; i <= 2 * n + 1; i++) {
ll a = dp[n][i], b = 1;
for(int j = 1; j <= 2 * n + 1; j++) {
if(i != j) {
(a *= k - j + MOD) %= MOD;
(b *= i - j + MOD) %= MOD;
}
}
(ans += a * qpow(b, MOD - 2) % MOD) %= MOD;
}
(ans *= facn) %= MOD;
}
write(ans);
return 0;
}
FFT
struct cp {
double x, y;
cp(double _x = 0, double _y = 0) : x(_x), y(_y) {}
cp operator + (const cp &t) const {return cp(x + t.x, y + t.y);}
cp operator - (const cp &t) const {return cp(x - t.x, y - t.y);}
cp operator * (const cp &t) const {return cp(x * t.x - y * t.y, x * t.y + y * t.x);}
};
void pre(){
while(len <= n + m)
len <<= 1, k++;
for(int i = 0; i <= len; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
void FFT(cp *P, int type) {
pre();
for(int i = 0; i < len; i++)
if(i < rev[i]) swap(P[i], P[rev[i]]);
for(int i = 1; i < len; i <<= 1) {
cp omega(cos(PI / i), type * sin(PI / i));
for(int k = 0; k < len; k += (i << 1)) {
cp now(1, 0);
for(int l = 0; l < i; l++) {
cp u = P[k + l], v = now * P[k + i + l];
P[k + l] = u + v;
P[k + i + l] = u - v;
now = now * omega;
}
}
}
}
线性代数
可爱的线代。
高斯消元
给定一个线性方程组，对其求解。
题解
高斯-约旦消元法。
我们依次枚举每一列，找到这一列中系数绝对值最大的一行，把这一行与当前行交换。然后对其他行加减消元。最终转化为阶梯矩阵。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAX_N = 105;
const double eps = 1e-5;
int n;
double a[MAX_N][MAX_N];
inline ll read() {...}
inline void write(ll x) {...}
bool Gauss() {
for(int i = 1; i <= n; i++) {
int mx = i;
for(int j = i + 1; j <= n; j++)
if(fabs(a[j][i]) > fabs(a[mx][i])) mx = j;
for(int j = 1; j <= n + 1; j++)
swap(a[i][j], a[mx][j]);
if(fabs(a[i][i]) <= eps) return 0;
for(int j = 1; j <= n; j++) {
if(j == i) continue;
double d = a[j][i] / a[i][i];
for(int k = i + 1; k <= n + 1; k++)
a[j][k] -= a[i][k] * d;
}
}
return 1;
}
int main() {
n = read();
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n + 1; j++)
a[i][j] = read();
if(Gauss()) {
for(int i = 1; i <= n; i++)
printf("%.2lf\n", a[i][n + 1] / a[i][i]);
}
else puts("No Solution");
return 0;
}
♥