前言补于 2025-6-14。然而,不幸的是,时间不足够,故开始对于数学笔记的补充。
好吧,这篇笔记会不定期更新,同步在洛谷博客和一位朋友 RyanLi 先生的博客之中。在此感谢他。
我们补充符号 $[p]$,该符号的意义是:若 $p$ 是真的,则 $[p]=1$;否则 $[p]=0$。
数论
在这个大篇的内容中,数都是整数,不论其为正数,负数或 $0$。如有不同,另会说明。
整除
若 $a=kb$,则 $b\backslash a$,即 $b$ 整除 $a$。它的性质显而易见,故不言明。这个符号是高德纳的,因为用竖线会把它和集合混一起。
同余
若 $a=mk+r,r\in(-m,m)\cup\mathbb Z$,我们可以得出两个式子:
- 取余:$a\bmod m=r$。
- 同余:$a\equiv r\pmod m$。
它们都很重要。
素数筛法
定义素数,即只能被 $1$ 或自身整除的自然数(且不是 $0$ 或者 $1$!),且所有不是素数的自然数我们都叫做合数。
我们有一个朴素的筛质数的方式:首先将 $2\sim n$ 的所有整数按顺序排列在纸上,首先划掉 $2$ 的所有倍数,随后划掉所有 $3$ 的所有倍数,如果某个数已经被划掉,那么我们不再管它,随后划掉 $5$ 的所有倍数……
这个东西其实叫做埃拉托斯克尼筛法,代码如下:
memset(vis,0,sizeof vis);
for(int i=2;i<=n;++i)
for(int j=i*2;j<=n;j+=i) vis[j]=1;
但是,这个算法有点慢,因为它是 $O(n\log n)$ 的(证明详见调和级数)。所以,我们改进一下就得到
memset(vis,0,sizeof vis);
pc=0;
for(int i=2;i<=n;++i){
if(vis[i]==0) pr[++pc]=i;
for(int j=1;j<=pc&&i*pr[j]<=n;++j){
vis[i*pr[j]]=0;
if(i%pr[j]==0) break;
}
}
这个算法是欧拉的线性筛法,复杂度几乎是 $O(n)$,足够使用。
唯一分解定理 是指对于某数 $n$,其可以分解成唯一的质因数分解式 $n=\prod_{p\in\mathbb P}p^{e_p}$。
最大公约数&最小公倍数
定义两数的 最大公约数 $a,b$ 为可以同时整除两数的最大数 $c$,形式化是 $\gcd(a,b)=(a,b)=\max\{c|(c\backslash a)\land(c\backslash b)\}$。
定义两数的 最小公倍数 $a,b$ 为两数可以整除该数 $c$ 的最大值,形式化式 $\mathrm{lcm}(a,b)=[a,b]=\min\{c|(a\backslash c)\land(b\backslash c)\}$。
不难通过唯一分解定理证明,$(a,b)[a,b]=ab$。
这里需介绍最大公约数的求法:$(a,b)=(b,a\bmod b)$。
证明。令 $a\bmod b=r$, 所以 $a=kb+r$, 所以 $r=a-kb$。假设 $g$ 为 $a,b$ 的一个公约数,那么 $g\backslash a,g\backslash b$,所以 $g\backslash (a-kb)$,即 $g\backslash r$,也就是 $g\backslash(a\bmod b)$。所以 $g$ 是 $a,b,a\bmod b$ 的公约数。证毕(本法被称为欧几里得算法)。
当 $b$ 到 $0$ 时,这个操作无法进行,得到的 $a$ 即最大公约数。这个方法我们可以把它压到一行里:
int gcd(int a,int b){
return b==0?a:gcd(b,a%b);
}
二元一次不定方程
众所周知,在实数范围内,$ax+by=c$ 一定有解。可是这是数论。故我们需要一些更实在的东西。
裴蜀定理是指:对于非零整数 $a,b$,满足对于任意整数 $x,y$,$(a,b)\backslash ax+by$,同时存在 $(x,y)$ 满足方程 $ax+by=(a,b)$。证明略。
我们可以导出一些重要结论。
众所周知,$(a,b)=(b,a\bmod b)$。所以,我们可以进行换元,把 $(a,b)\to(b,a\bmod b)$,这样我们得出两条方程:
$$ax+by=(a,b)$$
$$bx+(a\bmod b)y=(a,b)$$
消掉右边,得到 $ax+by=bx+(a\bmod b)y$。
众所周知,我们可以通过下取整除法得到余数,即
$$\mathrm{RHS}=bx+\left(a-b\left\lfloor\frac ab\right\rfloor\right)y$$
简单代换,得到
$$\mathrm{RHS}=ay+b\left(x-\left\lfloor\frac ab\right\rfloor y\right)$$
这样,我们就可以通过变量代换 $(x,y)\to(y,x-\lfloor a/b\rfloor y)$ 得到解。解的最下一层是 $(a,0)$,它的解是 $x=1,y=0$,然后递归即可。
代码如下。
void exgcd(int a,int b,int &x,int &y){
if(b==0) x=1,y=0;
else exgcd(b,a%b,y,x),y-=a/b*x;
}
如果解 $ax+by=c$ 时带入即可,当然,如果 $\lnot (a,b)\backslash c$,则无解。
一元模方程组
考虑如下方程组(满足 $a$ 两两互质)
$$ \begin{cases} x\equiv r_1\pmod{a_1}\\ x\equiv r_2\pmod{a_2}\\ \cdots\\ x\equiv r_n\pmod{a_n}\\ \end{cases} $$
我们要得出 $x$,就需要求助老智慧。
- 我们令 $A=\prod a_k$,并令 $\alpha_k=A/a_k$。
- 求出满足 $\alpha_k\alpha_k^{-1}\equiv 1\pmod{a_k}$。
- 设 $q_k=\alpha_k\alpha_k^{-1}$。没有取模。
则答案即为 $$\left(\sum_{k=1}^n r_kq_k\right)\bmod A$$。容易证明,对于任意答案 $C$,$C+kA$ 同样满足解。
我们来分析一下解的过程。这里我们只考虑某一个 $k$ 的情况。
对于方程组
$$ \begin{cases} x\equiv 1\pmod{a_k}\\ x\equiv 0\pmod{a_h}&h\neq k \end{cases} $$
容易证明, $x$ 的解是所有 $a_h$ 的乘积的常数倍,这是因为后一条方程。如果我们要同时满足前一条方程,我们需要在它的左端乘上 $\alpha_k=A/a_k$。即
$$x\equiv \alpha_k t\equiv1\pmod{a_k}$$
这相当于找出 $\alpha_k$ 对 $a_k$ 的逆元。当然,如果我们两端乘以 $r_k$ 时,等式依然成立。
现在我们删除后面的同余记号,并把得到的解累加起来,就满足了所有的方程,这是因为各个项满足了各个方程,故这个解满足所有方程。
代码如下。
int CRT(int n,int a[],int r[]){
int A=1,ans=0;
for(int i=1;i<=n;++i) A*=a[i];
for(int i=1;i<=n;++i){
int ah=A/a[i],ahi,tmp,g=gcd(ah,a[i]);
exgcd(ah,a[i],ahi,tmp);
ans+=r[i]*(ah*ahi)+A,ans%=A;
}
return (ans+A)%A;
}
当然,如果它们不互质,你就死定了。
不互质的情况
很不幸,不互质的时候,exCRT 事实上就是合并方程组。我不想这么做。所以我参考了一些别的书,写了下列的算法。
众所周知,对于质数 $p$,满足 $a\equiv n\pmod {p^k},a\equiv m\pmod{p^h},k<h$,那么当 $n\equiv m\pmod{p^k}$ 时,我们的要求即告满足,即 $m$ 为我们所要求的。
现在对于所有的 $a$ 进行唯一分解,然后将余数对 $p^k$ 取模,其中 $p$ 是其所有的质因子,$p^k$ 满足 $k$ 最大且可以被 $a$ 整除。按照上面的逻辑,保存对于 $p$ 的最大的 $p^k$ 以及最大的模数,最后对于这些模数如上所述的操作(如果不满足条件则直接无解),留下一份 $p^k$ 与模数的对应表。把表扔进上文所述的算法,这就是另类的 exCRT 了。
遗憾的是,它的时间复杂度是 $O(n\sqrt{N})$ 的,而 $N$ 在洛谷的板题很大,故而跑不过去。不过这个东西对于小数据应当够用。
和下文的多项式工业一样,不要滥用模版,因为它们不配套。
ll exCRT(vector<pair<ll,ll>> ar){
map<ll,ll> akh,pkh;
for(auto [a,r]:ar){
ll A=a,p=1;
while(A!=1&&p*p<=A){
++p;
ll pk=1;
while(A%p==0) pk*=p,A/=p;
if(pk==1) continue;
if(akh.count(p)){
if(pkh[p]>=pk&&akh[p]%pk!=r%pk) return -1;
else if(pkh[p]<pk&&r%pkh[p]!=akh[p]) return -1;
else if(pkh[p]<pk) pkh[p]=pk,akh[p]=r%pk;
}else pkh[p]=pk,akh[p]=r%pk;
}
if(A!=1){
if(akh.count(A)){
if(pkh[A]>=A&&akh[A]%A!=r%A) return -1;
}else pkh[A]=A,akh[A]=r%A;
}
}
vector<pair<ll,ll>> wyp;
for(auto ie:akh) wyp.push_back({pkh[ie.first],ie.second});
return CRT(wyp);
}
数论分块
现在,假设我们要对 $i=1\dots n$ 求函数 $f(i)\lfloor n/i\rfloor$ 的和。现在,我们考虑一下,这种 $\lfloor n/i\rfloor$ 的函数图像必然是一条直线。所以,可以考虑先求出 $f$ 的前缀和 $\sum f$,而后,我们可以搬上结论:
对于 $n$ 和整数 $i$,使得 $\lfloor n/i\rfloor=\lfloor n/j\rfloor$ 的最小 $j$ 为 $j=\lfloor n/(\lfloor n/i\rfloor)\rfloor$。
证明请参考 OI-wiki。代码大致是这样。
void NTsum(ll n){
ll ans=0,i=1,j;
while(i<=n){
j=min(n/(n/i),n);
ans+=(F[j]-F[i-1])*(n/i);
i=j+1;
}
return ans;
}
可以通过端点的减法来获得我们的解。
一些重要的函数及其性质
先介绍互质符号:$a\perp b\equiv(\gcd(a,b)=1)$。
首先,我们先介绍积性函数的概念:若 $a\perp b$,且若 $f(ab)=f(a)f(b)$,则 $f$ 被称为积性函数。
这里首先介绍一个性质:若函数 $f,g$ 满足
$$g(n)=\sum_{1\le k\le n}{f(n)}$$
且 $f$ 是积性函数,则 $g$ 也是积性函数。
欧拉函数
我们首先介绍 $\varphi$ 的概念:
$$\varphi(n)=\sum_{k=1}^n[k\perp n]$$
我们首先探究这个函数在质数及质数幂处的取值:$\varphi(p)=p-1$,这很显然;与此同时,对于 $\varphi(p^k)$,显然有 $p^{k-1}$ 个与之不互素的数字,则 $\varphi(p^k)=p^k-p^{k-1}$。所以,按照积性函数,我们可以得到
$$\varphi(n)=\prod_{p\backslash n}{p^{k-1}(p-1)}=n\prod_{p\backslash n}{\frac{p-1}{p}}=n\prod_{p\backslash n}{\left(1-\frac1p\right)}$$
当然这是单点处的取值。我们可以通过筛法来预处理小于等于 $n$ 的全部 $\varphi$。
int n,pr[N],vphi[N],tot;
bool vis[N];
void init(){
memset(vis,1,sizeof(vis));
vis[0]=vis[1]=0;
for(int i=2;i<=n;++i){
if(vis[i]) pr[++tot]=i;
for(int j=1;j<=tot&&i*pr[j]<=n;++j){
vis[i*pr[j]]=0;
if(i%pr[j]==0){
vphi[i*pr[j]]=i*pr[j];
break;
}else vphi[i*pr[j]]=i*(pr[j]-1);
}
}
}
现在,我们可以得到:
$$n=\sum_{d\backslash n}{\varphi(d)}$$
这个可以通过下一节中介绍的函数得出。
莫比乌斯函数
现在,我们要介绍一个取值范围较小的函数,即莫比乌斯函数。它满足这样一个性质:
$$\sum_{d\backslash n}\mu(d)=[n=1]$$
因而,我们可以得出这个 $\mu$ 函数的全体值,这是因为按照上文所述的 $[n=1]$ 也是一个积性函数。故而,我们将 $p^k$ 代入 $n$,得到
$$\mu(1)+\mu(p)+\mu(p^2)+\cdots+\mu(p^k)=0$$
因为 $\mu(1)=1$,且当 $k=1$ 时仅有 $\mu(p)+1=0$,故而 $\mu(p^k)=-[k=1]$。
我们补一下 $\mu$ 的欧拉筛预处理。
int tot,pr[N],n=1e6,mu[N],M[N];
bool vis[N];
void init(){
memset(vis,1,sizeof(vis));
vis[0]=vis[1]=0,mu[1]=M[1]=1;
for(int i=2;i<=n;++i){
if(vis[i]) pr[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*pr[j]<=n;++j){
vis[i*pr[j]]=0;
if(i%pr[j]==0) break;
mu[i*pr[j]]=-mu[i];
}
M[i]=M[i-1]+mu[i];
}
}
我们给出一个重要的式子,即莫比乌斯反演式:
$$g(m)=\sum_{d\backslash m}{f(d)}\Leftrightarrow f(m)=\sum_{d\backslash m}\mu\left(\frac md\right)g(d)$$
该式子可以让我们可以根据某个函数的取值来推测它的因数中各个因数的取值。
函数们的性质
首先,对于欧拉和莫比乌斯函数,有如下的两条值得注意的性质:
$$\sum_{d\backslash n}{\varphi(d)}=n$$
$$\sum_{d\backslash n}{\mu(d)}=[n=1]$$
这两条性质极端的重要,都可以让我们枚举因数并提到整个计算式的前方;第一个式子用来计算 $\gcd$,第二个用来计算我们能想到的很多东西。
技巧性内容
我们可以在这里补充几个技巧性的内容。
按照刚刚所述的莫比乌斯函数,我们知道,按照莫比乌斯函数的理论,我们可以简化下列内容的求和:
$$\sum_{i=1}^n\sum_{j=1}^m(\text{A function related to $\gcd(i,j)$})$$
我们以 $S=\sum_{i=1}^n\sum_{j=1}^m{\gcd(i,j)}$ 为例子。我们可以把它改写作
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{k}k[\gcd(i,j)=k]$$
我们可以将这个式子视作一个乘法,把 $\sum_k k$ 提到最前面:
$$=\sum_{k}k\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]$$
有一个经典的结论,它允许我们把关于 $\gcd$ 的求和指标都除以 $k$。
$$=\sum_kk\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor}[\gcd(i,j)=1]$$
这样的话,我们可以用经典的 $d$ 来做一个变换了,因为 $[n=1]=\sum_{d\backslash n}\mu(d)$。
$$=\sum_kk\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor}\sum_d\mu(d)[d\backslash i][d\backslash j]$$
再次应用“乘法”技巧,把关于 $d$ 的部分式子提到前面:
$$=\sum_kk\sum_d\mu(d)\sum_{i=1}^{\lfloor n/k\rfloor}[d\backslash i]\sum_{j=1}^{\lfloor m/k\rfloor}[d\backslash j]$$
好的,我们可以直接应用一个著名的方法,于是我们的和式就变成了这样:
$$=\sum_kk\sum_d\mu(d)\lfloor n/kd\rfloor\lfloor m/kd\rfloor$$
现在我们再来应用一次可爱的枚举法:我们令 $\vartheta=kd$,则这个式子可以改写为
$$=\sum_\vartheta\lfloor n/\vartheta\rfloor\lfloor m/\vartheta\rfloor\sum_d\mu(d)\sum_kk[dk=\vartheta]$$
由于 $k$ 关于 $d$ 和 $\vartheta$ 是唯一的,故而我们可以删掉关于 $k$ 的求和式:
$$=\sum_\vartheta\lfloor n/\vartheta\rfloor\lfloor m/\vartheta\rfloor\sum_d\mu(d)(\vartheta/d)[d\backslash\vartheta]$$
我们当然可以把它并入求和指标之中:
$$=\sum_\vartheta\lfloor n/\vartheta\rfloor\lfloor m/\vartheta\rfloor\sum_{d\backslash\vartheta}\mu(d)(\vartheta/d)$$
第二个 $\sum$ 可以用 $\mathrm O(n\log n)$ 的复杂度预处理,而第一个 $\sum$ 可以使用数论分块的 $\mathrm O(\log n)$ 来求得。不过第二个实际上就是 $\varphi(\vartheta)$。所以我们的结论还可以缩得更短:
$$\sum_{\vartheta}\varphi(\vartheta)\left\lfloor\frac n\vartheta\right\rfloor\left\lfloor\frac m\vartheta\right\rfloor$$
这种技巧可以处理较为广阔的东西。
杜教筛
嗯,这里需要先介绍符号,因为没它没法简洁的推出来我们的结论。
定义狄利克雷卷积满足对于两个数论函数 $f,g$,它们的狄利克雷卷积是
$$h(n)=\sum_{d\backslash n}f(d)g\left(\frac nd\right)$$
我们记作 $h=f\ast g$。所以,
$$\sum_{d\backslash n}{\varphi(d)}=n$$
$$\sum_{d\backslash n}{\mu(d)}=1$$
可以被简化成
$$\varphi\ast1=\mathrm{id},\mu\ast1=\epsilon$$
其中 $\mathrm{id}(n)=n,\epsilon(n)=[n=1]$。
杜教筛用于求数论函数的前缀和,而复杂度是 $O(n^{2/3})$ 的,面对很大很大的数据时有优势。定义 $f$ 的前缀和 $S$ 满足
$$S(n)=\sum_{k=1}^n{f(k)}$$
对于任意的数论函数满足
$$\sum_{k=1}^n{(f\ast g)(k)}=\sum_{k=1}^n{g(k)S\left(\left\lfloor\frac nk\right\rfloor\right)}$$
证明可以使用上面的技巧。所以,我们的最终目的是:构造 $g$,使得 $f\ast g$ 的前缀和容易被找出,且 $g$ 的前缀和也容易被找出。这样,我们要求的 $S(n)$,也就等于
$$g(1)S(n)=\sum_{k=1}^n{g(k)S(\lfloor n/k\rfloor)}-\sum_{k=2}^n{g(k)S(\lfloor n/k\rfloor)}=\sum_{k=1}^n{(f\ast g)(k)}-\sum_{k=2}^n{g(k)S(\lfloor n/k\rfloor)}$$
所以,我们可以用整数分块和递归来求我们想要的函数。
具体的,对于 $\varphi$ 的话,我们当然使用上面的结论,令 $f=\varphi,g=1$,则算式的左边等于 $S_\varphi(n)$,右边等于 $n(n+1)/2-\sum_{k=2}^n{S_\varphi(\lfloor n/k\rfloor)}$,还可以用数论分块再简化。代码如下。
// 对于较小的 x,已经用线性筛做了预处理。
ll vPhi(ll x){
if(x<N) return vphi[x];
if(svphi[x]!=0) return svphi[x]; // 记忆化
ll ans=x*(x+1)/2; //式子的左端
for(ll lf=2,rt;lf<=x;lf=rt+1){
rt=x/(x/lf); // 数论分块
ans-=(ll)(rt-lf+1)*vPhi(x/lf); // 递归
}
return svphi[x]=ans;
}
费马小定理和欧拉定理
之所以将它们放在一块,因为可以直接解。
若 $a\perp n$,则 $a^{\varphi(n)}\equiv 1\pmod n$
就这样。
生成函数和多项式
多项式和生成函数理论非常有用。可惜以我的学识,我可能不能提供什么。
生成函数简论
生成函数一般用来处理组合计数的问题。……嗯,我之前写的东西都丢了,所以只能重写了。
我们假设一个无穷数列为 $f_0,f_1,\cdots$,简单的记为 $\langle f\rangle$,现在我们可以写出三种生成函数,分别是普通生成函数(Ordinary Generating Function),指数生成函数(Exponential Generating Function)和狄利克雷生成函数(Dirichlet Generating Function):
$$F(z)=\sum_{k}{f_kz^k}$$
$$\hat F(z)=\sum_{k}{\frac{f_k}{k!}z^k}$$
$$\tilde F(n)=\sum_{k\geq 1}{\frac{f_k}{k^n}}$$
OGF 适用于统计计数,EGF 适用于组合,而 DGF 是给数论用的。我们现在还没有法子来同时辨析这么多函数,因此我们先专注于 OGF。注意前两个函数我们使用了 $z$,这意味着这两个都是复变函数。第三个我们选用字母 $n$ 作为自变量,这意味着它是一个整数函数。
对于 OGF 和 EGF,运算技巧可以直接列表。这里我们注明:对应的小写字母数列直接对应其大写字母的生成函数。
$$\alpha F(z)+\beta G(z)=\sum_n{(\alpha f_n+\beta g_n)z^n}$$
$$z^mF(z)=\sum_n{g_{n-m}z^n}$$
$$F(cz)=\sum_n{(c^nf_n)z^n}$$
$$F'(z)=\sum_n{[(n+1)g_{n+1}]z^n}$$
$$\int_0^zF(t)\mathrm dt=\sum_{n\geq 1}{(f_{n-1}/n)z^n}$$
$$F(z)G(z)=\sum_{n}\left(\sum_k{f_kg_{n-k}}\right)z^n$$
$$\hat F(z)\hat G(z)=\sum_{n}\left(\sum_k{\binom{n}{k}f_kg_{n-k}}\right)\frac{z^n}{n!}$$
最后一个式子可以让我们对于两个数列进行卷积。
熟悉一些无穷级数的推导仍然是有效的;例如,我们可以快速地证明,数列 $1,1,1,1,1,1,\cdots$ 的 ogf 是 $1/(1-z)$。
为此,我们引入广义二项式定理:
$$(x+y)^r=\sum_k{\binom{r}{k}x^ky^{r-k}}$$
这样就能解释下面的列表中出现的一些式子了。假定诸位可以通过函数的系数来推断它的样子。
$$\sum_{n\geq 0}z^n=\frac1{1-z}$$
$$\sum_{n\geq 0}{(-1)^nz^n}=\frac1{1+z}$$
$$\sum_{n\geq 0}[m\backslash n]z^n=\frac1{1-z^m}$$
$$\sum_{n\geq 0}k^nz^n=\frac1{1-kz}$$
$$\sum_{n\geq 0}\binom{k}{n}z^n=(1+z)^k$$
$$\sum_{n\geq 0}\binom{k+n-1}{n}z^n=(1+z)^{-k}$$
$$\sum_{n\geq 0}\frac{z^n}{n!}=e^z$$
$$\sum_{n\geq 0}\binom{m+n}{m}z^n=(1-z)^{-(m+1)}$$
$$\sum_{n\geq 1}\frac{z^n}{n}=\ln\frac1{1-z}$$
$$\sum_{n\geq 1}\frac{(-1)^{n+1}z^n}{n}=\ln(1+z)$$
不幸,OI 中我们是见不着多项式的封闭形式的:我们只能见到各个的系数在 FDFT 里跳跃。关于什么是 FDFT,请阅读下一节的内容。我们现在来看看生成函数的具体应用。
组合的乘积
我们来观赏一道例题。
在许多不同种类的食物中选出 n 个,每种食物的限制如下:
- 承德汉堡:偶数个
- 可乐:0 个或 1 个
- 鸡腿:0 个,1 个或 2 个
- 蜜桃多:奇数个
- 鸡块:4 的倍数个
- 包子:0 个,1 个,2 个或 3 个
- 土豆片炒肉:不超过一个。
- 面包:3 的倍数个
每种食物都是以「个」为单位,只要总数加起来是 n 就算一种方案。对于给出的 n 你需要计算出方案数,对 10007 取模。
我们只需要按照顺序列出对应的生成函数就解完了这道题。注意个数实际上是通过某一个项的系数来决定的。
- 偶数个承德汉堡:$1/(1-z^2)$
- 不多于一个可乐:$1+z$
- 不多于两个鸡腿:$1+z+z^2=(1-z^3)/(1-z)$
- 奇数个蜜桃多:$z/(1-z^2)$
- 四的倍数个鸡块:$1/(1-z^4)$
- 不多于三个包子:$1+z+z^2+z^3=(1-z^4)/(1-z)$
- 不多于一个土豆片炒肉:$1+z$
- 三的倍数个面包:$1/(1-z^3)$
乘在一起,得到 $z/(1-z)^4$。
由
$$\sum_{n\geq 0}\binom{k+n-1}{n}z^n=(1+z)^{-k}$$
和
$$z^mF(z)=\sum_n{g_{n-m}z^n}$$
得到我们的数列是 $\binom{n+2}{n-1}=\binom{n+2}{3}$。这样我们就可以迅速地解决问题。请注意,本题中有一个特殊的点,即 $n\geq 10^{500}$。这需要我们边读入边取模。与此同时,我们当然可以使用定义秒了这个题。
多项式工业
多项式乘积
求两个多项式的乘积。
一般来说,我们可以使用暴力,不过这样的复杂度是 $\mathrm O(n^2)$ 的,有些过于慢了。
故,我们在这里介绍一种算法——FDFT,即所谓「快速离散傅里叶变换」,用于解决这个问题。
复数
FDFT 的基础是复数,在这里我们探讨两样东西:
- 欧拉公式。我们直接给出这个公式:
$$e^{\mathrm i\theta}=\cos\theta+\mathrm i\sin\theta$$
证明可以通过将泰勒展开拓展到复数域上从而实现。
这个公式可以快速地求我们下面介绍的这个东西:
- $\mathbf n$ 次单位根。我们定义
$$\omega_{n}=e^{2\pi\mathrm i/n}$$
这样我们就可以快速的在单位圆上转圈。它有如下的数个性质:
$$\omega_{hn}^{hk}=\omega_{n}^{k}$$
$$\omega_{n}^{k+n/2}=-\omega_{n}^k$$
$$\omega_{n}^{k+n}=\omega_{n}^k$$
$$\omega_n^{-k}=(\omega^k_n)^*$$
这些个性质可以自己试着证明,这里不再赘述。
注意单位根的性质实际上是在复平面上进行旋转。
多项式
一般来说,确定一个多项式有多种方法,这里仅介绍其中的两种。
- 系数表示法:对于一个次数为 $n-1$ 的多项式,系数 $a_0,a_1,\cdots,a_n$ 确定了它。
- 点值表示法:对于一个次数为 $n-1$ 的多项式,可以通过带入 $n$ 个特殊的数来确定一个多项式。
如果多项式的次数一定,我们可以在这两种方法转化,从前者转为后者的方式(在普通的方法)是「秦九韶算法」,允许我们用 $O(n^2)$ 的时间复杂度代入 $n$ 个数并得到它们的对应值;而从后者转为前者的方式则是解方程,但是时间复杂度可能更高。
但是,点值表达式允许我们快速地进行乘法,所以我们现在想办法优化它们间的转化。
FDFT
现在我们设 $f(x)=\sum_{k=0}^{2^n-1}{a_kx^k}$。我们把系数按照奇偶进行分组,得到 $g(x)=\sum_{k=0}^{2^n-1}{a_{2k}x^{k}}$,$h(x)=\sum_{k=0}^{2^n-1}{a_{2k+1}x^{k}}$,我们发现它们满足这样的关系式:
$$f(x)=g(x^2)+xh(x^2)$$
我们带入 $\omega_{2^n}^{k}$ 和 $\omega_{2^n}^{k+2^n}$ 就得到了
$$f(\omega_{2^n}^{k})=g(\omega_{2^n}^{2k})+\omega_{2^n}^{k}h(\omega_{2^n}^{2k})=g(\omega_{2^n}^{k})+\omega_{2^n}^{k}h(\omega_{2^n}^{k})$$
$$f(\omega_{2^n}^{k+2^n})=g(\omega_{2^n}^{2k+2^{n}})+\omega_{2^n}^{k+2^n}h(\omega_{2^n}^{2k+2^n})=g(\omega_{2^n}^{k})-\omega_{2^n}^{k}h(\omega_{2^n}^{k})$$
并且,按照肉眼的观察,我们发现 $g,h$ 的次数都是 $f$ 的一半,所以可以用递归实现从系数形式向点值形式的转化。
IDFT
我们现在想要从点值形式转回系数形式。按照线性代数的知识,我们可以将代值的形式写作这样:
$$\begin{bmatrix}x^0&x^1\cdots&x^{2^n-1}\end{bmatrix}\begin{bmatrix}a_0\\ a_1\\\cdots\\a_{2^n-1}\end{bmatrix}=f(x)$$
令 $x=\omega_{2^n}^{0},\omega_{2^n}^{1},\cdots,\omega_{2^n}^{2^n-1}$ 并把它们依次带入式子的左侧,把结果堆叠起来我们即得到了矩阵形式:
$$\begin{bmatrix}(\omega_{2^{n}}^0)^0&(\omega_{2^n}^0)^1&\cdots&(\omega_{2^n}^0)^{2^n-1}\\(\omega_{2^n}^1)^0&(\omega_{2^n}^1)^1&\cdots&(\omega_{2^n}^1)^{2^n-1}\\\vdots&\vdots&\ddots&\vdots\\(\omega_{2^n}^{2^n-1})^0&(\omega_{2^n}^{2^n-1})^1&\cdots&(\omega_{2^n}^{2^n-1})^{2^n-1}\end{bmatrix}\begin{bmatrix}a_0\\ a_1\\\cdots\\a_{2^n-1}\end{bmatrix}=\begin{bmatrix}f(\omega_{2^n}^0)\\f(\omega_{2^n}^1)\\\vdots\\f(\omega_{2^n}^{2^n-1})\end{bmatrix}$$
令左侧矩阵为 $\mathbf F$,系数列向量为 $\mathbf a$,点值表示为 $\mathbf d$,则整个式子可以简写为 $\mathbf {Fa=d}$,现在考虑在左右两边相乘某个矩阵 $\mathbf E$ 使得 $$\mathbf{EFa=Ed=a}$$。
我们知道矩阵 $\mathbf F(i,j)=(\omega_{2^n}^{i})^j$(计数从 $0$),同时我们构造 $\mathbf E(i,j)=(\omega_{2^n}^{-i})^j$,那么
$$\mathbf{EF}(i,j)=\sum_{k=0}^{2^n-1}\mathbf E(i,k)\mathbf F(k,j)=\sum_{k=0}^{2^n-1}\omega_{2^n}^{k(i-j)}$$
它的结果有两种可能:当 $i=j$ 时,我们就是把 $2^n$ 个 $1$ 加起来,结果是 $2^n$;否则的话,我们要研究的就是某个等比数列的求和,结论是这玩意儿等于 $0$。所以
$$\mathbf{EF}(i,j)=2^n[i=j]$$
或者说
$$\mathbf{EF}=2^n\mathbf I$$
其中 $\mathbf I$ 是单位向量。这意味着
$$2^{-n}\mathbf E=\mathbf F^{-1}$$
我们就发现了 FDFT 的逆向结论:把所有单位根取成其复共轭,把结果除以 $2^n$,即从点值表示法还原了系数。
实现细节
首先端上来递归的代码(来源)。
void fft(int n, complex<double>* buffer, int offset, int step, complex<double>* epsilon)
{
if(n == 1) return;
int m = n >> 1;
fft(m, buffer, offset, step << 1, epsilon);
fft(m, buffer, offset + step, step << 1, epsilon);
for(int k = 0; k != m; ++k)
{
int pos = 2 * step * k;
temp[k] = buffer[pos + offset] + epsilon[k * step] * buffer[pos + offset + step];
temp[k + m] = buffer[pos + offset] - epsilon[k * step] * buffer[pos + offset + step];
}
for(int i = 0; i != n; ++i)
buffer[i * step + offset] = temp[i];
}
但是我们都知道递归的常数肥肠大。这时我们需要一点优化。
首先我们可以先确定每个元素在递归结束时都被分到了那里:
000 001 010 011 100 101 110 111
000 010 100 110|001 011 101 111
000 100|010 110|001 101|011 111
000|100|010|110|001|101|011|111
我们发现这玩意儿就是反转了一下二进制位。这种旋转表可以在 $\mathrm O(n)$ 的复杂度内进行递推($l$ 与 $\log n$ 同阶):
for(int i=0;i<(1<<l);++i) ce[i]=((i&1)<<(l-1))|(ce[i>>1]>>1);
然后,我们可以端上热乎的 FDFT 代码:
void fdft(clex *f,int n,bool inv){
for(int i=0;i<n;++i) if(i<ce[i]) swap(f[i],f[ce[i]]);
for(int l=2;l<=n;l<<=1){
double te=2*pi/(double)l;
clex step=clex(cos(te),sin(te)*(inv?-1:1));
for(int i=0;i<n;i+=l){
clex cur=1;
for(int j=0;j<(l>>1);++j){
clex g=f[i+j],h=f[i+j+(l>>1)];
f[i+j]=g+cur*h,f[i+j+(l>>1)]=g-cur*h;
cur*=step;
}
}
}
if(inv) for(int i=0;i<n;++i) f[i]/=n;
}
在 $\mathrm O(n\log n)$ 复杂度内解决了问题,通过将递归的操作留在循环内来实现。
最终,通过在点值内进行乘法再倒回来的方式,我们即实现了 $\mathrm O(n\log n)$ 的多项式乘法。完整实现请读者自己实现,
一些加速
为了让我们不用过于慢速的 double
,我们需要加速。
在素数域 $\mathbb Z(p=r\cdot 2^k+1)$ 里存在原根 $g$,满足对于 $i=0,1,2,\cdots,p-2$,$g^i$ 两两不同,且 $g^{p-1}\equiv 1\pmod p$。令 $g_n\equiv g^{(p-1)/n}\pmod p$,则
$$\omega_{hn}^{hk}=\omega_{n}^{k}\Leftrightarrow g^{hk}_{hn}\equiv g^k_n\pmod p$$
$$\omega_{n}^{k+n/2}=-\omega_{n}^k\Leftrightarrow g_n^{k+n/2}\equiv -g_n^k\pmod p$$
$$\omega_{n}^{k+n}=\omega_{n}^k\Leftrightarrow g_n^{n+k}\equiv g^k_n\pmod p$$
$$\sum_{k=0}^{2^n-1}\omega_{2^n}^{k(i-j)}=2^n[i=j]\Leftrightarrow \sum_{k=0}^{2^n-1}g_{2^n}^{k(i-j)}\equiv 2^n[i=j]\pmod p$$
证明留给读者自己完成。所以我们可以把复数换成素数域,这里给出三个素数:
$$g=3,p=119\times 2^{23}+1=998244353$$
$$g=3,p=479\times 2^{21}+1=1004535809$$
$$g=3,p=5\times 2^{25}+1=167772161$$
于是乎,我们可以这样改写 FDFT,改为 FNTT(快速数论变换):
typedef long long ll;
const ll Q=119<<23|1,g=3;
void ew(ll &x,ll Q){
x%=Q,x+=Q,x%=Q;
}
ll pw(ll a,ll n,ll Q){
ll ans=1;
for(a%=Q;n;n>>=1,a=a*a%Q) if(n&1) ans=ans*a%Q;
return ans;
}
void fntt(ll *f,int n,bool inv){
for(int i=0;i<n;++i) if(i<rv[i]) swap(f[i],f[rv[i]]);
for(int k=2;k<=n;k<<=1){
ll step=pw(pw(g,(inv?Q-2:1),Q),(Q-1)/k,Q);
for(int i=0;i<n;i+=k){
ll cur=1;
for(int j=0;j<(k>>1);++j,cur*=step,ew(cur,Q)){
ll g=f[i+j],h=f[i+j+(k>>1)];
f[i+j]=g+cur*h,ew(f[i+j],Q);
f[i+j+(k>>1)]=g-cur*h,ew(f[i+j+(k>>1)],Q);
}
}
}
if(inv){
ll iv=pw(n,Q-2,Q);
for(int i=0;i<n;++i) f[i]*=iv,ew(f[i],Q);
}
}
任意模数
现在我们假设这个素数域 $p$ 不能表示为 $r\cdot 2^k+1$(而且 $g$ 极难求)。所以,这里给出一种方案:即通过多种模数进行计算,最后通过中国剩余定理进行合并,而后模上该数。
由于这里的数很大,所以附上一个 __int128
的读写模板。
typedef __int128 _i;
const _i Q[3]={119<<23|1,479<<21|1,5<<25|1},g=3;
_i get(){
_i x=0,f=1;
char ch=getchar();
while(!isalnum(ch)) f=(ch=='-'?-1:1),ch=getchar();
while(isalnum(ch)) x=x*10+ch-'0',ch=getchar();
return x*f;
}
void put(_i x){
if(x<0) putchar('-'),x=-x;
if(x>=10) put(x/10);
putchar(x%10+'0');
}
void ew(_i &x,_i Q){
x%=Q,x+=Q,x%=Q;
}
_i pw(_i a,_i n,_i p){
_i ans=1;
for(a%=p;n;n>>=1,a=a*a%p) if(n&1) ans=ans*a%p;
return ans;
}
void fntt(_i *f,int n,bool inv,_i Q){
for(int i=0;i<n;++i) if(i<rv[i]) swap(f[i],f[rv[i]]);
for(int k=2;k<=n;k<<=1){
_i step=pw(pw(g,(inv?Q-2:1),Q),(Q-1)/k,Q);
for(int i=0;i<n;i+=k){
_i cur=1;
for(int j=0;j<(k>>1);++j,cur*=step,ew(cur,Q)){
_i g=f[i+j],h=f[i+j+(k>>1)];
f[i+j]=g+cur*h,ew(f[i+j],Q);
f[i+j+(k>>1)]=g-cur*h,ew(f[i+j+(k>>1)],Q);
}
}
}
if(inv){
_i iv=pw(n,Q-2,Q);
for(int i=0;i<n;++i) f[i]*=iv,ew(f[i],Q);
}
}
void exgcd(_i a,_i b,_i &x,_i &y){
if(b==0) x=1,y=0;
else exgcd(b,a%b,y,x),y-=a/b*x;
}
void p_prod(){
for(int i=0;i<=h;++i){
for(int j=1;j<=3;++j) af[j%3][i]=af[0][i]%Q[j%3];
for(int j=1;j<=3;++j) bf[j%3][i]=bf[0][i]%Q[j%3];
}
for(int i=0;i<h;++i) rv[i]=(rv[i>>1]>>1)|((i&1)*(h>>1));
for(int i=0;i<3;++i){
fntt(af[i],h,0,Q[i]),fntt(bf[i],h,0,Q[i]);
for(int j=0;j<h;++j) af[i][j]*=bf[i][j],ew(af[i][j],Q[i]);
fntt(af[i],h,1,Q[i]);
}
for(int i=0;i<h;++i){
_i qQ=Q[0]*Q[1]*Q[2];
for(int j=0;j<3;++j){
_i q=qQ/Q[j],tmp,ivq;
exgcd(q,Q[j],ivq,tmp);
f[i]+=af[j][i]*q*ivq,ew(f[i],qQ);
}
}
}
多项式牛顿迭代
可以先参考实数域上的牛顿迭代。现在我们规定符号 $\circ$ 代表函数符合,例如 $G\circ F$ 就等同于 $G(F(x))$。
设多项式 $G$ 以及需要求出的 $F$,以及需要进行取模的 $n$,我们想得到
$$G\circ F\equiv 0\pmod{x^n}$$
现在我们假设
$$F_0\equiv F\pmod{x^{\lceil n/2\rceil}}$$
和
$$F_1\equiv F\pmod{x^n}$$
则由泰勒公式,
$$G\circ F=\sum_{k\geq 0}\frac{G^{(k)}\circ F_0}{k!}(F-F_0)^k$$
因为
$$F-F_0\equiv 0\pmod{x^{\lceil n/2\rceil}}$$
所以当 $k\geq 2$ 时(两个最高项乘起来大于等于 $n$)
$$(F-F_0)^k\equiv 0\pmod{x^n}$$
故
$$G\circ F_1\equiv G\circ F_0+(G'\circ F_0)(F_1-F_0)\pmod{x^n}$$
由于
$$G\circ F_1\equiv 0\pmod {x^n}$$
则
$$F_1\equiv F_0-\frac{G\circ F_0}{G'\circ F_0}$$
这就是适用于多项式的牛顿迭代公式。
多项式求逆
设 $FG\equiv 1\pmod{x^n}$,变换一下就得到 $F-\frac{1}{G}\equiv 0\pmod {x^n}$,且令 $H(t)=F-1/t$。代入并化简就得到
$$G_1\equiv G_0(2-G_0F)\pmod{x^n}$$
注意,零次多项式的逆就是 $p$ 对零次项的逆元。别搞错了。
void inve(ll *a,ll *b,int n){
if(n==1) return b[0]=pw(a[0],Q-2),void();
int h=(n+1)>>1,L=1;
static ll c[N];
inve(a,b,h);
for(;L<(n<<1);L<<=1);
for(int i=0;i<L;++i) rv[i]=(rv[i>>1]>>1)|((i&1)*(L>>1));
fill(b+h,b+L,0),fntt(b,L,0);
copy(a,a+n,c),fill(c+n,c+L,0),fntt(c,L,0);
for(int i=0;i<L;++i) b[i]=ew(b[i]*ew(2-b[i]*c[i]));
fntt(b,L,1),fill(b+n,b+L,0);
}
多项式自然对数
设 $B\equiv \ln A\pmod{x^n}$,现在在 $\mathbb Z(x^n)$ 内 $\ln$ 也变成了一个可以应用 $(fg)'=f'g+g'f,f\circ g=(f'\circ g)g'$ 的函数,所以在式子两边取对数得到
$$B'\equiv \ln'(A)A'\equiv A'A^{-1}\pmod{x^n}$$
最后在式子两端积一下分就得到
$$B\equiv\int A'A^{-1}\mathrm dx\pmod{x^n}$$
这里附上了求导和积分和乘法的代码。都没有后效性。
void mult(ll *a,ll *b,ll *c,int n,int m){
int L;
for(L=1;L<=n+m;L<<=1);
for(int i=0;i<L;++i) rv[i]=(rv[i>>1]>>1)|((i&1)*(L>>1));
fntt(a,L,0),fntt(b,L,0);
for(int i=0;i<L;++i) c[i]=ew(a[i]*b[i]);
fntt(c,L,1),fntt(a,L,1),fntt(b,L,1);
}
void deriv(ll *f,ll *Df,int n){
for(int i=1;i<n;++i) Df[i-1]=ew(f[i]*i);
Df[n-1]=0;
}
void integ(ll *f,ll *Sfdx,int n){
Sfdx[0]=0;
for(int i=1;i<n;++i) Sfdx[i]=ew(f[i-1]*pw(i,Q-2));
}
void logar(ll *f,ll *g,int n){
static ll Df[N],invf[N],Dg[N];
deriv(f,Df,n),inve(f,invf,n),mult(Df,invf,Dg,n,n),integ(Dg,g,n);
}
多项式指数
现在我们试图求出 $B\equiv e^A\pmod{x^n}$ 了。套用逆元的技巧,我们设 $H(t)=\ln t-A$,则我们发现我们要求 $H(t)$ 的零点,那么代入多项式牛顿迭代公式,就有
$$F_1\equiv F_0-\frac{H\circ F_0}{H'\circ F_0}=F_0-F_0(\ln F_0-A)=F_0(\ln F_0-A+1)$$
把它塞进牛顿迭代的模板里面就彳亍。注意这些代码都是分着题写的,模板不是CC-CV的,所以不能直接拿来用。最完整的模板应该会在末尾。
void expon(ll *f,ll *g,int n){
if(n==1){
assert(f[0]==0);
g[0]=1;
return;
}
static ll h[N],k[N];
int m=(n+1)>>1,L;
expon(f,g,m),logar(g,k,n);
for(int i=0;i<n;++i) h[i]=ew(f[i]-k[i]);
h[0]=ew(h[0]+1),pre(L,n<<1);
fill(g+m,g+L,0),fntt(g,L,0);
fill(h+n,h+L,0),fntt(h,L,0);
for(int i=0;i<L;++i) g[i]=ew(g[i]*h[i]);
fntt(g,L,1),fill(g+n,g+L,0);
}
线性代数
一般来说,在线性代数里的内容,总会和线性代数里的东西沾点边——例如矩阵消元,线性基,等等。至于实际的数学证明,让它们见鬼去吧。
高斯消元
我们首先将我们的问题规范化:已知方阵 $\mathbf A$ 和列向量 $\mathbf v$,我们的任务是寻找列向量 $\mathbf x$ 满足 $\mathbf A\mathbf x=\mathbf v$。一般来说,$\mathbf A$ 和 $\mathbf v$ 通常被压在同一个矩阵中,即 $\mathbf M=(\mathbf A|\mathbf v)$,以便于我们的操作。最后我们一定会得到一个矩阵 $\mathbf M'=(\mathbf I|\mathbf x)$。
按照标准的想法,我们要把 $\mathbf A$ 削成一个单位矩阵 $\mathbf I$。我们按照各行来考虑这个问题。首先,对于第 $i$ 行,我们想把第 $i$ 行以外的,或者第 $i$ 行以上的第 $i$ 列系数削成 $0$。(前一种算法是高斯-约尔当消元法,后一种则是高斯消元法。)于是,我们有一个操作:把一行乘以某个系数,加到另一行上;这就是所谓的加减消元法了。我们依次对每行消元,最终得到我们想要的东西。
很不幸的是,这个问题有可能无解,也有可能多解。 在消元过程中,我们常常会把第 $i$ 列的所有系数都变成 $0$,不管它,消完元后,我们只需要判断一件事儿:对于某条方程,如果系数等于 $0$ ,如果等式的右侧等于 $0$,那么多解;否则无解。
以下是代码。
// 洛谷 p2455 - [SDOI2006] 线性方程组
// Alea
#include <bits/stdc++.h>
#define fs(x) fixed<<setprecision(x)
using namespace std;
const double eps=1e-3;
int n;
double M[107][107];
bool vis[107];
int main(){
cin>>n;
for(int i=1;i<=n;++i) for(int j=1;j<=n+1;++j) cin>>M[i][j],vis[i]=0;
bool ok=1,inc=0;
// Gauss-Jordan Elimination.
for(int i=1;i<=n;++i){
int ii=i;
for(int j=1;j<=n;++j){
if(!vis[j]&&fabs(M[ii][i])<fabs(M[j][i])) ii=j;
}
for(int j=1;j<=n+1;++j) swap(M[ii][j],M[i][j]);
if(fabs(M[i][i])<eps){
ok=0;
continue;
}
vis[i]=1;
// We've chosen the main unknown. In the normal case.
double mii=M[i][i];
for(int j=1;j<=n+1;++j) M[i][j]/=mii;
for(int j=1;j<=n;++j){
if(j==i) continue;
double c=M[j][i];
for(int k=1;k<=n+1;++k) M[j][k]-=M[i][k]*c;
}
}
if(ok){
for(int i=1;i<=n;++i) cout<<"x"<<i<<"="<<fs(2)<<M[i][n+1]<<endl;
}else{
for(int i=1;i<=n;++i) if(fabs(M[i][i])<eps&&fabs(M[i][n+1])>eps)inc=1;
cout<<(inc?"-1":"0")<<endl;
}
return 0;
}
高斯消元法同时适用于多种数的空间,甚至异或空间!(但是异或空间就是 $\mathbb Z_2$。)这给了我们极大的扩展性来完成一些线性方程组的工作,包括求解神奇的马尔科夫链。
线性基
我们首先来介绍一下我们的任务:我们将每个正整数视作一个向量,分量为它各个二进制上的值,同时仅允许在线性叠加的时候使用 $0$ 和 $1$ 做系数,异或代替加法,这样我们就得到了一个向量空间 $V(\mathbb Z_2)$。
首先的问题是——如何通过给定的一组数,计算线性基。我们采用贪心算法,因为前一节介绍的高斯消元法会亿点点麻烦。
对于第 $i$ 个基向量 $v_i$,我们希望它最高位上的数为 $1$,于是乎我们可以从大到小便利各个存在的基向量,如果要插入的数 $x$ 在第 $i$ 位上同样是 $1$,那么就消掉它。最终我们的算法一定会在 $x$ 的某个为 $1$ 的位上停下来,或者不停下来——前者可以直接作为一个基向量来使用,而后者就表明我们已经满员了。
可以证明,对于二进制数长度为 $n$ 的,它的基向量最多不会超过 $n$ 个。证明请参考一般的线性代数教科书。
最后,我们补充一下插入的代码。注意那些看上去奇怪的地方。
void inse(ll x,int y){
for(ll i=63;i+1;--i){
if(!(x>>i)) continue;
if(p[i]) x^=p[i];
else return p[i]=x,void();
}
}
你可以将各种的存在性定则或者需要对二取余的存在性定则的最优化补齐视作线性基(虽然本人并不知道这是什么鬼意思。)例如,你甚至可以把完全平方数的一些可能性内容塞进线性基,前提是要保证亲爱的素数的大小不要太大——个数也不要太多。
至于最大值,我们可以挨个儿的异或,然后通过比较更大值找到线性基。
ll maxe(){
ll ans=0;
for(ll i=63;~i;--i) ans=max(ans,ans^p[i]);
return ans;
}
但是,对于第 $k$ 大值,我们需要首先对线性基处理一下:首先,我们需要将线性基的各个数削成一条近似的线:
void init(){
for(ll i=1;i<=63;++i) for(ll j=0;j<i;++j)
if(p[i]&(1<<j)) p[i]^=p[j];
for(ll i=0;i<64;++i) S+=(p[i]!=0);
}
这样,我们就可以使用“二进制对应法”求出第 $k$ 大的值了。这里我们需要假定 $k$ 从 $0$ 开始计。与此同时我们要确保异或的东西确乎是一个数值。
ll kth(ll k){
if(k>=(1<<S)) return -1; // we cannot calculate it.
ll ans=0;
for(ll i=0;i<64;++i) if(p[i]!=0) ans^=(p[i]*(k&1)),k>>=1;
return ans;
}
对于存在性和求并问题,我们直接往里插入进去就好。
博弈论
在这里,我们研究的是组合博弈中的公平组合博弈。它们通常非常抽象,因为状态很多。
公平博弈(impartial game)指满足如下条件的组合博弈:
- 在任意确定状态下,所有参与者可选择的行动完全相同,仅取决于当前状态,与身份无关;
- 博弈中的同一个状态不可能多次抵达,博弈以参与者无法行动为结束,且博弈一定会在有限步后以非平局结束。
公平博弈总是对称博弈。
——OI Wiki
博弈图
通常来说,公平组合游戏都涉及到一些状态;我们可以将这些状态放在一张图上,称为状态图。
我们定义必胜客必胜点和必败点(相对于这些节点自己)。我们可以轻松地发现:
- 出度为 $0$ 的点为必败点。因为在这一步中,该操作者没法进行转移。
- 对于一个点,如果它的所有后继点都是必胜点,那么它是一个必败点;否则,它是一个必胜点。因为,点在这一步进行的转移可以取对于自己最优的转移(也就是对方的必败点);否则就死掉了。
下面我们直接定义 SG 函数,这里 $u$ 和 $v$ 是状态图上的点,而 $G$ 指的就是状态图,而 $\mathrm{mex} S$ 代表从集合 $S$ 中取不属于该集合的最小自然数。
$$\mathrm{SG}(u)=\mathop{\mathrm{mex}}\limits_{(u,v)\in G}\{\mathrm{SG}(v)\}$$
于是乎, $\mathrm{SG}$ 就用一种较为另类的方式描述了必胜或必败点。因为
- 对于出度为 $0$ 的点 $u$,$S$ 是一个空集,所以 $\mathrm{SG}(u)=0$。
- 对于出度不为 $0$ 的点 $u$,倘若它直接相连的点 $v$ 中存在 $\mathrm{SG}(v)=0$,则其 $\mathrm{SG}(u)\neq 0$;否则 $\mathrm{SG}(u)=0$。
容易发现,$\mathrm{SG}$ 函数值是否为 $0$ 与该点是否必败相一致。
可是它有什么用呢?
Bash 和 Nim
首先,我们来研究 Bash 游戏,即对于 $n$ 个石头,每次取最多 $m$ 个石头,判断是否必胜。
我们当然可以先画出一个图,但是画图太抽象了,所以我们可以用 dp 进行解决。
首先,对于 $1\leq n\leq m$,一次性就拿完了,所以 $\mathrm{SG}(n)\neq 0$。分析以后,你会发现 $\mathrm{SG}(n)=n$。
而后,对于 $n=m+1$,我们发现无论拿多少个都会剩下让乙一次性拿完的方案,所以 $\mathrm{SG}(m+1)=0$。
手推一下,我们发现 $\mathrm{SG}(n)=n\bmod(m+1)$。
现在我们研究 Nim 游戏,它包含 $n$ 堆石头 $a_1\cdots a_n$,每次可以不限量的取石头,问先手能否必胜。
如果 $n=1$,那么很简单:一次性取完,我们就赢了,即 $\mathrm{SG}(a_k)=a_k$。然而,不幸,$n$ 会很大,而完整的画一个图又不怎么可能。怎么办呢?
我们介绍这样一个定理。
一个游戏的 SG 函数值等于各个不相交游戏的 SG 函数值的异或和。
用数学来写是这样:
$$\mathrm{SG}\left(\bigcup_{k}X_k\right)=\bigoplus_k{\mathrm{SG}(X_k)}$$
所以,我们直接把这些 $SG$ 异或起来,看看能不能得到 $0$,如果 $0$,就死了,因为按照上面的定义,$\mathrm{SG}$ 与是否能必胜相一致。
杂项
主定理
对于 $T(n)=aT(n/b)+O(n^d)$,有如下结论:
$$T(n)=\begin{cases}O(n^d)&b^d>a\\ O(n^d\log n)&b^d=a\\ O(n^{\log_b a})&b^d<a\end{cases}$$
证明。我们令 $k=\lceil\log_b{n}\rceil$,此时 $O$ 会忽略掉向上取整带来的常数影响,则
- 若 $b^d=a$,则 $T(b^k)=aT(b^{k-1})+O(b^{kd})$,等式两端同除以 $a^k$ 得到 $T(b^k)/a^k=T(b^{k-1})/(a^{k-1})+O(1)$,显而易见 $T(b^k)/a^k=O(k)$,即 $T(n)=T(b^k)=O(ka^k)=O(nk)=O(n\log n)$,抽出最左和最右端即得。
- 若 $b^d>a$,用归纳法,对于 $k=0$,显然成立;假定 $T(b^k)=O(b^{kd})$ 在 $k\leq K$ 时成立,对于 $K+1$,代入得到 $T(b^{K+1})=O(ab^{Kd}+b^{(K+1)d})$,由于 $a<b^d$,那么显然 $ab^{Kd}<b^{(K+1)d}$,所以 $O$ 里面的第二项主导整个 $O$,即 $T(b^{K+1})=O(b^{(K+1)d})$,所以按照归纳法,$T(n)=O(n^d)$ 成立。
- 对于 $b^d<a$,仍然用归纳法(略去基底),$T(b^{K+1})=O(a^{K+1}+b^{(k+1)d})$ 的第一项主导整个 $O$,即 $T(b^{K+1})=O(a^{K+1})$。更改一下形式,就是 $T(n)=O(n^{\log_b{a}})$。
初赛场上可以现推。