板子
目录
1.数学
1.1 高精度计算
我只会java，这里贴一个java常用板子。
设此时有 BigInteger a,b，那么读入函数是 sc.nextBigInteger()。
int 类型的读入函数是 sc.nextInt()。
a+b：a.add(b)
a-b：a.substract(b)
a*b：a.multiply(b)
a/b：a.divide(b)
a%b：a.remainder(b)
gcd(a,b)：a.gcd(b)
下面是 a-b 的高精度运算。
import java.math.*;
import java.util.*;
import java.io.*;
import java.text.*;
public class Main{
public static void main(String[] args){
Scanner sc = new Scanner(System.in);
BigInteger a,b;
a = sc.nextBigInteger(); b = sc.nextBigInteger();
System.out.println(a.subtract(b));
}
}
windows 下 cmd 先 javac Main.java 再 java Main 。
1.2 数论
1.2.1 素数
有素数计数函数 \(\pi(x) \sim \frac{x}{\ln x}\)
一个数的质因数约有 \(\frac{\log n}{\log \log n}\) 个
1.2.1.1 判断单个数是否为质数
时间 \(O(\sqrt n)\)，空间 \(O(1)\)
bool isprime(int x) {
if(x == 1) return false;
int m = sqrt(x + 0.5);
for(int i=2;i<=m;i++) {
if(x % i == 0) return false;
}
return true;
}
1.2.1.2 欧拉筛：质数、欧拉函数、莫比乌斯函数、约数个数、约数个数和
时间 \(O(n)\) ，空间 \(O(n)\)
理论上只要是积性函数就可以使用欧拉筛筛出来。
#include<cstdio>
using namespace std;
#define IL inline
typedef long long LL;
const int N = 1e6 + 3;
bool np[N];
int pri[N],phi[N],mu[N],di[N],si[N];
IL int Euler_sieve(int n) {
int tot = 0;
phi[0] = mu[0] = di[0] = si[0] = 0;
np[1] = phi[1] = mu[1] = di[1] = si[1] = 1;
for(int i=2;i<=n;i++) {
if(!np[i]) {
pri[++tot] = i;
phi[i] = i-1;
mu[i] = -1;
di[i] = 2;
si[i] = i + 1;
}
for(int j=1;j<=tot&&i*pri[j]<=n;j++) {
np[i*pri[j]] = true;
if(i % pri[j] == 0) {
phi[i*pri[j]] = pri[j] * phi[i];
mu[i*pri[j]] = 0;
di[i*pri[j]] = (di[i] << 1) - di[i/pri[j]];
si[i*pri[j]] = si[i] + pri[j] * (si[i]-si[i/pri[j]]);
break;
}
phi[i*pri[j]] = phi[i] * (pri[j]-1);
mu[i*pri[j]] = -mu[i];
di[i*pri[j]] = di[i] << 1;
si[i*pri[j]] = si[i] * (1+pri[j]);
}
}
return tot;
}
变量名称解释：
- tot:质数个数，函数Euler_sieve(int n)也会返回质数个数。
- pri[i]:意为prime,这个数组里存储小于n的质数。
- np[i]:意为not prime。这个bool数组是判断i是否为质数的。
- phi[i]:即 \(\varphi_i\)。即为\(i\)的欧拉函数
- mu[i]:即 \(\mu_i\)。即为\(i\)的莫比乌斯函数
- di[i]:意为divisor。即为\(i\)的约数个数函数
- si[i]:即 \(\sigma_i\)。即为\(i\)的约数和函数
输入：n，意义见输出。
输出：小于等于n的所有质数以及质数个数。小于等于n的所有自然数的欧拉函数、莫比乌斯函数、约数个数函数、约数和函数。
1.2.1.3 Miller-Rabin 素性测试
时间复杂度 \(O(k \log^3 n)\) ，空间复杂度 \(O(1)\)
long long 范围内都可以测出来。由于快速幂中使用 x 作为模数，且 long long 的极限为 9e18，若 x 大于 3e9，为避免 TLE 则开 long long，为避免 WA 则开 __int128。
bool isprime(LL x) //判断 x 是不是质数。
LL gcd(LL a,LL b) { return b == 0 ? a : gcd(b,a%b);}
IL LL ksm(LL a,LL b,LL mod) {
LL res = 1LL;
while(b) {
if(b&1LL) res = (__int128)res * a % mod;
a = (__int128)a * a % mod;
b >>= 1LL;
}
return res;
}
IL bool mr(LL x,LL b) {
LL d=x-1, p=0;
while((d&1)==0) d>>=1,++p;
int i;
LL cur=ksm(b,d,x);
if(cur==1 || cur==x-1) return true;
for(i=1;i<=p;++i) {
cur=(__int128)cur*cur%x;
if(cur==x-1) return true;
}
return false;
}
IL bool isprime(LL x) {
if(x == 46856248255981 || x < 2) return false;
if(x==2 || x==3 || x==7 || x==61 || x==24251) return true;
return mr(x,2) && mr(x,3) && mr(x,7) && mr(x,61) && mr(x,24251);
}
1.2.1.4 Pollard-Rho算法随机找一个 x 的因数
时间复杂度 \(O(n^{\frac{1}{4}})\) ，空间复杂度 \(O(1)\)
LL PR(LL x) //输入一个 x ，返回 x 的一个因数。
IL LL f(LL x,LL c,LL mod) { return ((__int128)x*x+c) % mod;}
IL LL PR(LL x) {
LL s=0,t=0,c=1LL*rand()%(x-1) + 1;
LL val = 1LL;
for(ri goal=1;;goal<<=1,s=t,val=1LL) {
for(ri stp=1;stp<=goal;stp++) {
t = f(t,c,x);
val = (__int128)val * abs(t-s) % x;
if(stp%127 == 0) {
LL d = gcd(val,x);
if(d > 1) return d;
}
}
LL d = gcd(val,x);
if(d > 1) return d;
}
}
1.2.1.5 结合Miller-Rabin 和 Pollard-Rho 求出一个数的所有质因子
朴素做法大约是 \(O(q\frac{\sqrt n}{\ln n} + \sqrt n)\) ，需要打 \(O(\sqrt n)\) 的表。
此做法时间复杂度 \(O(k\log^3 n + n^\frac{1}{4}\frac{\log n}{\log \log n})\) ，空间复杂度 \(O(1)\)
需要使用上面的 Miller-Rabin 和 Pollard-Rho 板子。
变量解释：
- faccnt：质因数个数（包括相同的）。
- fac[N]：存质因数的数组。
输入：x
输出：x 的所有质因数。
int faccnt = 0;
LL fac[N];
IL void getfac(LL x) { // x == 1 时最好在函数外面特判一下，不要丢进函数里。
if(x == 1) return ;
if(isprime(x)) fac[++faccnt] = x;
else {
LL d = x; while(d >= x) d = PR(x);
getfac(d); getfac(x/d);
}
}
int gaofac(LL x) {faccnt = 0; getfac(x); return getfac;}
1.2.2 欧拉函数
欧拉函数（Euler's totient function），即 \(\varphi(n)\) ，表示的是小于等于 \(n\) 和 \(n\) 互质的数的个数。
比如说 \(\varphi(1) = 1\) 。
当 n 是质数的时候，显然有 \(\varphi(n) = n - 1\) 。
1.2.2.1 求欧拉函数
\(φ(n)=n\prod_{i=1}^s (1−\frac{1}{p_i})\)，此过程可以使用 1.2.1.5 算法优化。
1.2.2.2 欧拉定理与扩展欧拉定理
欧拉定理：若 \(\gcd (a,b)=1\)，则 $a^{\varphi(m)} \equiv 1 \pmod{m} $
扩展欧拉定理（欧拉降幂）：
\[a^b \equiv
\begin{cases}
a^{b \mod \varphi(p)}~~~~~~~ ,\gcd(a,b)=1 \\
a^b ~~~~~~~~~~~~~~~~~~~~~~~,\gcd(a,b) \neq 1,b < \varphi(p) \pmod{p} \\
a^{b \mod \varphi(p)+\varphi(p)} ,\gcd(a,p) \neq 1, b \geq \varphi(p)
\end{cases}
\]
1.2.3 乘法逆元
设 \(ax \equiv 1 \pmod{b}\)，求出 x。
快速幂：x = ksm(a,b-2,b)
扩展欧几里得算法：exgcd(a,b,d,x,y)
打表：求出 1~n 中所有数对于 p 的逆元。
inv[1] = 1;
for (int i = 2; i <= n; ++i) inv[i] = (long long)(p - p / i) * inv[p % i] % p;
线性求任意 n 个数的逆元：设有 n 个数，\(1 \leq a_i < p\) ，求出所有 \(a_i\) 的逆元。
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2,p);
// 当然这里也可以用 exgcd 来求逆元,视个人喜好而定.
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;
1.2.4 中国剩余定理（CRT）
如果模数很大而且不是质数的情况，可以将模数拆分，然后若模数唯一分解式中每个质因数次数都为 1 ，则可以分别求出表达式在不同质因数下的数值，再使用 CRT 求解。
1.2.5 卢卡斯定理（Lucas）
对于质数 p ,有
\[{n \choose m} \mod p = {\lfloor n/p \rfloor \choose \lfloor m/p \rfloor} \cdot {n \mod p \choose m \mod p} \mod p
\]
如果递归求解的话，p 的范围不能太大，一般在 \(10^5\) 左右。
IL LL C(LL n,LL m,LL p) {
if(n<m) return 0;
return jc[n]*ksm(jc[m],p-2,p)%p*ksm(jc[n-m],p-2,p)%p; //可以预处理阶乘的逆元
}
IL LL Lucas(LL n,LL m,LL p) {
if(n<m) return 0; if(!n) return 1; if(!m) return 1;
return Lucas(n/p,m/p,p) * C(n%p,m%p,p) % p;
}
1.2.6 莫比乌斯反演
1.2.6.1 积性函数
定义：若 \(x \bot y\) 且 \(f(xy)=f(x)f(y)\) ，则\(f(n)\)为积性函数。（\(x \bot y \iff gcd(x,y)=1\)）
性质：若\(f(x)\)和\(g(x)\)都为积性函数，则以下函数也为积性函数。
\[h(x)=f(x^p)
\]
\[h(x)=f^p(x)
\]
\[h(x)=f(x)g(x)
\]
\[h(x)=\sum_{d|x}f(d)g(\frac{x}{d})
\]
例子：
- 单位函数：\(\epsilon(n)=1[n=1]\)。
- 恒等函数：\(id_k(n)=n^k\)。\(id_1(n)\)通常简记作\(id(n)\)。
- 常数函数：\(1(n)=1\)。
- 除数函数：\(\sigma_k(n)=\sum_{d|n}d^k\)。\(\sigma_0(n)\)通常简记作\(d(n)\)或\(\tau(n)\)，这是约数个数函数。\(\sigma_1(n)\)通常简记作\(\sigma(n)\)，这是约数和函数。
- 欧拉函数：\(\varphi(n)=\sum_{i=1}^n[gcd(i,n)=1]\)
- 莫比乌斯函数：\[\mu(n)= \begin{cases} 1& \text{$n=1$}\\ 0& \text{$\exists d:d^2|n$}\\ (-1)^{\omega(n)}& \text{otherwise} \end{cases}\]其中\(\omega(n)\)表示n的不同质因子个数，是一个加性函数。
1.2.6.2 Dirichlet 卷积
定义：定义两个数论函数 \(f,g\) 的 Dirichlet 卷积为
\[(f\ast g)(n)=\sum_{d\mid n}f(d)g(\frac{n}{d})
\]
性质：Dirichlet 卷积满足以下运算规律：
- 交换律 \((f * g=g * f)\) ；
- 结合律 \((f * g) * h=f * (g * h)\) ；
- 分配律 \(f * (g+h)=f * g+f * h\) ；
- \(f*\varepsilon=f\) ，其中 \(\varepsilon\) 为 Dirichlet 卷积的单位元（任何函数卷 \(\varepsilon\) 都为其本身）
例子：
\[\begin{aligned}
\varepsilon=\mu \ast 1&\iff\varepsilon(n)=\sum_{d\mid n}\mu(d)\ \\
d=1 \ast 1&\iff d(n)=\sum_{d\mid n}1\ \\
\sigma=\operatorname{id} \ast 1&\iff\sigma(n)=\sum_{d\mid n}d\ \\
\varphi=\mu \ast \operatorname{id}&\iff\varphi(n)=\sum_{d\mid n}d\cdot\mu(\frac{n}{d}) \\
id = \varphi \ast 1 &\iff id(n) = \sum_{d \mid n} \varphi(d)
\end{aligned}
\]
1.2.6.3 莫比乌斯函数
定义：\(\mu\) 为莫比乌斯函数，定义为
\[\mu(n)= \begin{cases} 1&n=1\ 0&n\text{ 含有平方因子}\ (-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数}\ \end{cases}
\]
详细解释一下：
令 \(n=\prod_{i=1}^kp_i^{c_i}\) ，其中 \(p_i\) 为质因子， \(c_i\ge 1\) 。上述定义表示：
- \(n=1\) 时， \(\mu(n)=1\) ；
- 对于 \(n\not= 1\) 时：
- 当存在 \(i\in [1,k]\) ，使得 \(c_i > 1\) 时， \(\mu(n)=0\) ，也就是说只要某个质因子出现的次数超过一次， \(\mu(n)\) 就等于 \(0\) ；
- 当任意 \(i\in[1,k]\) ，都有 \(c_i=1\) 时， \(\mu(n)=(-1)^k\) ，也就是说每个质因子都仅仅只出现过一次时，即 \(n=\prod_{i=1}^kp_i\) ， \({p_i}_{i=1}^k\) 中个元素唯一时， \(\mu(n)\) 等于 \(-1\) 的 \(k\) 次幂，此处 \(k\) 指的便是仅仅只出现过一次的质因子的总个数。
性质：莫比乌斯函数不但是积性函数，还有如下性质：
\[\sum_{d\mid n}\mu(d)= \begin{cases} 1&n=1\ 0&n\neq 1\ \end{cases}
\]
即 \(\sum_{d\mid n}\mu(d)=\varepsilon(n)\) ， \(\mu * 1 =\varepsilon\)
1.2.6.4 数论分块
若有可能存在 \(k < n\) 时，求出 \(\sum_{i=1}^n \lfloor \frac{k}{i} \rfloor\)
LL ans = 0;
for (LL l = 1, r; l <= n;l = r + 1) { //此处l意同i,r意同j,下个计算区间的l应为上个区间的r+1
if (k / l != 0) r = min(k / (k / l), n);
else r = n; // l大于k时
ans += (k / l) * (r - l + 1); //这个区间内k/i均相等,对i求和是等差数列求和
}
若不存在 \(k < n\) (这里 k = n || k = m)时，求出 \(\sum_{i=1}^{min(n,m)} \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor\)
此时可将代码中 r = n/(n/i) 替换成 r = min(n/(n/i), m/(m/i)) ，并且不需要 k / l != 0 一句判断。
1.2.6.5 莫比乌斯反演
公式：设 \(f(n),g(n)\) 为两个数论函数。
如果有 \(f(n)=\sum_{d\mid n}g(d)\) ，那么有 \(g(n)=\sum_{d\mid n}\mu(d)f(\frac{n}{d})\) 。
如果有 \(f(n)=\sum_{n|d}g(d)\) ，那么有 \(g(n)=\sum_{n|d}\mu(\frac{d}{n})f(d)\) 。
1.2.6.6 莫比乌斯反演扩展
结尾补充一个莫比乌斯反演的非卷积形式。
对于数论函数 \(f,g\) 和完全积性函数 \(t\) 且 \(t(1)=1\) ：
\[ f(n)=\sum_{i=1}^nt(i)g\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\ \iff g(n)=\sum_{i=1}^n\mu(i)t(i)f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)
\]
1.2.7 杜教筛
杜教筛被用来处理数论函数的前缀和问题。对于求解一个前缀和，杜教筛可以在低于线性时间的复杂度内求解
对于数论函数 \(f\) ，要求我们计算 \(S(n)=\sum_{i=1}^{n}f(i)\) .
我们想办法构造一个 \(S(n)\) 关于 \(S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\) 的递推式
对于任意一个数论函数 \(g\) ，必满足
\[\sum_{i=1}^{n}\sum_{d \mid i}f(d)g\left(\frac{i}{d}\right)=\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\ \Leftrightarrow \sum_{i=1}^{n}(f\ast g)(i)=\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)
\]
可以得到递推式：
\[g(1)S(n)=\sum_{i=1}^n(f\ast g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)
\]
那么假如我们可以快速对 \(\sum_{i=1}^n(f \ast g)(i)\) 求和，并用数论分块求解 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\) 就可以在较短时间内求得 \(g(1)S(n)\) .
1.2.7.1 问题一
题目大意：求 \(S_1(n)= \sum_{i=1}^{n} \mu(i)\) 和 \(S_2(n)= \sum_{i=1}^{n} \varphi(i)\) 的值， \(n\le 2^{31} -1\) 。
莫比乌斯函数前缀和
由狄利克雷卷积 ，我们知道：
\(\because \epsilon =\mu \ast 1\) （ \(\epsilon(n)=~[n=1]\) ）
\(\therefore \epsilon (n)=\sum_{d \mid n} \mu(d)\)
\(S_1(n)=\sum_{i=1}^n \epsilon (i)-\sum_{i=2}^n S_1(\lfloor \frac n i \rfloor)\)
\(= 1-\sum_{i=2}^n S_1(\lfloor \frac n i \rfloor)\)
观察到 \(\lfloor \frac n i \rfloor\) 最多只有 \(O(\sqrt n)\) 种取值，我们就可以应用 整除分块 （或称数论分块）来计算每一项的值了。
直接计算的时间复杂度为 \(O(n^{\frac 3 4})\) 。考虑先线性筛预处理出前 \(n^{\frac 2 3}\) 项，剩余部分的时间复杂度为
\(O(\int_{0}^{n^{\frac 1 3}} \sqrt{\frac{n}{x}} ~ dx)=O(n^{\frac 2 3})\)
对于较大的值，需要用 map 存下其对应的值，方便以后使用时直接使用之前计算的结果。
欧拉函数前缀和
当然也可以用杜教筛求出 \(\varphi (x)\) 的前缀和，但是更好的方法是应用莫比乌斯反演：
\(\sum_{i=1}^n \sum_{j=1}^n 1[\gcd(i,j)=1]=\sum_{i=1}^n \sum_{j=1}^n \sum_{d \mid i,d \mid j} \mu(d)\)
\(=\sum_{d=1}^n \mu(d) {\lfloor \frac n d \rfloor}^2\)
由于题目所求的是 \(\sum_{i=1}^n \sum_{j=1}^i 1[\gcd(i,j)=1]\) ，所以我们排除掉 \(i=1,j=1\) 的情况，并将结果除以 \(2\) 即可。
观察到，只需求出莫比乌斯函数的前缀和，就可以快速计算出欧拉函数的前缀和了。时间复杂度 \(O(n^{\frac 2 3})\) 。
使用杜教筛求解
求 \(S(i)=\sum_{i=1}^n\varphi(i)\) .
同样的， \(\varphi\ast 1=ID\)
\[\begin{split} &\sum_{i=1}^n(\varphi\ast 1)(i)=\sum_{i=1}^n1\cdot S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\ &\sum_{i=1}^nID(i)=\sum_{i=1}^n1\cdot S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\ &\frac{1}{2}n(n+1)=\sum_{i=1}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\ &S(n)=\frac{1}{2}n(n+1)-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\ \end{split}
\]
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
using namespace std;
const int maxn = 2000010;
typedef long long ll;
ll T, n, pri[maxn], cur, mu[maxn], sum_mu[maxn];
bool vis[maxn];
map<ll, ll> mp_mu;
ll S_mu(ll x) {
if (x < maxn) return sum_mu[x];
if (mp_mu[x]) return mp_mu[x];
ll ret = 1ll;
for (ll i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
ret -= S_mu(x / i) * (j - i + 1);
}
return mp_mu[x] = ret;
}
ll S_phi(ll x) {
ll ret = 0ll;
for (ll i = 1, j; i <= x; i = j + 1) {
j = x / (x / i);
ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
}
return ((ret - 1) >> 1) + 1;
}
int main() {
scanf("%lld", &T);
mu[1] = 1;
for (int i = 2; i < maxn; i++) {
if (!vis[i]) {
pri[++cur] = i;
mu[i] = -1;
}
for (int j = 1; j <= cur && i * pri[j] < maxn; j++) {
vis[i * pri[j]] = true;
if (i % pri[j])
mu[i * pri[j]] = -mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
for (int i = 1; i < maxn; i++) sum_mu[i] = sum_mu[i - 1] + mu[i];
while (T--) {
scanf("%lld", &n);
printf("%lld %lld\n", S_phi(n), S_mu(n));
}
return 0;
}
1.2.8 原根
1.2.8.1 求一个数最小原根（不完美）
\(O(X \frac{\log X}{\log \log X} \log X)\) 求一个数的原根。辣鸡做法，以后会有更好的
int tot = 0;
int fac[105];
IL int getroot(int x) {
tot = 0;
int phix = x - 1,m = sqrt(phix+0.5);
for(int i=2;i<=m;i++) {
if(phix % i == 0) {
fac[tot++] = i;
while(phix % i == 0) phix /= i;
}
}
if(phix > 1) fac[tot++] = phix;
phix = x-1;
for(int i=2;i<=phix;i++) {
bool flag = true;
for(int j=0;j<tot&&flag;j++)
if(ksm(i,phix/fac[j],x)==1) flag = false;
if(flag) return i;
}
return -1;
}
1.3 线性代数
1.3.1 线性基
给定 n 个整数（数字可能重复），求在这些数中选取任意个，使得他们的异或和最大。
void insert(LL x) // 插入一个数
LL qmax() // 查询最大的异或和
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
#define ri register int
#define IL inline
const int N = 100;
int n;
LL b;
LL a[N];
void insert(LL x) {
for(ri i=60;~i;i--) {
if(x&(1LL<<i)) {
if(!a[i]) {a[i]=x; return;}
else x ^= a[i];
}
}
}
LL qmax() {
LL ans = 0;
for(ri i=60;~i;i--) {
ans = max(ans,ans^a[i]);
}
return ans;
}
int main() {
scanf("%d",&n);
for(ri i=1;i<=n;i++) {
scanf("%lld",&b);
insert(b);
}
printf("%lld\n",qmax());
return 0;
}
1.3.2 单纯形法求线性规划（完全不懂怎么改的板子）
有 \(n\) 个实数变量 \(x_1,x_2,…,x_n\) 和 \(m\) 条约束，其中第 \(i\) 条约束形如 \(∑^n_{j=1}a_{ij}x_j≤b_i\)。
此外这 \(n\) 个变量需要满足非负性限制，即 \(x_j≥0\)。
在满足上述所有条件的情况下，你需要指定每个变量 \(x_j\) 的取值，使得目标函数 \(F=∑^n_{j=1}c_jx_j\) 的值最大。
输入格式
第一行三个正整数 \(n,m,t\)。其中 \(t∈\{0,1\}\)。
第二行有 \(n\) 个整数 \(c_1,c_2,⋯,c_n\)，整数间均用一个空格分隔。
接下来 \(m\) 行，每行代表一条约束，其中第 \(i\) 行有 \(n+1\) 个整数 \(a_{i1},a_{i2},⋯,a_{in},b_i\)，整数间均用一个空格分隔。
输出格式
如果不存在满足所有约束的解，仅输出一行 "Infeasible"。
如果对于任意的 \(M\)，都存在一组解使得目标函数的值大于 \(M\)，仅输出一行"Unbounded"。
否则，第一行输出一个实数，表示目标函数的最大值 \(F\)。
如果 \(t=1\)，那么你还需要输出第二行，用空格隔开的 \(n\) 个非负实数，表示此时 \(x_1,x_2,…,x_n\) 的取值，如有多组方案请任意输出其中一个。
如果 \(t=0\)，或者出现 Infeasible 或 Unbounded 时，不需要输出第二行。
此算法最坏情况下是指数级别的复杂度。
#include <bits/stdc++.h>
#ifdef __WIN32
#define LLFORMAT "I64"
#else
#define LLFORMAT "ll"
#endif
using namespace std;
#define eps 1e-7
int simplex(vector<vector<double> > &a, vector<double> &b, vector<double> &c, vector<int> &basic) {
/*
* Solving the LP problem of:
* minimize c^Tx
* with Ax=b
* basic[i] denotes the index of basic variables for a[i][...]
* Return:
* -1 if solution is unbounded, 0 otherwise
* basic variables for Ax=b will denoted by ''basic'', the corresponding solution will be basic variables taking corresponding b value and nonbasic variables taking 0
*/
int m = b.size(), n = c.size();
for (int i = 0; i < m; ++i) {
assert(b[i] > -eps);
int k = basic[i];
assert(-eps < a[i][k] - 1 && a[i][k] - 1 < eps);
for (int l = 0; l < m; ++l) if(l != i) assert(-eps < a[l][k] && a[l][k] < eps);
assert(-eps < c[k] && c[k] < eps);
}
while(true) {
int k = -1;
for (int j = 0; j < n; ++j) if(c[j] < -eps) {
k = j;
break;
}
if(k == -1) {
double ans = 0;
for (int i = 0; i < m; ++i) ans += c[basic[i]] * b[i];
return 0;
}
int l = -1;
for (int i = 0; i < m; ++i) if(a[i][k] > eps) {
if(l == -1) l = i;
else {
double ti = b[i] / a[i][k], tl = b[l] / a[l][k];
if(ti < tl - eps || (ti < tl + eps && basic[i] < basic[l])) l = i;
}
}
if(l == -1) return -1;
basic[l] = k;
double tmp = 1 / a[l][k];
for (int j = 0; j < n; ++j) a[l][j] *= tmp;
b[l] *= tmp;
for (int i = 0; i < m; ++i) if(i != l) {
tmp = a[i][k];
for (int j = 0; j < n; ++j) a[i][j] -= tmp * a[l][j];
b[i] -= tmp * b[l];
}
tmp = c[k];
for (int j = 0; j < n; ++j) c[j] -= tmp * a[l][j];
}
}
int main() {
ios::sync_with_stdio(false);
int n, m, T;
cin >> n >> m >> T;
vector<double> c(n + m, 0);
for (int i = 0; i < n; ++i) {
cin >> c[i];
c[i] *= -1;
}
auto C = c;
vector<vector<double> > a(m, vector<double>(n + m, 0));
vector<double> b(m);
vector<int> basic(m, -1), tmp;
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) cin >> a[i][j];
a[i][i + n] = 1;
cin >> b[i];
if(b[i] > -eps) basic[i] = i + n;
else tmp.push_back(i);
}
if(!tmp.empty()) {
sort(tmp.begin(), tmp.end(), [&](int i, int j) { return b[i] > b[j]; });
vector<vector<double> > A;
vector<double> B, C(n + m + 1, 0);
vector<int> Basic;
for (int i: tmp) {
vector<double> foo;
for (int j = 0; j < n + m; ++j) foo.push_back(-a[i][j]);
foo.push_back(1);
double bar = -b[i];
for (int i = 0; i < A.size(); ++i) {
double tmp = foo[Basic[i]];
for (int j = 0; j <= n + m; ++j) foo[j] -= tmp * A[i][j];
bar -= tmp * B[i];
}
for (int j = n + m; j >= 0; --j) if(-eps < foo[j] - 1 && foo[j] - 1 < eps) {
Basic.push_back(j);
break;
}
for (int i = 0; i < A.size(); ++i) {
double tmp = A[i][Basic.back()];
for (int j = 0; j <= n + m; ++j) A[i][j] -= tmp * foo[j];
B[i] -= tmp * bar;
}
A.push_back(foo);
B.push_back(bar);
assert(bar > -eps);
assert(A.size() == Basic.size());
}
for (int i = 0; i < A.size(); ++i) if(Basic[i] == n + m) {
for (int j = 0; j < n + m; ++j) C[j] = -A[i][j];
}
for (int i = 0; i < m; ++i) if(b[i] > -eps) {
A.push_back(a[i]);
A[A.size() - 1].push_back(0);
B.push_back(b[i]);
Basic.push_back(basic[i]);
}
assert(simplex(A, B, C, Basic) == 0);
bool flag = true;
for (int i = 0; i < m; ++i) if(Basic[i] == n + m) {
if(B[i] > eps) {
cout << "Infeasible\n";
return 0;
}
int k = -1;
for (int j = 0; j < n + m; ++j) if(A[i][j] > eps || A[i][j] < -eps) {
k = j;
break;
}
if(k != -1) {
double tmp = 1 / A[i][k];
Basic[i] = k;
for (int j = 0; j <= n + m; ++j) A[i][j] *= tmp;
B[i] *= tmp;
for (int l = 0; l < m; ++l) if(l != i) {
tmp = A[l][k];
for (int j = 0; j <= n + m; ++j) A[l][j] -= tmp * A[i][j];
B[l] -= tmp * B[i];
}
}
else flag = false;
break;
}
if(flag) {
A.push_back(vector<double>(n + m, 0));
A[A.size() - 1].push_back(1);
B.push_back(0);
Basic.push_back(n + m);
for (int i = 0; i < A.size() - 1; ++i) {
double tmp = A[i].back();
for (int j = 0; j <= n + m; ++j) A[i][j] -= tmp * A[A.size() - 1][j];
B[i] -= tmp * B.back();
}
}
a = A;
b = B;
basic = Basic;
c.push_back(0);
for (int i = 0; i < a.size(); ++i) {
double tmp = c[basic[i]];
for (int j = 0; j <= n + m; ++j) c[j] -= tmp * a[i][j];
}
}
auto foo = simplex(a, b, c, basic);
if(foo == -1) cout << "Unbounded" << endl;
else {
double res = 0;
vector<double> ans(n, 0);
for (int i = 0; i < basic.size(); ++i) if(basic[i] < n) ans[basic[i]] = b[i];
for (int j = 0; j < n; ++j) res -= C[j] * ans[j];
cout << setprecision(8) << res << endl;
if(T == 1) {
for (int i = 0; i < n; ++i) cout << setprecision(8) << ans[i] << ' ';
cout << endl;
}
}
return 0;
}
1.4 多项式
1.4.1 lagrange 插值
以下两份代码都是从 0 开始存的，且 n 为多项式最高次幂次数。（注意 n 不是项数）
1.4.1.1 x 连续的 lagrange 插值
坐标连续的拉格朗日插值法：
设
\[pre_i=\prod_{j=0}^i k-j \\ suf_i = \prod_{j=i}^n k-j \\ f(k) = \sum_{i=0}^n y_i \prod_{i=0,i\not=j}^n \frac{k-j}{i-j} = \sum_{i=0}^n y_i \frac{pre_{i-1}\cdot suf_{i+1}}{i! \cdot (n-i)!}
\]
//(x0,y0),(x1,y1),..,(xn,yn) -> (xi,y_xi)
LL pre[N],suf[N];
LL lagrange(int n,int *x,LL *y,LL xi) {
pre[0] = (xi-x[0]) % mod; suf[n+1] = 1;
for(int i=1;i<=n;i++) pre[i] = (xi-x[i])%mod*pre[i-1]%mod;
for(int i=n;i>=0;i--) suf[i] = (xi-x[i])%mod*suf[i+1]%mod;
LL ans = y[0]%mod*suf[1]%mod*ijc[n]%mod;
if(n&1) ans = -ans;
for(int i=1;i<=n;i++) {
LL now = y[i]*pre[i-1]%mod*suf[i+1]%mod*ijc[i]%mod*ijc[n-i]%mod;
if((n-i)&1) ans = (ans - now + mod) % mod;
else ans = (ans + now) % mod;
}
return ans;
}
1.4.1.2 n 个不连续的坐标的 lagrange 插值法
\[f(k) = \sum_{i=0}^n y_i \prod_{i=0,i\not=j}^n \frac{k-x_j}{x_i-x_j}
\]
注：这里 lagrange 函数里的 n 是多项式最高项次数。
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
#define IL inline
#define ri register int
const int mod = 998244353;
int ksm(int a,int b,int p) {
int res = 1;
while(b) {
if(b&1) res = 1LL * res * a % p;
a = 1LL * a * a % p;
b >>= 1;
}
return res;
}
int lagrange(int n,int *x,int *y,int xi) {
int ans = 0;
for(int i=0;i<=n;i++) {
int fz = 1, fm = 1;
for(int j=0;j<=n;j++) if(i != j) {
fz = 1LL * fz * (xi - x[j]) % mod;
fm = 1LL * fm * (x[i]-x[j]) % mod;
}
ans = (ans + 1LL*y[i]*fz%mod*ksm(fm,mod-2,mod) % mod) % mod;
}
return (ans + mod) % mod;
}
const int N = 2e3 + 3;
int n,m;
int a[N],b[N];
int main() {
scanf("%d%d",&n,&m);
for(int i=0;i<n;i++) scanf("%d%d",&a[i],&b[i]);
printf("%d\n",lagrange(n-1,a,b,m));
return 0;
}
1.4.2 快速傅里叶变换
给定一个 n 次多项式 F(x)，和一个 m 次多项式 G(x)。
请求出 F(x) 和 G(x) 的卷积。
即求出
\[h_i = \sum_{j+k=i} f_jg_k
\]
注意循环卷积代码未经测试，谨慎使用。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef double db;
#define IL inline
#define ri register int
const db PI = acos(-1.0);
// 虚数部分
struct com {
db x,y;
com(db x=0.0,db y=0.0):x(x),y(y) {}
};
IL com operator + (const com& a,const com& b) { return com(a.x+b.x,a.y+b.y);}
IL com operator - (const com& a,const com& b) { return com(a.x-b.x,a.y-b.y);}
IL com operator * (const com& a,const com& b) { return com(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);}
IL com operator / (const com& a,db p) { return com(a.x/p, a.y/p);}
IL com polar(db rho,db theta) { return com(rho*cos(theta), rho*sin(theta));}
// 快速傅里叶变换部分
const int N = 2097152; // 要开到给定的 2 倍以上，且必须是 2 的整数次幂的形式。
int rev[N];
com eps[N], ieps[N];
IL void initeps() {
for(ri i=0;i<N;i++) eps[i] = polar(1.0,2*PI*i/N);
ieps[0] = eps[0] = 1;
for(ri i=1;i<N;i++) ieps[i] = eps[N-i];
}
IL void cal_rev(int nA, int nB, int& lim, int& p) {
lim = 1; p = 0;
while(lim <= (nA+nB)) {lim <<= 1; ++p;}
for(ri i=0;i<lim;i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(p-1));
}
IL void trans(com *x, com* w, int n) {
for(ri i=0; i<n; i++) if(i < rev[i]) swap(x[i],x[rev[i]]);
for(ri i=2;i<=n;i<<=1) {
int l = i>>1, d = N / i;
for(ri j=0;j<n;j+=i) {
for(ri k=0;k<l;k++) {
com t = x[j+k+l] * w[d*k];
x[j+k+l] = x[j+k] - t;
x[j+k] = x[j+k] + t;
}
}
}
}
IL void dft(com* x, int n) { trans(x,eps,n);}
IL void idft(com* x, int n) { trans(x,ieps,n); for(ri i=0;i<n;i++) x[i] = x[i] / n;}
com fft_a[N], fft_b[N];
IL int mul_normal(db *C, db *A, db *B, int nA, int nB) { // 正常卷积 C = A * B
nA++; nB++; // nA 与 nB , 它们分别是 A 和 B 的最高次幂的次数.
// 如果它们是长度的话，请不要 ++.
int lim, p;
cal_rev(nA, nB, lim, p); // 如果多项式长度都相等的话，这句可以写在外面
for(ri i=0;i<nA;i++) fft_a[i] = com(A[i],0);
for(ri i=0;i<nB;i++) fft_b[i] = com(B[i],0);
dft(fft_a, lim); dft(fft_b, lim);
for(ri i=0;i<lim;i++) fft_a[i] = fft_a[i] * fft_b[i];
idft(fft_a, lim);
for(ri i=0;i<nA+nB-1;i++) C[i] = fft_a[i].x;
return nA + nB - 1; // 返回多项式的长度.
}
IL int mul_circle(db *C, db* A, db* B, int n) { // 循环卷积 C = A * B，A 和 B 长度相等
// n 是多项式的长度
// n 不是多项式最高次幂的次数.
//n++;// 多项式的最高次幂次数是 n-1.
int lim,p;
cal_rev(n, n, lim, p);
for(ri i=0;i<n;i++) fft_a[i] = com(A[i],0);
for(ri i=0;i<n;i++) fft_b[i] = com(B[i],0);
dft(fft_a, lim); dft(fft_b,lim);
for(ri i=0;i<lim;i++) fft_a[i] = fft_a[i] * fft_b[i];
idft(fft_a, lim);
for(ri i=0;i<n;i++) C[i] = fft_a[i].x + fft_a[i+n].x;
return n;
}
db f[N], g[N], ans[N];
int main() {
initeps(); // 一定要提前预处理 eps
int n,m; scanf("%d%d",&n,&m);
for(ri i=0;i<=n;i++) scanf("%lf",&f[i]);
for(ri i=0;i<=m;i++) scanf("%lf",&g[i]);
mul_normal(ans,f,g,n,m);
for(ri i=0;i<=n+m;i++) printf("%d%c",(int)(ans[i]+0.5)," \n"[i==n+m]);
return 0;
}
1.4.3 快速数论变换
给定一个 n 次多项式 F(x)，和一个 m 次多项式 G(x)。
请求出 F(x) 和 G(x) 的卷积。
即求出
\[h_i = \sum_{j+k=i} f_jg_k
\]
常见模数有 998244353, 1004535809, 469762049 它们的原根都是 3.
注意循环卷积代码未经测试，谨慎使用。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef double db;
#define IL inline
#define ri register int
const int MOD = 998244353;
const int N = 2097152; // 要开到给定的 2 倍以上，且必须是 2 的整数次幂的形式。
const int G = 3;
IL int ksm(int a,int b,int m) {
int res = 1;
while(b) {
if(b&1) res = 1LL * res * a % m;
a = 1LL * a * a % m;
b >>= 1;
}
return res;
}
IL int inv(int x) { return ksm(x,MOD-2,MOD);}
int rev[N];
int eps[N], ieps[N];
IL void initeps() {
int g = ksm(G, (MOD-1) / N, MOD), ig = inv(g);
ieps[0] = eps[0] = 1;
for(int i=1;i<N;i++) eps[i] = 1LL * eps[i-1] * g % MOD;
for(int i=1;i<N;i++) ieps[i] = 1LL * ieps[i-1] * ig % MOD;
}
IL void cal_rev(int nA, int nB, int& lim, int& p) {
lim = 1; p = 0;
while(lim <= (nA+nB)) {lim <<= 1; ++p;}
for(int i=0;i<lim;i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(p-1));
}
IL void trans(int *x, int* w, int n) {
for(ri i=0; i<n; i++) if(i < rev[i]) swap(x[i],x[rev[i]]);
for(ri i=2;i<=n;i<<=1) {
int l = i>>1, d = N / i;
for(ri j=0;j<n;j+=i) {
for(ri k=0;k<l;k++) {
int t = 1LL * x[j+k+l] * w[d*k] % MOD;
x[j+k+l] = (1LL * x[j+k] - t + MOD) % MOD;
x[j+k] = (x[j+k] + t) % MOD;
}
}
}
}
IL void dft(int* x, int n) { trans(x,eps,n);}
IL void idft(int* x, int n) {
trans(x,ieps,n);
int in = inv(n);
for(ri i=0;i<n;i++) x[i] = 1LL * x[i] * in % MOD;
}
int ntt_a[N], ntt_b[N];
IL int mul_normal(int *C, int *A, int *B, int nA, int nB) {// 正常卷积 C = A * B
nA++; nB++; // nA 与 nB , 它们分别是 A 和 B 的最高次幂的次数.
// 如果它们是长度的话，请不要 ++.
int lim, p;
cal_rev(nA, nB, lim, p); // 如果多项式长度都相等的话，这句可以写在外面
for(ri i=0;i<nA;i++) ntt_a[i] = A[i];
for(ri i=0;i<nB;i++) ntt_b[i] = B[i];
dft(ntt_a, lim); dft(ntt_b, lim);
for(ri i=0;i<lim;i++) ntt_a[i] = 1LL * ntt_a[i] * ntt_b[i] % MOD;
idft(ntt_a, lim);
for(ri i=0;i<nA+nB-1;i++) C[i] = ntt_a[i];
return nA + nB - 1; // 返回多项式的长度.
}
IL int mul_circle(int *C, int* A, int* B, int n) {// 循环卷积 C = A * B，A 和 B 长度相等
// n 是多项式的长度
// n 不是多项式最高次幂的次数.
//n++;// 多项式的最高次幂次数是 n-1.
int lim,p;
cal_rev(n, n, lim, p);
for(ri i=0;i<n;i++) ntt_a[i] = A[i];
for(ri i=0;i<n;i++) ntt_b[i] = B[i];
dft(ntt_a, lim); dft(ntt_b,lim);
for(ri i=0;i<lim;i++) ntt_a[i] = 1LL * ntt_a[i] * ntt_b[i] % MOD;
idft(ntt_a, lim);
for(ri i=0;i<n;i++) C[i] = (ntt_a[i] + ntt_a[i+n]) % MOD;
return n;
}
int f[N], g[N], ans[N];
int main() {
initeps();// 一定要提前预处理 eps
int n,m; scanf("%d%d",&n,&m);
for(ri i=0;i<=n;i++) scanf("%d",&f[i]);
for(ri i=0;i<=m;i++) scanf("%d",&g[i]);
mul_normal(ans,f,g,n,m);
for(ri i=0;i<=n+m;i++) printf("%d%c",ans[i]," \n"[i==n+m]);
return 0;
}
1.4.3 快速沃尔什变换
求出给定长度为 \(2^n\) 两个序列 A,B，设
\[C_i=\sum_{j\oplus k = i}A_j \times B_k
\]
分别当 ⊕ 是 or,and,xor 时求出 C.
注意这里的数组是从下标 0 开始存的。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
#define IL inline
#define ri register int
#define mk make_pair
const int N = 19;
const int MOD = 998244353;
const int inv2 = 499122177;
int n;
int a[1<<N],b[1<<N];
int ora[1<<N],orb[1<<N],orc[1<<N];
int anda[1<<N],andb[1<<N],andc[1<<N];
int xora[1<<N],xorb[1<<N],xorc[1<<N];
IL void OR(int *A,int n,int x) {
for(ri i=1;i<(1<<n);i<<=1) {
for(ri p=i<<1,j=0;j<(1<<n);j+=p) {
for(ri k=0;k<i;k++) {
A[i+j+k] = (1LL * A[i+j+k] + A[j+k] * x + MOD) % MOD;
}
}
}
}
IL void AND(int *A,int n,int x) {
for(ri i=1;i<(1<<n);i<<=1) {
for(ri p=i<<1,j=0;j<(1<<n);j+=p) {
for(ri k=0;k<i;k++) {
A[j+k] = (1LL * A[j+k] + A[i+j+k]*x + MOD) % MOD;
}
}
}
}
IL void XOR(int *A,int n,int t) {
for(ri i=1;i<(1<<n);i<<=1) {
for(ri p=i<<1,j=0;j<(1<<n);j+=p) {
for(ri k=0;k<i;k++) {
int x = A[j+k], y = A[i+j+k];
A[j+k] = (1LL*x+y) % MOD;
A[i+j+k] = (1LL*x-y+MOD) % MOD;
if(t == -1) {
A[j+k] = 1LL * A[j+k] * inv2 % MOD;
A[i+j+k] = 1LL * A[i+j+k] * inv2 % MOD;
}
}
}
}
}
int main() {
scanf("%d",&n);
for(ri i=0;i<(1<<n);i++) scanf("%d",a+i);
for(ri i=0;i<(1<<n);i++) scanf("%d",b+i);
memcpy(ora,a,sizeof(int)*(1<<n));
memcpy(orb,b,sizeof(int)*(1<<n));
memcpy(anda,a,sizeof(int)*(1<<n));
memcpy(andb,b,sizeof(int)*(1<<n));
memcpy(xora,a,sizeof(int)*(1<<n));
memcpy(xorb,b,sizeof(int)*(1<<n));
OR(ora,n,1); OR(orb,n,1);
for(ri i=0;i<(1<<n);i++) orc[i] = 1LL*ora[i]*orb[i]%MOD;
OR(orc,n,-1);
AND(anda,n,1); AND(andb,n,1);
for(ri i=0;i<(1<<n);i++) andc[i] = 1LL*anda[i]*andb[i]%MOD;
AND(andc,n,-1);
XOR(xora,n,1); XOR(xorb,n,1);
for(ri i=0;i<(1<<n);i++) xorc[i] = 1LL*xora[i]*xorb[i]%MOD;
XOR(xorc,n,-1);
for(ri i=0;i<(1<<n);i++) printf("%d%c",orc[i]," \n"[i==(1<<n)-1]);
for(ri i=0;i<(1<<n);i++) printf("%d%c",andc[i]," \n"[i==(1<<n)-1]);
for(ri i=0;i<(1<<n);i++) printf("%d%c",xorc[i]," \n"[i==(1<<n)-1]);
return 0;
}
1.4.4 分治 FFT
给定序列 \(g_{1\dots n - 1}\)，求序列 \(f_{0\dots n - 1}\)。
其中 \(f_i=\sum_{j=1}^if_{i-j}g_j\)，边界为 \(f_0=1\)。
答案对 \(998244353\) 取模。
做法：
不妨设 \(F(x)=\sum_{i=0}^{\infty}f_ix^i,G(x)=\sum_{i=0}^{\infty}g_ix^i\)，且 \(g_0=0\).
那么有 \(F(x)G(x)=\sum_{i=0}^{\infty}x^i\sum_{j+k=i}f_jg_k=F(x)-f_0x^0\)
即 \(F(x)G(x) \equiv F(x)-f_0 (\bmod x^n)\)
即 \(F(x) \equiv \frac{f_0}{1-G(x)} (\bmod x^n)\)
即 \(F(X) \equiv (1-G(x))^{-1} (\bmod x^n)\)
多项式求逆即可。
1.4.5 多项式全家桶
不包括任意模数NTT。
/*
Polynomial Template 多项式模板
By: NaCly_Fish
基本定义:
poly F; 声明一个多项式F
poly F.t 多项式F的次数
poly F.a[i] 多项式F的i次项系数
目前支持操作:
void print_poly(poly f) 从低到高次输出多项式f的系数
poly inverse(poly f) 求多项式f的逆
poly integral(poly f) 求多项式f的不定积分(常数项为0)
poly derivative(poly f) 求多项式f的导数
poly log(poly f) 求多项式f的对数
poly devide(poly f,poly g) 求多项式f除以g
poly mod(poly f,poly g) 计算f模g
poly exp(poly f) 求e^f (多项式指数函数)
poly power(poly f,int k) 计算f^k 要求常数项为1(多项式快速幂)
poly EXpower(poly f,int k) 计算f^k 无特殊要求
poly sqrt(poly f) 对f计算多项式开根
poly sin(poly f) 计算sin f
poly cos(poly f) 计算cos f
poly asin(poly f) 计算asin f
poly atan(poly f) 计算atan f
void reverse(poly &f) 翻转多项式f的系数
void clear(poly &f) 清空多项式f大于其次数的系数
*/
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cmath>
#include<ctime>
#define ll long long
#define N 262147
#define reg register
#define p 998244353
using namespace std;
const int LIM = 1<<18;
inline void read(int &x){
x = 0;
char c = getchar();
while(c<'0'||c>'9') c = getchar();
while(c>='0'&&c<='9'){
x = (x<<3)+(x<<1)+(c^48);
c = getchar();
}
}
void print(int x){
if(x>9) print(x/10);
putchar(x%10+'0');
}
namespace polynomial{
#define img 86583718
#define random(a,b) ((ll)rand()*rand()%(b-a+1)+a)
struct poly{
int t;
int a[N];
};
int lg2[N],rt[N],irt[N],inv[N],rev[N];
inline int power(int a,int t);
inline void NTT(poly &f,int type,int lim);
inline void print_poly(poly f);
inline void clear(poly &f);
inline void reverse(poly &f);
poly inverse(poly f);
poly log(poly f);
inline poly integral(poly f);
inline poly derivative(poly f);
poly exp(poly f);
poly divide(poly f,poly g);
poly mod(poly f,poly g);
poly power(poly f,int k);
poly sqrt(poly f);
poly sin(poly f);
poly cos(poly f);
poly asin(poly f);
poly atan(poly f);
poly EXpower(poly f,int k);
inline int mod_sqrt(int a);
inline void swap(int &x,int &y){
if(x==y) return;
x ^= y ^= x ^= y;
}
inline int add(int a,int b){
return a+b>=p?a+b-p:a+b;
}
inline int dec(int a,int b){
return a<b?a-b+p:a-b;
}
poly mod_power(poly f,int t,poly G){
poly h,g = f;
int lim;
t--;
while(t){
if(t&1){
lim = 1;
h = f;
while(lim<=h.t+g.t) lim <<= 1;
NTT(g,1,lim),NTT(h,1,lim);
for(reg int i=0;i<=lim;++i)
g.a[i] = (ll)g.a[i]*h.a[i]%p;
NTT(g,-1,lim);
g.t += h.t;
for(reg int i=g.t+1;i<=lim;++i) g.a[i] = 0;
g = mod(g,G);
}
lim = 1;
while(lim<=(f.t<<1)) lim <<= 1;
NTT(f,1,lim);
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*f.a[i]%p;
NTT(f,-1,lim);
f.t <<= 1;
for(reg int i=f.t+1;i<=lim;++i) f.a[i] = 0;
f = mod(f,G);
t >>= 1;
}
return g;
}
inline int mod_sqrt(int a){
srand((ll)a*time(0)%p);
int pw,s = p-1,t = 0;
while(!(s&1)){
++t;
s >>= 1;
}
int b = rand();
while(power(b,(p-1)/2)==1) b = rand();
int inv = power(a,p-2),x = power(a,(s+1)>>1);
for(reg int k=0;k<t;++k){
if(power((ll)x*x%p*inv%p,power(2,t-k-1))==1) continue;
pw = (ll)s*power(2,k-1)%p;
x = (ll)x*power(b,pw)%p;
}
return min(x,p-x);
}
inline void init(){
inv[1] = 1;
for(reg int i=2;i<N;++i)
inv[i] = (ll)(p-p/i)*inv[p%i]%p;
for(reg int i=2;i<N;++i) lg2[i] = lg2[i>>1]+1;
rt[0] = irt[0] = 1;
rt[1] = power(3,(p-1)/LIM);
irt[1] = power(rt[1],p-2);
for(reg int i=2;i<=LIM;++i){
rt[i] = (ll)rt[i-1]*rt[1]%p;
irt[i] = (ll)irt[i-1]*irt[1]%p;
}
}
inline poly EXpower(poly f,int k){
int u,v,low = 0,n = f.t;
for(reg int i=0;i<=n;++i){
if(!f.a[i]) continue;
low = i;
break;
}
if(low){
for(reg int i=low;i<=n;++i){
f.a[i-low] = f.a[i];
f.a[i] = 0;
}
}
v = f.a[0],u = power(v,p-2);
for(reg int i=0;i<=n;++i) f.a[i] = (ll)f.a[i]*u%p;
f = power(f,k);
u = power(v,k);
for(reg int i=0;i<=n;++i) f.a[i] = (ll)f.a[i]*u%p;
if(low){
if((ll)low*k>p) low = p;
else low *= k;
for(reg int i=n;i>=0;--i){
if(i+low<=n) f.a[i+low] = f.a[i];
f.a[i] = 0;
}
}
return f;
}
poly atan(poly f){
int lim = 1,n = f.t;
static poly g = f;
while(lim<=(n<<1)) lim <<= 1;
NTT(f,1,lim);
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*f.a[i]%p;
NTT(f,-1,lim);
for(reg int i=n+1;i<=lim;++i) f.a[i] = 0;
f.a[0] = (f.a[0]+1)%p;
f = inverse(f),g = derivative(g);
NTT(f,1,lim),NTT(g,1,lim);
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*g.a[i]%p;
NTT(f,-1,lim);
for(reg int i=n+1;i<=lim;++i) f.a[i] = 0;
return integral(f);
}
poly asin(poly f){
int lim = 1,n = f.t;
poly g = f;
while(lim<=(n<<1)) lim <<= 1;
NTT(f,1,lim);
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*f.a[i]%p;
NTT(f,-1,lim);
for(reg int i=n+1;i<=lim;++i) f.a[i] = 0;
for(int i=0;i<=n;++i) f.a[i] = dec(0,f.a[i]);
f.a[0] = add(f.a[0],1);
f = inverse(sqrt(f));
g = derivative(g);
NTT(f,1,lim),NTT(g,1,lim);
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*g.a[i]%p;
NTT(f,-1,lim);
for(reg int i=n+1;i<=lim;++i) f.a[i] = 0;
return integral(f);
}
poly log(poly f){
clear(f);
int lim = 1,n = f.t;
poly g = derivative(f);
f = inverse(f);
while(lim<=(n<<1)) lim <<= 1;
NTT(f,1,lim),NTT(g,1,lim);
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*g.a[i]%p;
NTT(f,-1,lim);
for(reg int i=n+1;i<=lim;++i) f.a[i] = 0;
return integral(f);
}
poly cos(poly f){
for(reg int i=0;i<=f.t;++i)
f.a[i] = (ll)f.a[i]*img%p;
poly g = exp(f);
f = inverse(g);
for(reg int i=0;i<=f.t;++i) f.a[i] = add(f.a[i],g.a[i]);
for(reg int i=0;i<=f.t;++i)
f.a[i] = (ll)f.a[i]*inv[2]%p;
return f;
}
poly sin(poly f){
for(reg int i=0;i<=f.t;++i)
f.a[i] = (ll)f.a[i]*img%p;
poly g = exp(f);
f = inverse(g);
for(reg int i=0;i<=f.t;++i) f.a[i] = dec(g.a[i],f.a[i]);
int w = power(img<<1,p-2);
for(reg int i=0;i<=f.t;++i)
f.a[i] = (ll)f.a[i]*w%p;
return f;
}
poly sqrt(poly f){
int n = f.t,top = 0,lim = 1;
int s[70];
poly g,h;
memset(g.a,0,sizeof(g.a));
while(n){
s[++top] = n;
n >>= 1;
}
if(f.a[0]==1) g.a[0] = 1;
else g.a[0] = mod_sqrt(f.a[0]);
while(top--){
g.t = n = s[top+1];
while(lim<=(n<<1)) lim <<= 1;
h = g;
for(reg int i=0;i<=n;++i) h.a[i] = add(h.a[i],h.a[i]);
h = inverse(h);
NTT(g,1,lim);
for(reg int i=0;i!=lim;++i)
g.a[i] = (ll)g.a[i]*g.a[i]%p;
NTT(g,-1,lim);
for(reg int i=n+1;i!=lim;++i) g.a[i] = 0;
for(reg int i=0;i<=n;++i) g.a[i] = add(g.a[i],f.a[i]);
NTT(g,1,lim),NTT(h,1,lim);
for(reg int i=0;i!=lim;++i)
g.a[i] = (ll)g.a[i]*h.a[i]%p;
NTT(g,-1,lim);
for(reg int i=n+1;i!=lim;++i) g.a[i] = 0;
}
return g;
}
poly exp(poly f){
int lim = 1,n = f.t,top = 0;
int s[70];
poly g,h;
memset(g.a,0,sizeof(g.a));
while(n){
s[++top] = n;
n >>= 1;
}
g.a[0] = 1;
while(top--){
g.t = n = s[top+1];
h = g,g = log(g);
for(reg int i=0;i<=n;++i) g.a[i] = dec(f.a[i],g.a[i]);
g.a[0] = add(g.a[0],1);
while(lim<=(n<<1)) lim <<= 1;
NTT(h,1,lim),NTT(g,1,lim);
for(reg int i=0;i!=lim;++i)
g.a[i] = (ll)g.a[i]*h.a[i]%p;
NTT(g,-1,lim);
for(reg int i=n+1;i<=lim;++i) g.a[i] = 0;
}
return g;
}
poly power(poly f,int k){
f = log(f);
for(int i=0;i<=f.t;++i)
f.a[i] = (ll)f.a[i]*k%p;
return exp(f);
}
inline void clear(poly &f){
for(reg int i=f.t+1;i<N;++i) f.a[i] = 0;
}
poly inverse(poly f){
int n = f.t,top = 0,lim = 1;
int s[70];
poly g,h,q;
memset(g.a,0,sizeof(g.a));
while(n){
s[++top] = n;
n >>= 1;
}
g.a[0] = power(f.a[0],p-2);
while(top--){
g.t = n = s[top+1];
q = g,h = f;
for(reg int i=n+1;i<=f.t;++i) h.a[i] = 0;
while(lim<=(n<<1)) lim <<= 1;
NTT(g,1,lim),NTT(h,1,lim);
for(reg int i=0;i!=lim;++i)
g.a[i] = (ll)g.a[i]*g.a[i]%p*h.a[i]%p;
NTT(g,-1,lim);
for(reg int i=n+1;i!=lim;++i) g.a[i] = 0;
for(reg int i=0;i<=n;++i)
g.a[i] = dec(add(q.a[i],q.a[i]),g.a[i]);
}
return g;
}
poly mod(poly f,poly g){
if(f.t<g.t) return f;
poly q = divide(f,g);
int n = f.t,m = g.t,lim = 1;
while(lim<=n) lim <<= 1;
NTT(g,1,lim),NTT(q,1,lim);
for(reg int i=0;i!=lim;++i) g.a[i] = (ll)g.a[i]*q.a[i]%p;
NTT(g,-1,lim);
f.t = m-1;
for(reg int i=0;i<m;++i) f.a[i] = dec(f.a[i],g.a[i]);
for(reg int i=m;i<=n;++i) f.a[i] = 0;
return f;
}
poly divide(poly f,poly g){
int n = f.t,m = g.t,lim = 1;
reverse(f),reverse(g);
f.t = g.t = n-m;
while(lim<=(f.t<<1)) lim <<= 1;
for(reg int i=f.t+1;i<=n;++i) f.a[i] = 0;
for(reg int i=g.t+1;i<=m;++i) g.a[i] = 0;
g = inverse(g);
NTT(f,1,lim),NTT(g,1,lim);
for(reg int i=0;i!=lim;++i) f.a[i] = (ll)f.a[i]*g.a[i]%p;
NTT(f,-1,lim);
for(reg int i=f.t+1;i!=lim;++i) f.a[i] = 0;
reverse(f);
return f;
}
inline void reverse(poly &f){
int l = f.t>>1;
for(reg int i=0;i<=l;++i)
swap(f.a[i],f.a[f.t-i]);
}
poly integral(poly f){
for(reg int i=f.t;i>=1;--i)
f.a[i] = (ll)f.a[i-1]*inv[i]%p;
f.a[0] = 0;
return f;
}
poly derivative(poly f){
for(reg int i=0;i<f.t;++i)
f.a[i] = (ll)f.a[i+1]*(i+1)%p;
f.a[f.t] = 0;
return f;
}
inline void print_poly(poly f){
for(int i=0;i<=f.t;++i){
print(f.a[i]);
putchar(' ');
}
putchar('\n');
}
inline void NTT(poly &f,int type,int lim){
int l = lg2[lim]-1;
for(reg int i=0;i<=lim;++i){
rev[i] = (rev[i>>1]>>1)|((i&1)<<l);
if(i>=rev[i]) continue;
swap(f.a[i],f.a[rev[i]]);
}
reg int r,w,x,y;
l = LIM>>1;
for(reg int mid=1;mid<lim;mid<<=1){
r = mid<<1;
for(reg int j=0;j<lim;j+=r){
for(reg int k=0;k<mid;++k){
w = type==1?rt[l*k]:irt[l*k];
x = f.a[j|k];
y = (ll)w*f.a[j|k|mid]%p;
f.a[j|k] = add(x,y);
f.a[j|k|mid] = dec(x,y);
}
}
l >>= 1;
}
if(type==1) return;
w = inv[lim];
for(reg int i=0;i<=lim;++i)
f.a[i] = (ll)f.a[i]*w%p;
}
inline int power(int a,int t){
int res = 1;
while(t){
if(t&1) res = (ll)res*a%p;
a = (ll)a*a%p;
t >>= 1;
}
return res;
}
#undef img
#undef random
};
using namespace polynomial;
int main(){
init();
return 0;
}
2.数据结构
2.1 平衡树
2.1.1 treap
可实现
1.插入 x 数
2.删除 x 数(若有多个相同的数，只删除一个)
3.查询 x 数的排名(排名定义为比当前数小的数的个数 +1)
4.查询排名为 x 的数
5.求 x 的前驱(前驱定义为小于 x，且最大的数)
6.求 x 的后继(后继定义为大于 x，且最小的数)
所有指针都要初始化为null，不初始化为null无法insert进去.
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
#define IL inline
typedef long long LL;
const int N = 1e5 + 3;
struct Node *null;
struct Node {
Node *ch[2];
int v,sz,r,cnt;
IL Node() {}
IL Node(int v):v(v) {ch[0]=ch[1]=null; r=rand(); sz=cnt=1;}
IL bool operator < (const Node& rhs) const { return r < rhs.r;}
IL int cmp(int x) const {
if(x == v) return -1;
return v < x;
}
IL void upd() { sz = ch[0] -> sz + ch[1] -> sz + cnt;}
};
IL void initnull() {null = new Node(); null->v=null->sz=null->r=null->cnt=0;}
IL void rotate(Node* &o,int d) {
Node *k = o->ch[d^1]; o->ch[d^1] = k->ch[d]; k->ch[d] = o;
o->upd(); k->upd(); o = k;
}
void insert(Node* &o,int x) {
if(o == null) {o = new Node(x); return;}
o->sz++;
int d = o->cmp(x);
if(d == -1) {o->cnt++; return;}
insert(o->ch[d],x);
if(o->r < o->ch[d]->r) rotate(o,d^1);
o -> upd();
}
void erase(Node* &o,int x) {
if(o == null) return;
int d = o->cmp(x);
if(d == -1) {
Node* u = o;
if(o->cnt > 1) {o->cnt--; o->sz--; return;}
if(o->ch[0] != null && o->ch[1] != null) {
int d2 = o->ch[0]->r > o->ch[1]->r;
rotate(o,d2); erase(o->ch[d2],x);
}
else {
if(o->ch[0] == null) o = o->ch[1]; else o = o->ch[0];
delete u;
}
}
else erase(o->ch[d],x);
if(o != null) o->upd();
}
int query_rank(Node *o, int x) {
if(o == null) return 1;
if(o->v == x) return o->ch[0]->sz + 1;
else if(o->v < x) return o->ch[0]->sz + o->cnt + query_rank(o->ch[1],x);
else return query_rank(o->ch[0],x);
}
int query_kth(Node* o, int k) {
if(o == null) return 0;
if(k <= o->ch[0]->sz) return query_kth(o->ch[0],k);
else if(k > o->ch[0]->sz + o->cnt)
return query_kth(o->ch[1],k - o->ch[0]->sz - o->cnt);
return o->v;
}
Node *ans,*root;
void query_pre(Node* o,int x) {
if(o == null) return;
if(o->v < x) { ans = o; query_pre(o->ch[1],x);}
else query_pre(o->ch[0],x);
}
void query_sub(Node* o,int x) {
if(o == null) return;
if(x < o->v) { ans = o; query_sub(o->ch[0],x);}
else query_sub(o->ch[1],x);
}
IL void ins(int x) {insert(root,x);}
IL void del(int x) {erase(root,x);}
IL int rank(int x) {return query_rank(root,x);}
IL int kth(int x) {return query_kth(root,x);}
IL int pre(int x) {ans=null;query_pre(root,x);return ans->v;}
IL int sub(int x) {ans=null;query_sub(root,x);return ans->v;}
IL void init() { initnull(); ans = root = null;}
int main() {
int n; scanf("%d",&n);
int op,x;
init();
for(int i=0;i<n;i++) {
scanf("%d%d",&op,&x);
switch(op) {
case 1: ins(x); break;
case 2: del(x); break;
case 3: printf("%d\n",rank(x)); break;
case 4: printf("%d\n",kth(x)); break;
case 5: printf("%d\n",pre(x)); break;
case 6: printf("%d\n",sub(x)); break;
}
}
return 0;
}
2.1.2 Splay
这个 splay 有两个哨兵节点，分别是1和n+2。其实只需要一个哨兵节点就行，但是，为了美观和对称……
本题代码是维护可以翻转的区间序列。
void splay(Node* &o,int x) 把左数第 x 个元素移动到根节点。
Node* merge(Node* left,Node* right) 把left树与right树合并。
void split(Node *o,int k,Node *&left,Node *&right) 把o这棵树前k个元素分裂到left树中，剩下的分裂到right树中。
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
using namespace std;
#define IL inline
#define ri register int
typedef long long LL;
struct Node *null;
struct Node {
Node *ch[2];
int v,sz,cnt,flip;
IL Node() {}
IL Node(int v):v(v){ch[0]=ch[1]=null;flip=0;sz=cnt=1;}
IL int cmp(int k) const {
int d = k - ch[0]->sz;
if(d == 1) return -1;
return d > 0;
}
IL void upd() { sz = cnt + ch[0]->sz + ch[1]->sz;}
IL void pushdown() {
if(flip) {
flip = 0;
swap(ch[0],ch[1]);
ch[0]->flip ^= 1;
ch[1]->flip ^= 1;
}
}
};
IL void initnull() { null = new Node(); null->sz = null->v = null->cnt = null->flip = 0;}
IL void rotate(Node *&o,int d) {
Node* k = o->ch[d^1]; o->ch[d^1] = k->ch[d]; k->ch[d] = o;
o->upd(); k->upd(); o = k;
}
void splay(Node* &o,int k) {
o->pushdown();
int d = o->cmp(k);
if(d == 1) k -= o->ch[0]->sz + o->cnt;
if(d != -1) {
Node* p = o->ch[d];
p->pushdown();
int d2 = p->cmp(k);
int k2 = (d2 == 0 ? k : k - p->ch[0]->sz - p->cnt);
if(d2 != -1) {
splay(p->ch[d2],k2);
if(d == d2) rotate(o,d^1); else rotate(o->ch[d],d);
}
rotate(o,d^1);
}
}
IL Node* merge(Node* left,Node* right) {
splay(left,left->sz);
left->ch[1] = right;
left->upd();
return left;
}
IL void split(Node* o,int k,Node *&left, Node *&right) {
splay(o,k);
left = o;
right = o->ch[1];
o->ch[1] = null;
left->upd();
}
const int N = 1e5 + 9;
int n,m,valcnt;
int val[N];
Node *root;
IL Node* build(int sz) {
if(!sz) return null;
Node *l = build(sz/2);
Node *o = new Node(val[++valcnt]);
o->ch[0] = l; o->ch[1] = build(sz-sz/2-1);
o->upd();
return o;
}
IL void init(int sz) {
initnull(); root = null;
for(int i=0;i<=n+2;i++) val[i] = i;
valcnt = 0;
root = build(sz);
}
vector<int> ans;
void print(Node *o) {
if(o == null) return;
o->pushdown();
print(o->ch[0]);
ans.push_back(o->v);
print(o->ch[1]);
}
int main() {
scanf("%d%d",&n,&m);
init(n+2);
while(m--) {
int a,b; scanf("%d%d",&a,&b);
Node *left,*mid,*right,*o;
split(root,a,left,o);
split(o,b-a+1,mid,right);
mid->flip ^= 1;
root = merge(merge(left,mid),right);
}
ans.clear();
print(root);
for(int i=1;i<ans.size()-1;i++) printf("%d ",ans[i]-1);
return 0;
}
2.2 线段树
以下代码可以完成：区间乘、区间加、以及求出区间和。
#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
using namespace std;
typedef long long LL;
const int N = 1e5 + 3;
template<typename T>
inline void read(T& x) {
x=0;T f=1;char ch = getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
x *= f;
}
LL MOD;
struct Intervaltree {
int n; LL* val;
int ql,qr,v;
LL sumv[N<<2],addv[N<<2],mulv[N<<2];
void build(int o,int L,int R) {
mulv[o]=1;
if(L == R) {sumv[o]=val[L]%MOD;return;}//sumv[o]=addv[o]=val[L]%MOD;
int M = L + (R-L) / 2;
build(o<<1,L,M); build(o<<1|1,M+1,R);
sumv[o] = (sumv[o<<1]+sumv[o<<1|1]) % MOD;
}
inline void build(int n,LL* val) {
this->n = n; this->val = val;
build(1,1,n);
}
inline void pushdown(int o,int L,int R) {
int M = L + (R-L) / 2;
sumv[o<<1]=(mulv[o]*sumv[o<<1]%MOD+(M-L+1)*addv[o]%MOD)%MOD;
sumv[o<<1|1]=(mulv[o]*sumv[o<<1|1]%MOD+(R-M)*addv[o]%MOD)%MOD;
mulv[o<<1]=mulv[o<<1]*mulv[o]%MOD;
mulv[o<<1|1]=mulv[o<<1|1]*mulv[o]%MOD;
addv[o<<1]=(addv[o<<1]*mulv[o]+addv[o])%MOD;
addv[o<<1|1]=(addv[o<<1|1]*mulv[o]+addv[o])%MOD;
mulv[o]=1; addv[o] = 0;
}
void upd_add(int o,int L,int R) {
if(ql<=L&&R<=qr) {
addv[o] = (addv[o]+v) % MOD;
sumv[o] = (sumv[o]+v*(R-L+1)) % MOD;
return ;
}
pushdown(o,L,R);
int M = L + (R-L) / 2;
if(ql <= M) upd_add(o<<1,L,M);
if(M < qr) upd_add(o<<1|1,M+1,R);
sumv[o] = (sumv[o<<1]+sumv[o<<1|1]) % MOD;
}
void upd_mul(int o,int L,int R) {
if(ql<=L&&R<=qr) {
mulv[o] = (mulv[o]*v) % MOD;
addv[o] = (addv[o]*v) % MOD;
sumv[o] = (sumv[o]*v) % MOD;
return ;
}
pushdown(o,L,R);
int M = L + (R-L) / 2;
if(ql <= M) upd_mul(o<<1,L,M);
if(M < qr) upd_mul(o<<1|1,M+1,R);
sumv[o] = (sumv[o<<1]+sumv[o<<1|1]) % MOD;
}
LL que(int o,int L,int R) {
if(ql <= L && R <= qr) return sumv[o];
pushdown(o,L,R);
LL ans = 0;
int M = L + (R-L) / 2;
if(ql <= M) ans = (ans + que(o<<1,L,M)) % MOD;
if(M < qr) ans = (ans + que(o<<1|1,M+1,R)) % MOD;
return ans;
}
};
Intervaltree solver;
int n,m;
LL a[N];
int main() {
read(n); read(m); read(MOD);
for(int i=1;i<=n;i++) read(a[i]);
for(int i=1;i<=n;i++) a[i] %= MOD;
solver.build(n,a);
while(m--) {
int op,x,y,v; read(op); read(x); read(y);
solver.ql=x,solver.qr=y;
if(op!=3) read(v),solver.v=v;
switch(op) {
case 1: solver.upd_mul(1,1,n); break;
case 2: solver.upd_add(1,1,n); break;
case 3: printf("%lld\n",solver.que(1,1,n)); break;
}
}
return 0;
}
2.2.1 李超线段树
本题代码强制在线了。
要求在平面直角坐标系下维护两个操作：
- 在平面上加入一条线段。记第 i 条被插入的线段的标号为 i。
- 给定一个数 k，询问与直线 x=k 相交的线段中，交点纵坐标最大的线段的编号。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
#define IL inline
#define ri register int
const int N = 1e5 + 3;
const int M1 = 39989;
const int M2 = 1e9;
int n,lastans;
struct LINE {
double k,b;
LINE(double k=0.0,double b=0.0):k(k),b(b) {}
}l[N];
IL double f(int i,int x) { return l[i].k * x + l[i].b;}
IL double inter(int a,int b) { return (l[b].b-l[a].b) / (l[a].k-l[b].k);}
IL int cmp(int x,int a,int b) { return f(a,x) > f(b,x) ? a : b;}
struct LiChaoSegmentTree {
int lazy[M1<<2];
void build(int o,int L,int R) {
lazy[o] = 0;
if(L == R) return;
int M = L + (R-L) / 2;
build(o<<1,L,M); build(o<<1|1,M+1,R);
}
int qL,qR;
void upd(int o,int L,int R,int v) {
if(qL <= L && R <= qR) {
if(!lazy[o]) {lazy[o] = v; return;}
double ly0=f(v,L),ry0=f(v,R),ly1=f(lazy[o],L),ry1=f(lazy[o],R);
if(ly0<=ly1 && ry0<=ry1) return;
if(ly0>=ly1 && ry0>=ry1) {lazy[o]=v; return;}
double x = inter(lazy[o],v);
int M = L + (R-L) / 2;
if(ly0 >= ly1) {
if(x <= M) upd(o<<1,L,M,v);
else upd(o<<1|1,M+1,R,lazy[o]), lazy[o] = v;
}
else {
if(x > M) upd(o<<1|1,M+1,R,v);
else upd(o<<1,L,M,lazy[o]), lazy[o] = v;
}
return;
}
int M = L + (R-L) / 2;
if(qL <= M) upd(o<<1,L,M,v);
if(M < qR) upd(o<<1|1,M+1,R,v);
}
int pos;
int query(int o,int L,int R) {
int ans = 0;
if(lazy[o]) ans = lazy[o];
if(L == R) return ans;
int M = L + (R-L) / 2;
if(pos <= M) ans = cmp(pos,ans,query(o<<1,L,M));
if(M < pos) ans = cmp(pos,ans,query(o<<1|1,M+1,R));
return ans;
}
};
LiChaoSegmentTree lpr;
int main() {
lastans = 0;
scanf("%d",&n);
int cnt = 0;
for(int i=1;i<=n;i++) {
int op; scanf("%d",&op);
if(op == 0) {
int x1; scanf("%d",&x1);
x1 = (x1 + lastans - 1) % M1 + 1; lpr.pos = x1;
printf("%d\n",lastans=lpr.query(1,1,M1));
}
else {
int x0,y0,x1,y1; scanf("%d%d%d%d",&x0,&y0,&x1,&y1);
x0 = (x0+lastans-1)%M1+1; x1 = (x1+lastans-1)%M1+1;
y0 = (y0+lastans-1)%M2+1; y1 = (y1+lastans-1)%M2+1;
if(x0 == x1) l[++cnt] = LINE(0.0,max(y0,y1));
else l[++cnt] = LINE((double)(y1-y0)/(x1-x0),y0-(double)(y1-y0)/(x1-x0)*x0);
lpr.qL = min(x0,x1); lpr.qR = max(x0,x1);
lpr.upd(1,1,M1,cnt);
}
}
return 0;
}
2.2.2 主席树
静态区间第 k 小。
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
#define IL inline
#define ri register int
void read(int &x) {
x = 0;
int f=1; char ch = getchar();
while(!isdigit(ch)) {if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
x *= f;
}
const int N = 2e5 + 3;
int n,q,m;
int a[N],b[N];
struct PresidentSegmentTree {
int tot;
int rt[N];
int numv[N<<5],ls[N<<5],rs[N<<5];
int build(int L,int R) {
int root = ++tot;
numv[root] = 0;
if(L == R) return root;
int M = L + R >> 1;
ls[root] = build(L,M);
rs[root] = build(M+1,R);
return root;
}
int pos;
int upd(int pre,int L,int R) {
int root = ++tot;
ls[root]=ls[pre]; rs[root]=rs[pre]; numv[root]=numv[pre]+1;
if(L == R) return root;
int M = L+R >> 1;
if(pos <= M) ls[root] = upd(ls[pre],L,M);
else rs[root] = upd(rs[pre],M+1,R);
return root;
}
int query(int u,int v,int L,int R,int k) {
if(L == R) return L;
int x = numv[ls[v]] - numv[ls[u]];
int M = L+R >> 1;
if(x >= k) return query(ls[u],ls[v],L,M,k);
else return query(rs[u],rs[v],M+1,R,k-x);
}
}lpr;
int main() {
read(n); read(q);
for(ri i=1;i<=n;i++) { read(a[i]); b[i] = a[i];}
sort(b+1,b+1+n);
m = unique(b+1,b+1+n) - b - 1;
lpr.tot = 0;
lpr.rt[0] = lpr.build(1,m);
for(ri i=1;i<=n;i++) {
lpr.pos = lower_bound(b+1,b+1+m,a[i]) - b;
lpr.rt[i] = lpr.upd(lpr.rt[i-1],1,m);
}
while(q--) {
int l,r,k; read(l); read(r); read(k);
printf("%d\n",b[lpr.query(lpr.rt[l-1],lpr.rt[r],1,m,k)]);
}
return 0;
}
2.2.3 线段树分裂 & 合并
给出一个可重集 a（编号为 1），它支持以下操作：
0 p x y：将可重集 p 中大于等于 x 且小于等于 y 的值放入一个新的可重集中（新可重集编号为从 2 开始的正整数，是上一次产生的新可重集的编号+1）。
1 p t：将可重集 tt 中的数放入可重集 p，且清空可重集 t（数据保证在此后的操作中不会出现可重集 tt）。
2 p x q：在 p 这个可重集中加入 x 个数字 q。
3 p x y：查询可重集 p 中大于等于 x 且小于等于 y 的值的个数。
4 p k：查询在 p 这个可重集中第 k 小的数，不存在时输出 -1。
#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
using namespace std;
typedef long long LL;
#define IL inline
#define ri register int
IL void read(int &x) {
int f=1; x=0; char ch = getchar();
while(!isdigit(ch)) {if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)) {x=x*10+ch-'0'; ch=getchar();}
x *= f;
}
const int N = 2e5 + 3;
int n,m,totrt = 0, totrab = 0, totnode = 0;
int rt[N];
int rab[N<<5],ls[N<<5],rs[N<<5];
LL numv[N<<5];
IL int newnode() { return totrab ? rab[totrab--] : ++totnode;}
IL void del(int o) { rab[++totrab] = o; ls[o]=rs[o]=numv[o]=0; return;}
int upd(int o,int L,int R,int pos,int v) {
if(!o) o = newnode();
numv[o] += v;
if(L == R) { return o;}
int M = L+R >> 1;
if(pos <= M) ls[o] = upd(ls[o],L,M,pos,v);
else rs[o] = upd(rs[o],M+1,R,pos,v);
return o;
}
LL query(int o,int L,int R,int qL,int qR) {
if(!o) return 0;
if(qL <= L && R <= qR) return numv[o];
int M = L+R >> 1;
LL ans = 0;
if(qL <= M) ans += query(ls[o],L,M,qL,qR);
if(M < qR) ans += query(rs[o],M+1,R,qL,qR);
return ans;
}
int merge(int u,int v,int L,int R) {
if(!u||!v) return u|v;
numv[u] += numv[v];
if(L == R) return u;
int M = L+R >> 1;
ls[u] = merge(ls[u],ls[v],L,M);
rs[u] = merge(rs[u],rs[v],M+1,R);
del(v);
return u;
}
int kth(int o,int L,int R,int k) {
if(L == R) return L;
int M = L+R >> 1;
if(numv[ls[o]] >= k) return kth(ls[o],L,M,k);
else return kth(rs[o],M+1,R,k-numv[ls[o]]);
}
int split(int u,LL k,int L,int R) {
int root = newnode();
int M = L+R >> 1;
if(k > numv[ls[u]]) rs[root] = split(rs[u],k-numv[ls[u]],L,M);
else swap(rs[u],rs[root]);
if(k < numv[ls[u]]) ls[root] = split(ls[u],k,M+1,R);
numv[root] = numv[u] - k;
numv[u] = k;
return root;
} // <= k -> the tree of u, > k -> the tree of root.
int main() {
totnode = totrab = 0; totrt = 1;
read(n); read(m);
for(ri i=1;i<=n;i++) {
int x; read(x);
rt[1] = upd(rt[1],1,n,i,x);
}
while(m--) {
int op,x,y,z;
read(op);
if(op == 0) {
read(x); read(y); read(z);
LL k1 = query(rt[x],1,n,1,z), k2 = query(rt[x],1,n,y,z);
int tmp;
rt[++totrt] = split(rt[x],k1-k2,1,n);
tmp = split(rt[totrt],k2,1,n);
rt[x] = merge(rt[x],tmp,1,n);
}
else if(op == 1) {
read(x); read(y);
rt[x] = merge(rt[x],rt[y],1,n);
}
else if(op == 2) {
read(x); read(y); read(z);
rt[x] = upd(rt[x],1,n,z,y);
}
else if(op == 3) {
read(x); read(y); read(z);
printf("%lld\n",query(rt[x],1,n,y,z));
}
else if(op == 4) {
read(x); read(y);
if(numv[rt[x]] < y) printf("-1\n");
else printf("%d\n",kth(rt[x],1,n,y));
}
}
return 0;
}
3.图论
3.1 树上问题
3.1.1 树链剖分
如题，已知一棵包含 N 个结点的树（连通且无环），每个节点上包含一个数值，需要支持以下操作：
操作 1： 格式： 1 x y z 表示将树从 x 到 y 结点最短路径上所有节点的值都加上 z。
操作 2： 格式： 2 x y 表示求树从 x 到 y 结点最短路径上所有节点的值之和。
操作 3： 格式： 3 x z 表示将以 x 为根节点的子树内所有节点值都加上 z。
操作 4： 格式： 4 x 表示求以 x 为根节点的子树内所有节点值之和
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long LL;
#define IL inline
#define ri register int
const int N = 1e5 + 3;
int n,m,root;
LL MOD;
LL a[N],b[N];
struct SegmentTree {
LL sumv[N<<2],addv[N<<2];
void build(int o,int L,int R) {
addv[o] = sumv[o] = 0;
if(L == R) { sumv[o] = addv[o] = a[L]; return;}
int M = L + (R-L) / 2;
build(o<<1,L,M); build(o<<1|1,M+1,R);
sumv[o] = sumv[o<<1] + sumv[o<<1|1];
}
IL void pushdown(int o,int L,int R) {
int M = L + (R-L) / 2;
addv[o<<1] = (addv[o<<1] + addv[o]) % MOD;
addv[o<<1|1] = (addv[o<<1|1] + addv[o]) % MOD;
sumv[o<<1] = (sumv[o<<1] + addv[o]*(M-L+1)) % MOD;
sumv[o<<1|1] = (sumv[o<<1|1] + addv[o]*(R-M)) % MOD;
addv[o] = 0;
}
int qL,qR;
LL v;
void upd(int o,int L,int R) {
if(qL <= L && R <= qR) { addv[o] += v; sumv[o] = (sumv[o] + v*(R-L+1)) % MOD; return;}
pushdown(o,L,R);
int M = L + (R-L) / 2;
if(qL <= M) upd(o<<1,L,M);
if(M < qR) upd(o<<1|1,M+1,R);
sumv[o] = sumv[o<<1] + sumv[o<<1|1];
}
LL query(int o,int L,int R) {
if(qL <= L && R <= qR) { return sumv[o];}
pushdown(o,L,R);
int M = L + (R-L) / 2;
LL ans = 0;
if(qL <= M) ans = (ans + query(o<<1,L,M)) % MOD;
if(M < qR) ans = (ans + query(o<<1|1,M+1,R)) % MOD;
return ans;
}
}lpr;
int dfs_clock=0;
int dep[N],sz[N],son[N],pa[N];
int dfn[N],top[N];
vector<int> G[N];
void dfs1(int u,int fa,int nowd) {
dep[u] = nowd;
sz[u] = 1;
son[u] = 0;
pa[u] = fa;
for(ri i=0;i<G[u].size();i++) {
int v = G[u][i];
if(v == fa) continue;
dfs1(v,u,nowd+1);
sz[u] += sz[v];
if(sz[v] > sz[son[u]]) son[u] = v;
}
}
// get dep,sz,max_son,pa;
void dfs2(int u,int nowtop) {
dfn[u] = ++dfs_clock;
a[dfs_clock] = b[u];
top[u] = nowtop;
if(son[u]) dfs2(son[u],nowtop);
for(ri i=0;i<G[u].size();i++) {
int v = G[u][i];
if(v == pa[u] || v == son[u]) continue;
dfs2(v,v);
}
}
//get dfn,top
IL void _1() {
int x,y; scanf("%d%d%lld",&x,&y,&lpr.v);
lpr.v %= MOD;
while(top[x] != top[y]) {
if(dep[top[x]] < dep[top[y]]) swap(x,y);
lpr.qL = dfn[top[x]], lpr.qR = dfn[x];
lpr.upd(1,1,n);
x = pa[top[x]];
}
if(dep[x] > dep[y]) swap(x,y);
lpr.qL = dfn[x], lpr.qR = dfn[y];
lpr.upd(1,1,n);
}
IL LL _2() {
LL ans = 0;
int x,y; scanf("%d%d",&x,&y);
while(top[x] != top[y]) {
if(dep[top[x]] < dep[top[y]]) swap(x,y);
lpr.qL = dfn[top[x]], lpr.qR = dfn[x];
ans = (ans + lpr.query(1,1,n)) % MOD;
x = pa[top[x]];
}
if(dep[x] > dep[y]) swap(x,y);
lpr.qL = dfn[x], lpr.qR = dfn[y];
ans = (ans + lpr.query(1,1,n)) % MOD;
return ans;
}
IL void _3() {
int x; scanf("%d%lld",&x,&lpr.v);
lpr.qL = dfn[x]; lpr.qR = dfn[x] + sz[x] - 1;
lpr.v %= MOD;
lpr.upd(1,1,n);
}
IL LL _4() {
int x; scanf("%d",&x);
lpr.qL = dfn[x]; lpr.qR = dfn[x] + sz[x] - 1;
return lpr.query(1,1,n);
}
int main() {
dfs_clock = 0;
scanf("%d%d%d%lld",&n,&m,&root,&MOD);
for(ri i=1;i<=n;i++) scanf("%lld",&b[i]);
for(ri i=1;i<n;i++) {
int u,v; scanf("%d%d",&u,&v);
G[u].push_back(v);
G[v].push_back(u);
}
dfs1(root,0,1);
dfs2(root,root);
lpr.build(1,1,n);
while(m--) {
int op; scanf("%d",&op);
if(op == 1) _1();
if(op == 2) printf("%lld\n",_2());
if(op == 3) _3();
if(op == 4) printf("%lld\n",_4());
}
return 0;
}
3.1.2 虚树
建虚树需要先树链剖分再建虚树。
虚树中只有关键点和 dfs 序相邻的关键点的 LCA 。
虚树的边权需要仔细考虑清楚。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long LL;
typedef pair<int,int> pii;
#define se second
#define fi first
#define pb push_back
#define IL inline
#define ri register int
struct Edge {
int u,v,w;
Edge(int u=0,int v=0,int w=0):u(u),v(v),w(w) {}
};
const int N = 250000 + 3;
const LL INFL = 1e18 + 7;
int n,m;
int dfs_clock = 0;
int dep[N],sz[N],pa[N],son[N];
int dfn[N],top[N];
LL minv[N];
vector<Edge> G[N];
void dfs1(int u,int f,int nowd,LL noww) {
dep[u] = nowd;
sz[u] = 1;
pa[u] = f;
son[u] = 0;
minv[u] = noww;
for(auto e : G[u]) {
int v = e.v, w = e.w;
if(v == f) continue;
dfs1(v,u,nowd+1,min(noww,(LL)w));
sz[u] += sz[v];
if(sz[v] > sz[son[u]]) son[u] = v;
}
}
void dfs2(int u,int nowtop) {
dfn[u] = ++dfs_clock;
top[u] = nowtop;
if(son[u]) dfs2(son[u],nowtop);
for(auto e : G[u]) {
int v = e.v;
if(v == pa[u] || v == son[u]) continue;
dfs2(v,v);
}
}
int LCA(int u,int v) {
int t1 = top[u], t2 = top[v];
while(top[u] != top[v]) {
if(dep[top[u]] < dep[top[v]]) {
swap(u,v);
}
u = pa[top[u]];
}
return dep[u] > dep[v] ? v : u;
}
int h[N];
int sta[N];
IL bool cmp(int a,int b) { return dfn[a] < dfn[b];}
IL void build_virtual_tree() {
int statop;
sort(h+1,h+1+m,cmp);
sta[statop = 1] = 1; G[1].clear();
for(ri i=1,l;i<=m;i++) {
if(h[i] != 1) {
l = LCA(h[i],sta[statop]);
if(l != sta[statop]) {
while(dfn[l] < dfn[sta[statop-1]]) {
G[sta[statop-1]].pb(Edge(sta[statop-1],sta[statop]));
statop--;
}
if(dfn[l] > dfn[sta[statop-1]]) {
G[l].clear();
G[l].pb(Edge(l,sta[statop]));
sta[statop] = l;
}
else G[l].pb(Edge(l,sta[statop--]));
}
G[h[i]].clear(); sta[++statop] = h[i];
}
}
for(ri i=1;i<statop;i++) G[sta[i]].pb(Edge(sta[i],sta[i+1]));
return;
}
bool vip[N];
LL dp(int u) {
LL ans = 0, sum = 0;
for(auto e : G[u]) {
int v = e.v;
sum += dp(v);
}
if(vip[u]) ans = minv[u];
else ans = min(minv[u],sum);
return ans;
}
int main() {
int n; scanf("%d",&n);
for(ri i=1;i<n;i++) {
int u,v,w; scanf("%d%d%d",&u,&v,&w);
G[u].pb(Edge(u,v,w));
G[v].pb(Edge(v,u,w));
}
dfs1(1,0,1,INFL); dfs2(1,1);
int T; scanf("%d",&T);
while(T--) {
scanf("%d",&m);
for(ri i=1;i<=m;i++) scanf("%d",h+i);
for(ri i=1;i<=m;i++) vip[h[i]] = true;
build_virtual_tree();
printf("%lld\n",dp(1));
for(ri i=1;i<=m;i++) vip[h[i]] = false;
}
return 0;
}
3.1.3 求出每个子树重心
从这棵树的叶子节点开始往上推出各个子树的重心。关于重心，我们有一条性质：以树的重心为根时，所有子树的大小都不超过整棵树大小的一半。假设我们当前要开始找以 u 节点为根的子树重心，已经把 u 的儿子们的重心都找出来了，那么以 u 节点为根的子树重心一定在重儿子子树重心到 u 的路径上。重心可能有两个。
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
using namespace std;
const int N = 2e5 + 3;
vector<int> G[N];
int pa[N],ans[N],sz[N];
void dfs(int u,int fa) {
pa[u] = fa; ans[u] = u; sz[u] = 1;
for(int i=0;i<G[u].size();i++) {
int v = G[u][i];
if(v == fa) continue;
dfs(v,u);
sz[u] += sz[v];
}
for(int i=0;i<G[u].size();i++) {
int v = G[u][i], ansv = ans[v];
if(v == fa) continue;
while(sz[ansv]*2 < sz[u]) ansv = pa[ansv];
if(sz[ansv] < sz[ans[u]]) ans[u] = ansv;
}
}
int main() {
int n; scanf("%d",&n);
for(int i=0,u,v;i<n-1;i++) {
scanf("%d%d",&u,&v); u--; v--;
G[u].push_back(v); G[v].push_back(u);
}
dfs(0,-1);
for(int i=0;i<n;i++) {
if(sz[ans[i]]*2 == sz[i]) {
int a=ans[i], b=pa[ans[i]];
if(a > b) swap(a,b);
printf("%d %d\n",a+1,b+1);
}
else printf("%d\n",ans[i]+1);
}
return 0;
}
3.1.4 点分治
时间复杂度 \(O(n \log^2 n)\) ，点分治本身一个 log ，套了树状数组一个 log。
下面的代码是求出给定一棵 n 个节点的树，每条边有边权，求出树上两点距离小于等于 k 的点对数量。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<queue>
using namespace std;
typedef long long LL;
#define IL inline
#define ri register int
#define mk make_pair
#define pb push_back
#define fi first
#define se second
const int N = 4e4 + 5;
int n,K;
int BIT[N];
IL int lb(int x) {return x & (-x);}
IL int sum(int x) {
int ans = 0;
for(;x;x-=lb(x)) ans += BIT[x];
return ans;
}
IL void add(int x,int d) {
for(;x<=K;x+=lb(x)) BIT[x] += d;
}
struct Edge {
int u,v,w;
Edge(int u=0,int v=0,int w=0):u(u),v(v),w(w){}
};
vector<Edge> G[N];
int sz[N];
bool vis[N];
void dfs_cen(int u,int f,int tot,int &cen) {
sz[u] = 1;
int now = 0;
for(auto e : G[u]) {
int v = e.v;
if(v == f || vis[v]) continue;
dfs_cen(v,u,tot,cen);
sz[u] += sz[v];
now = max(now,sz[v]);
}
now = max(now,tot - sz[u]);
if(now * 2 <= tot) cen = u;
}
IL int get_cen(int u,int tot) { int cen; dfs_cen(u,0,tot,cen); return cen;}
void dfs_sz(int u,int f) {
sz[u] = 1;
for(auto e : G[u]) {
int v = e.v, w = e.w;
if(v == f || vis[v]) continue;
dfs_sz(v,u);
sz[u] += sz[v];
}
}
int d[N];
void dfs_dis(int u,int f,int dis,int &cnt) {
d[++cnt] = dis;
for(auto e : G[u]) {
int v = e.v, w = e.w;
if(v == f || vis[v]) continue;
dfs_dis(v,u,dis+w,cnt);
}
}
queue<int> tag;
int work(int u,int tot) {
u = get_cen(u,tot);
dfs_sz(u,0);
vis[u] = true;
//printf("work : %d %d\n",u,tot);
//solve
int ans = 0;
for(auto e : G[u]) {
int v = e.v, w = e.w;
if(vis[v]) continue;
int cnt = 0;
dfs_dis(v,u,w,cnt);
for(ri i=1;i<=cnt;i++) {
if(d[i] > K) continue;
ans += sum(K-d[i]) + 1;
}
for(ri i=1;i<=cnt;i++) {
add(d[i],+1);
tag.push(d[i]);
}
}
while(!tag.empty()) {
int x = tag.front(); tag.pop();
add(x,-1);
}
for(auto e : G[u]) {
int v = e.v;
if(vis[v]) continue;
ans += work(v,sz[v]);
}
return ans;
}
int main() {
scanf("%d",&n);
for(ri i=1;i<n;i++) {
int u,v,w; scanf("%d%d%d",&u,&v,&w);
G[u].pb(Edge(u,v,w));
G[v].pb(Edge(v,u,w));
}
scanf("%d",&K);
printf("%d\n",work(1,n));
return 0;
}
3.2 图的匹配问题
3.2.1 二分图最大匹配
3.2.1.1 使用邻接矩阵实现的匈牙利算法
时间复杂度最坏 \(O(|V|^3)\)
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int N = 1000 + 3;
struct Hungary {
int n,m;
int Left[N];
bool S[N],T[N],G[N][N];
void init(int n,int m) {
this->n = n; this->m = m;
fill(Left,Left+m+1,0);
fill(S,S+n+1,0);
fill(T,T+m+1,0);
for(int i=1;i<=n;i++) fill(G[i],G[i]+m+1,0);
}
void AddEdge(int u,int v) { G[u][v] = 1;}
bool match(int i) {
S[i] = true;
for(int j=1;j<=m;j++) if(G[i][j] && !T[j]) {
T[j] = true;
if(Left[j] == 0 || match(Left[j])) {
Left[j] = i;
return true;
}
}
return false;
}
int solve() {
int ans = 0;
fill(Left,Left+m+1,0);
for(int i=1;i<=n;i++) {
fill(T,T+m+1,0);
if(match(i)) ans++;
}
return ans;
}
};
Hungary solver;
int main() {
int n,m,e; scanf("%d%d%d",&n,&m,&e);
solver.init(n,m);
while(e--) {
int u,v; scanf("%d%d",&u,&v);
if(u > n || u<1 || v > m || v < 1) continue;
solver.AddEdge(u,v);
}
printf("%d\n",solver.solve());
return 0;
}
3.2.1.2 使用邻接表实现的匈牙利算法
时间复杂度 \(O(VE)\)
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
typedef vector<int> vi;
#define IL inline
#define ri register int
#define pb push_back
#define mk make_pair
#define fi first
#define se second
const int N = 1000 + 3;
struct Hungary {
int n,m;
int Left[N];
bool S[N],T[N];
vector<int> G[N];
void init(int n,int m) {
this->n = n; this->m = m;
fill(Left,Left+m+1,0);
fill(S,S+n+1,0);
fill(T,T+m+1,0);
for(int i=1;i<=n;i++) G[i].clear();
}
void AddEdge(int u,int v) { G[u].pb(v);}
bool match(int u) {
S[u] = true;
for(auto v : G[u]) if(!T[v]) {
T[v] = true;
if(Left[v] == 0 || match(Left[v])) {
Left[v] = u;
return true;
}
}
return false;
}
int solve() {
int ans = 0;
fill(Left,Left+m+1,0);
for(int i=1;i<=n;i++) {
fill(T,T+m+1,0);
if(match(i)) ans++;
}
return ans;
}
};
Hungary solver;
int main() {
int n,m,e; scanf("%d%d%d",&n,&m,&e);
solver.init(n,m);
while(e--) {
int u,v; scanf("%d%d",&u,&v);
if(u > n || u<1 || v > m || v < 1) continue;
solver.AddEdge(u,v);
}
printf("%d\n",solver.solve());
return 0;
}
3.2.1.3 使用 Dinic 算法的建模方法
将左边所有点接上源点，右边所有点接上汇点，容量皆为 1 。原来的每条边从左往右连边，容量也皆为 1 ，最大流即最大匹配。
如果使用 Dinic 算法求该网络的最大流，可在 \(O(\sqrt n m)\) 求出。
3.2.1.4 二分图最大独立集
选最多的点，满足两两之间没有边相连。
二分图中，最大独立集 = n - 最大匹配。
3.2.1.5 二分图最小点覆盖
选最少的点，满足每条边至少有一个端点被选，不难发现补集是独立集。
二分图中，最小点覆盖 = n - 最大独立集 = 最大匹配。
3.2.2 二分图最佳完美匹配
\(O(n^3)\) bfs-KM 算法。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define IL inline
typedef long long LL;
const int N = 400 + 3;
const int INF = 0x3f3f3f3f;
struct Kuhn_Munkers {
int n;
int W[N][N];
int Lx[N],Ly[N];
int left[N];
int slack[N];
int pre[N];
bool T[N];
IL void init(int n) {
this->n = n;
for(int i=1;i<=n;i++) fill(W[i],W[i]+1+n,INF);
}
IL void bfs(int u) {
fill(slack,slack+1+n,INF);
fill(pre,pre+1+n,0);
int x,y=0,yy=0,a;
left[y] = u;
for(;;) {
x = left[y]; a = INF, T[y] = true;
for(int i=1;i<=n;i++) if(!T[i]){
if(slack[i] > Lx[x]+Ly[i]-W[x][i]) {
slack[i] = Lx[x] + Ly[i] - W[x][i];
pre[i] = y;
}
if(slack[i] < a) a = slack[i],yy = i;
}
for(int i=0;i<=n;i++) {
if(T[i]) { Lx[left[i]] -= a; Ly[i] += a;}
else slack[i] -= a;
}
y = yy;
if(!left[y]) break;
}
while(y) left[y] = left[pre[y]], y = pre[y];
}
IL int KM() {
fill(Lx,Lx+1+n,0);
fill(Ly,Ly+1+n,0);
fill(left,left+1+n,0);
for(int i=1;i<=n;i++) {
fill(T,T+1+n,false);
bfs(i);
}
int ans = 0LL;
for(int j=1;j<=n;j++) ans += W[left[j]][j];
return ans;
}
}solver;
int n;
int p[N];
LL a[N],b[N],c[N];
int main() {
scanf("%d",&n); solver.init(n);
for(int i=1;i<=n;i++) scanf("%lld",&a[i]);
for(int i=1;i<=n;i++) scanf("%d",&p[i]);
for(int i=1;i<=n;i++) scanf("%lld",&b[i]);
for(int i=1;i<=n;i++) scanf("%lld",&c[i]);
for(int i=1;i<=n;i++) {
for(int j=1;j<=n;j++) {
int v = 0;
for(int k=1;k<=n;k++) if(b[i]+c[j] > a[k]) v += p[k];
solver.W[i][j] = v;
}
}
printf("%d\n",solver.KM());
return 0;
}
3.2.3 一般图最大匹配(完全不懂怎么改的板子)
带花树算法 \(2 \leq|V| \leq 500, 1 \leq |E| \leq 124750\) 无法分析复杂度。
输出一般图最大匹配的边数。
输出一般图匹配中每个点匹配后对应的另一个点。
#include <cstdio>
#include <cstring>
#include<cstdlib>
#include <algorithm>
#define M 250010
using namespace std;
char inp[33554432], *inpch = inp;
int Head[M], Next[M], Go[M], Pre[510], Nxt[510], F[510], S[510], Q[510], Vis[510], *Top = Q, Cnt = 0, Tim = 0, n, m, x, y;
inline void addedge(int x, int y)
{
Go[++Cnt] = y;
Next[Cnt] = Head[x];
Head[x] = Cnt;
}
int find(int x)
{
return x == F[x] ? x : F[x] = find(F[x]);
}
int lca(int x, int y)
{
for(Tim++, x = find(x), y = find(y); ; x ^= y ^= x ^= y)
if(x)
{
if(Vis[x] == Tim) return x;
Vis[x] = Tim;
x = find(Pre[Nxt[x]]);
}
}
void blossom(int x, int y, int l)
{
while(find(x) != l)
{
Pre[x] = y, y = Nxt[x];
if (S[y]==-1)exit(233);
S[*Top = y] = 0, *Top++;
if(F[x] == x) F[x] = l;
if(F[y] == y) F[y] = l;
x = Pre[y];
}
}
int Match(int x)
{
for(int i = 1; i <= n; i++) F[i] = i;
memset(S, -1, sizeof S);
S[*(Top = Q) = x] = 0, Top++;
for(int *i = Q; i != Top; *i++)
for(int T = Head[*i]; T; T = Next[T])
{
int g = Go[T];
if(S[g] == -1)
{
Pre[g] = *i, S[g] = 1;
if(!Nxt[g])
{
for(int u = g, v = *i, lst; v; u = lst, v = Pre[u])
lst = Nxt[v], Nxt[v] = u, Nxt[u] = v;
return 1;
}
S[*Top = Nxt[g]] = 0, Top++;
}
else if(!S[g] && find(g) != find(*i))
{
int l = lca(g, *i);
blossom(g, *i, l);
blossom(*i, g, l);
}
}
return 0;
}
inline void Read(int& x)
{
x = 0;
while(*inpch < '0') *inpch++;
while(*inpch >= '0') x = x * 10 + *inpch++ - '0';
}
int main()
{
fread(inp, 1, 33554432, stdin);
Read(n), Read(m);
for(int i = 1; i <= m; i++)
{
Read(x), Read(y);
addedge(x, y);
addedge(y, x);
}
int ans = 0;
for(int i = n; i >= 1; i--)
if(!Nxt[i]) ans += Match(i);
printf("%d\n", ans);
for(int i = 1; i <= n; i++) printf("%d ", Nxt[i]);
putchar('\n');
return 0;
}
3.2.4 一般图最大权匹配(完全不懂怎么改的板子)
带花树算法 \(2 \leq|V| \leq 400, 1 \leq |E| \leq 79800\) 无法分析复杂度。
第一行输出总共的权值和。
接下来输出一般图最大权匹配中每个点匹配后对应的另一个点。
#include <bits/stdc++.h>
#define pb emplace_back
using namespace std;
bool chkmax(int &a,int b){ return a<b?a=b,true:false; }
bool chkmin(int &a,int b){ return b<a?a=b,true:false; }
const int N=605,inf=1<<30,NUL=1,IN=2,OUT=3; // 1.5n,2w
struct solver {
struct edge_t {
int u,v,w;
edge_t(){ u=v=w=0; }
edge_t(int u,int v,int w):u(u),v(v),w(w){}
} e[N][N];
int nc,n,m,h[N],lk[N],fl[N],fc,match[N],in[N],mpos[N],col[N],from[N];
vector<int> flower[N];
void init(int n0){
nc=n0,n=m=n0+(n0&1);
for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) e[i][j]=edge_t(i,j,0);
}
void add(int u,int v,int c){ if(chkmax(e[u][v].w,c)) e[v][u]=e[u][v]; }
int len(edge_t e){ return h[e.u]+h[e.v]-e.w*2; }
int other(edge_t e,int x){ return in[e.u]+in[e.v]-x; }
int slack(int x){ return mpos[x]?len(e[x][in[mpos[x]]])/(1+(col[x]==OUT)):inf; }
void upd(int v,int u){ if(!mpos[u]||len(e[u][mpos[u]])>len(e[u][v])) mpos[u]=v; }
void set(int x,int b){ in[x]=b; if(x>n) for(int y:flower[x]) set(y,b); }
int lca(int u,int v){
for(u=in[u],v=in[v],++fc;u||v;u=in[from[u]],v=in[from[v]]){
if(fl[u]==fc) return u; if(u) fl[u]=fc;
if(fl[v]==fc) return v; if(v) fl[v]=fc;
}
return 0;
}
void color(int x,int c){
col[x]=c,mpos[x]=0;
if(c!=IN){
for(int y=1;y<=m;++y) if(y!=x&&in[y]==y&&col[y]!=IN){
if(col[y]==OUT) upd(y,x);
if(c==OUT) upd(x,y);
}
}
}
void contract(int u,int v,int p){
int x=n+1; while(in[x]) ++x; chkmax(m,x);
match[from[x]=match[x]=match[p]]=in[x]=x;
vector<int> &c=flower[x];
for(u=in[u];u!=p;u=in[from[u]]) c.pb(u);
c.pb(p),reverse(c.begin(),c.end());
for(v=in[v];v!=p;v=in[from[v]]) c.pb(v);
for(int y:c) for(int k=1;k<=m;++k) if(!e[x][k].u||len(e[y][k])<len(e[x][k])) e[x][k]=e[k][x]=e[y][k];
set(x,x); for(int y:c) mpos[y]=from[y]=match[y]=col[y]=0; color(x,OUT);
}
void decompose(int b){
vector<int> &c=flower[b];
for(int x:c) set(x,x);
int len=c.size(),p,q; vector<int> cyc(len);
int fr=other(e[in[from[b]]][b],in[from[b]]),bk=other(e[b][match[b]],match[b]);
for(p=0;c[p]!=fr;++p); for(q=0;q<len;++q) cyc[q]=c[(p+q)%len]; for(p=0;cyc[p]!=bk;++p);
if(p&1) reverse(cyc.begin()+1,cyc.end()),p=len-p; from[cyc[0]]=from[b];
for(q=0;q<p;++q) color(cyc[q],(q&1)?OUT:IN),match[cyc[q]]=cyc[q^1],from[cyc[q+1]]=cyc[q];
for(q=p+1;q<len;++q) color(cyc[q],NUL),match[cyc[((q-1)^1)+1]]=cyc[q];
from[match[b]]=match[match[cyc[p]]=match[b]]=cyc[p],color(cyc[p],IN);
match[b]=in[b]=mpos[b]=from[b]=match[b]=0; c.clear();
for(int k=1;k<=m;++k) e[b][k]=e[k][b]=edge_t();
}
bool detect(int u,int v){
u=in[u],v=in[v];
vector<int> path;
if(col[v]==OUT){
int p=lca(u,v);
if(!p){
while(u) path.pb(u),u=in[from[u]];
reverse(path.begin(),path.end());
while(v) path.pb(v),v=in[from[v]];
}else return contract(u,v,p),false;
}else if(match[v]){
return from[v]=u,color(v,IN),from[match[v]]=v,color(match[v],OUT),false;
}else for(path.pb(v);u;u=in[from[u]]) path.pb(u);
for(int i=0;i<path.size();++i) match[path[i]]=path[i^1];
return true;
}
bool augment(){
for(int x=1;x<=m;++x) mpos[x]=from[x]=0,col[x]=in[x]?NUL:0;
for(int x=1;x<=m;++x) if(in[x]==x) color(x,match[x]?NUL:OUT);
while(true){
for(int x=1,c=1;c<=m;x=x%m+1,++c) if(in[x]==x&&col[x]!=IN&&!slack(x)) if(c=detect(mpos[x],x)) return true;
for(int x=n+1,c=1;c<=m-n;x=(x-n)%(m-n)+n+1,++c) if(in[x]==x&&col[x]==IN&&!h[x]) decompose(x),c=0;
int d=inf;
for(int x=1;x<=m;++x) if(in[x]==x) {
if(col[x]!=IN) chkmin(d,slack(x));
else if(x>n) chkmin(d,h[x]/2);
}
if(d==inf) return false;
for(int x=1;x<=n;++x){
if(col[in[x]]==OUT) h[x]-=d;
if(col[in[x]]==IN) h[x]+=d;
}
for(int x=n+1;x<=m;++x) if(in[x]==x) {
if(col[x]==OUT) h[x]+=d*2;
if(col[x]==IN) h[x]-=d*2;
}
}
return false;
}
void untie(vector<int>&c,int v){
int p=0,l=c.size();
for(int x:c) set(x,x);
while(c[p]!=in[v]) ++p;
if(in[v]>n) untie(flower[in[v]],v);
for(int i=1;i<l;i+=2) untie(c[(p+i)%l],c[(p+i+1)%l]);
c.clear();
}
void untie(int x,int y){
int u=e[x][y].u,v=e[x][y].v;
if(in[u]!=x) swap(u,v);
if(e[x][y].w) lk[lk[v]=u]=v;
if(x>n) untie(flower[x],u);
if(y>n) untie(flower[y],v);
}
void mwpm(){
for(int u=1;u<=n;++u){
for(int v=1;v<=n;++v) chkmax(h[u],e[u][v].w);
h[u]+=h[u]&1,in[u]=u;
}
while(augment());
for(int x=1;x<=m;++x) if(in[x]==x&&match[x]<x) untie(match[x],x);
long long ret=0;
for(int x=1;x<=n;++x) if(x<lk[x]) ret+=e[x][lk[x]].w;
cout<<ret<<'\n';
for(int x=1;x<=nc;++x) cout<<lk[x]<<" \n"[x==nc];
}
} ;
int main(){
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
solver a;
int n, m;
cin >> n >> m;
a.init(n);
for (int i = 1; i <= m; ++i) {
int u, v, c;
cin >> u >> v >> c;
a.add(u, v, c);
}
a.mwpm();
}
3.3 网络流
3.3.1 Dinic（未完成）
3.3.2 最小费用最大流
常规做法请阅读蓝书。
输出两个数，第一个数是最大流，第二个数是最小费用。
3.3.2.1 基于 Capacity Scaling 的弱多项式复杂度 dijkstra 做法
\(O(m^2 \log U \log m)\), \(U\) 为边的最大容量。
#include <iostream>
#include <queue>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long ll;
typedef pair<ll, int> pli;
const ll INF = 1e18;
const ll LARGE = 1e12;
int n, m;
vector<bool> vis;
vector<int> head, nxt, from, to, pre;
vector<ll> raw_cap, cap, cost, p, dis;
priority_queue<pli, vector<pli>, greater<pli> > q;
void add(int u, int v, ll f, ll w)
{
nxt.push_back(head[u]);
head[u] = to.size();
from.push_back(u);
to.push_back(v);
raw_cap.push_back(f);
cap.push_back(0);
cost.push_back(w);
}
void add_edge(int u, int v, ll f, ll w)
{
add(u, v, f, w);
add(v, u, 0, -w);
}
ll c(int id)
{
return p[from[id]] + cost[id] - p[to[id]];
}
void dijkstra(int s)
{
vis.assign(n + 2, false);
dis.assign(n + 2, INF);
pre.assign(n + 2, -1);
dis[s] = 0;
q.push(pli(0, s));
while (!q.empty())
{
int u = q.top().second;
ll w = q.top().first;
q.pop();
if (vis[u]) continue;
vis[u] = true;
for (int i = head[u]; ~i; i = nxt[i])
{
int v = to[i];
if (cap[i] && dis[v] > w + c(i))
{
dis[v] = w + c(i);
pre[v] = i;
q.push(pli(dis[v], v));
}
}
}
}
void add_one_cap(int id)
{
int u = from[id];
int v = to[id];
if (cap[id])
{
++cap[id];
return;
}
dijkstra(v);
if (dis[u] < INF && dis[u] + c(id) < 0)
{
++cap[id ^ 1];
while (u != v)
{
int x = pre[u];
--cap[x];
++cap[x ^ 1];
u = from[x];
}
}
else ++cap[id];
ll max_dis = 0;
ll cur_len = c(id);
for (int i = 1; i <= n; ++i) if (dis[i] < INF) max_dis = max(max_dis, dis[i]);
for (int i = 1; i <= n; ++i) p[i] += dis[i] < INF ? dis[i] : max_dis + max(0ll, -cur_len);
dijkstra(n + 1);
for (int i = 1; i <= n; ++i) p[i] += dis[i];
}
int main()
{
int s, t;
cin >> n >> m >> s >> t;
head.resize(n + 2, -1);
p.resize(n + 2, 0);
for (int i = 1; i <= m; ++i)
{
ll u, v, f, w;
cin >> u >> v >> f >> w;
add_edge(u, v, f, w);
}
add_edge(t, s, LARGE, -LARGE);
for (int i = 1; i <= n; ++i)
{
add_edge(n + 1, i, 0, 0);
cap[to.size() - 2] = 1;
}
for (int i = 40; i >= 0; --i)
{
for (int j = 0; j <= m * 2 + 1; ++j) cap[j] <<= 1;
for (int j = 0; j <= m * 2; j += 2)
{
if ((raw_cap[j] >> i) & 1)
{
add_one_cap(j);
}
}
}
ll min_cost = 0;
for (int i = 0; i < m; ++i) min_cost += cap[i << 1 | 1] * cost[i << 1];
cout << cap[m << 1 | 1] << ' ' << min_cost;
return 0;
}
3.3.2.2 基于基于 Capacity Scaling 的弱多项式复杂度 SPFA 做法
复杂度为 \(O(nm^2 \log U)\) ，但很难卡满。
#include <iostream>
#include <queue>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long ll;
typedef pair<ll, int> pli;
const ll INF = 1e18;
const ll LARGE = 1e12;
int n, m;
queue<int> q;
vector<bool> inq;
vector<ll> raw_cap, cap, cost, dis;
vector<int> head, nxt, from, to, pre;
void add(int u, int v, ll f, ll w)
{
nxt.push_back(head[u]);
head[u] = to.size();
from.push_back(u);
to.push_back(v);
raw_cap.push_back(f);
cap.push_back(0);
cost.push_back(w);
}
void add_edge(int u, int v, ll f, ll w)
{
add(u, v, f, w);
add(v, u, 0, -w);
}
void spfa(int s)
{
inq.assign(n + 1, false);
dis.assign(n + 1, INF);
pre.assign(n + 1, -1);
dis[s] = 0;
q.push(s);
while (!q.empty())
{
int u = q.front();
inq[u] = false;
q.pop();
for (int i = head[u]; ~i; i = nxt[i])
{
int v = to[i];
ll w = cost[i];
if (cap[i] && dis[v] > dis[u] + w)
{
dis[v] = dis[u] + w;
pre[v] = i;
if (!inq[v])
{
inq[v] = true;
q.push(v);
}
}
}
}
}
void add_one_cap(int id)
{
int u = from[id];
int v = to[id];
if (cap[id])
{
++cap[id];
return;
}
spfa(v);
if (dis[u] < INF && dis[u] + cost[id] < 0)
{
++cap[id ^ 1];
while (u != v)
{
int x = pre[u];
--cap[x];
++cap[x ^ 1];
u = from[x];
}
}
else ++cap[id];
}
int main()
{
int s, t;
cin >> n >> m >> s >> t;
head.resize(n + 1, -1);
for (int i = 1; i <= m; ++i)
{
ll u, v, f, w;
cin >> u >> v >> f >> w;
add_edge(u, v, f, w);
}
add_edge(t, s, LARGE, -LARGE);
for (int i = 40; i >= 0; --i)
{
for (int j = 0; j <= m * 2 + 1; ++j) cap[j] <<= 1;
for (int j = 0; j <= m * 2; j += 2)
{
if ((raw_cap[j] >> i) & 1)
{
add_one_cap(j);
}
}
}
ll min_cost = 0;
for (int i = 0; i < m; ++i) min_cost += cap[i << 1 | 1] * cost[i << 1];
cout << cap[m << 1 | 1] << ' ' << min_cost;
return 0;
}
4.字符串
4.1 Manacher 算法
p[i]-1 是以 i 为中心的回文串长度。
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const int N = 11000000;
char c[N+5],s[2*N+5];
int p[2*N+5];
int Manacher(char *s) {
memset(p,0,sizeof(p));
int max_len=0,mid=0,max_r=0,len=strlen(s);
for(int i=0;i<len;i++) {
p[i]=min(p[2*mid-i],max_r-i);
while(s[i+p[i]]==s[i-p[i]]) p[i]++;
if(max_r<i+p[i]) max_r=i+p[i],mid=i;
max_len=max(max_len,p[i]-1);
}
return max_len;
}
int main() {
while(scanf("%s",c)!=EOF) {
int len=strlen(c);
s[0]='$',s[1]='#';
for(int i=0;i<len;i++)
s[2*i+2]=c[i],s[2*i+3]='#';
s[2*len+2]='\0';
printf("%d\n",Manacher(s));
}
return 0;
}
5.计算几何
5.1 极角排序
此题 SWERC2018 F.Paris by Night 需要计算一条过两个给定点的直线，使得直线分开的两边权值之和的差值最小。
这里的极角排序是利用判断象限后求叉积，以辐角 \([0,2 \pi )\) 从小到大排序。
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
#define IL inline
#define ri register int
const db eps = 1e-10;
IL int dcmp(db x) {if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;}
struct P {
db x,y;
IL P(db x=0.0, db y=0.0, db w=0):x(x),y(y) {}
};
typedef P V;
IL V operator + (const V& a,const V& b) {return V(a.x+b.x, a.y+b.y);}
IL V operator - (const V& a,const V& b) {return V(a.x-b.x, a.y-b.y);}
IL V operator * (const V& a,db p) {return V(a.x*p,a.y*p);}
IL V operator / (const V& a,db p) {return V(a.x/p,a.y/p);}
IL bool operator < (const P& a, const P& b) {return a.x < b.x || (a.x == b.x && a.y < b.y);}
IL bool operator == (const P& a, const P& b) {return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0;}
IL db Dot(const V& a,const V& b) {return a.x*b.x+a.y*b.y;}
IL db Length(const V& a) {return sqrt(Dot(a,a));}
IL db Angle(const V& a,const V& b) {return acos(Dot(a,b) / Length(a) / Length(b));}
IL db Cross(const V& a,const V& b) {return a.x*b.y-a.y*b.x;}
IL int Quadrant(const P& a) {
if(dcmp(a.x) > 0 && dcmp(a.y) >= 0) return 1;
if(dcmp(a.x) <= 0 && dcmp(a.y) > 0) return 2;
if(dcmp(a.x) < 0 && dcmp(a.y) <= 0) return 3;
if(dcmp(a.x) >= 0 && dcmp(a.y) < 0) return 4;
return 0;
}
IL bool cmp(const P& a, const P& b) {
return Quadrant(a) < Quadrant(b) || (Quadrant(a) == Quadrant(b) && dcmp(Cross(a,b)) > 0);
}
const int N = 4000 + 5;
const LL INFL = 1e18 + 5;
struct PP {
P p; int w;
PP (P p=P(0,0), int w=0):p(p),w(w){}
};
PP p[N], tmp[N];
IL bool OnLeft(const P& a, const V& vec) { return dcmp(Cross(vec, a)) > 0;}
LL sum = 0, suml = 0, sumr = 0, ans = INFL;
int main() {
sum = suml = sumr = 0;
ans = INFL;
int n; scanf("%d",&n);
for(ri i=0;i<n;i++) {
int x,y,w; scanf("%d%d%d",&x,&y,&w);
p[i] = PP(P(x,y),w);
sum += w;
}
for(ri i=0;i<n;i++) {
int m = 0;
suml = sumr = 0;
for(ri j=0;j<n;j++) {
if(i == j) continue;
tmp[m++] = PP(p[j].p-p[i].p, p[j].w);
}
sort(tmp, tmp+m, [](PP a, PP b) {return cmp(a.p, b.p);});
for(ri l=0,r=1;l<m;l++) {
for( ;OnLeft(tmp[r].p, tmp[l].p); ) { // the point of (l,r) are on left.
suml += tmp[r].w;
r = (r+1) % m;
}
sumr = sum - p[i].w - tmp[l].w - suml;
ans = min(ans, abs(sumr-suml));
if(l != m-1) suml -= tmp[l+1].w;
}
}
printf("%lld\n",ans);
return 0;
}
6.动态规划
7.杂项
7.1 卡常火车头
#define fastcall __attribute__((optimize("-O3")))
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
7.2 防止 unordered_map 被卡方法
struct custom_hash {
static uint64_t splitmix64(uint64_t x) {
x += 0x9e3779b97f4a7c15;
x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
return x ^ (x >> 31);
}
size_t operator()(uint64_t x) const {
static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x + FIXED_RANDOM);
}
};
unordered_map<long long, int, custom_hash> safe_map;
gp_hash_table<long long, int, custom_hash> safe_hash_table;
7.3 emacs
7.3.1 emacs 配置
(global-linum-mode t) ;;line num 显示行号
(setq-default indent-tabs-mode nil) ;;将tab转变为空格
(setq c-basic-offset 4)
(setq default-tab-width 4) ;;tab width 设置缩进以及tab键
(electric-pair-mode t)
(electric-layout-mode t)
(electric-indent-mode t)
;;electric pair 括号补全
(show-paren-mode t) ;;show paren 括号配对
;;(global-hl-line-mode t) ;;high light line 高亮当前行
(global-set-key [f6] 'gdb) ;;摁F6进入gdb调试
(global-set-key [f7] 'compile) ;;摁F7编译
(global-set-key [f8] 'shell) ;;摁F8进入shell
;;windows下的含中文的文件emacs过于慢的问题。
(dolist (charset '(kana han symbol cjk-misc bopomofo))
(set-fontset-font (frame-parameter nil 'font)
charset
(font-spec :family "Microsoft Yahei" :size 18)))
