Tiw_Air_OAO 数学杂题选讲

伯努利数求自然数幂前缀和

$$n^k=[x^k]\frac{1}{1-nx}=k![x^k]e^{nx}$$

$\text{EGF}$（指数型生成函数）形式更为常用，不如说 $\text{OGF}$（普通型生成函数）几乎没用。

例：

$$\sum_{i=0}^{n}i^k=k![x^k]\sum_{i=0}^n e^{ix}=k![x^k]\frac{e^{(n+1)x}-1}{e^{x}-1}$$

后面的分式可以拆成 $\dfrac{x}{e^{x}-1}$（只有常数项不为零的形式幂级数才可以求逆）与 $\dfrac{e^{(n+1)x}-1}{x}$ 的乘积。

前者与 $n$ 无关，可以预先求出。事实上，这就是伯努利数的 $\text{EGF}$ 。

与斯特林数的异同

斯特林数也可以用于求自然数幂和。

$$\sum_{i=0}^n i^k=\sum_{j=0}^k{k\brace j}{n+1\choose j+1}j!$$

但其实它们是差不多的！回忆斯特林数的 $\text{EGF}$ ：

$$\sum_{i=0}^{+\infty}{n\brace k}\frac{x^n}{n!}=\frac{(e^x-1)^k}{k!}$$

那么：

$$\begin{aligned} \sum_{i=0}^{n}i^k&=k![x^k]\sum_{i=0}^n e^{ix}\\ &=k![x^k]\sum_{i=0}^n (e^{x}-1+1)^{i}\\ &=k![x^k]\sum_{i=0}^n \sum_{j=0}^{i}{i\choose j}(e^x-1)^{j}\\ &=k![x^k]\sum_{j=0}^n {n+1\choose j+1}(e^x-1)^{j}\\ &=k![x^k]\sum_{j=0}^n \frac{(e^x-1)^{j}}{j!}{n+1\choose j+1}j!\\ &=k![x^k]\sum_{j=0}^k \frac{(e^x-1)^{j}}{j!}{n+1\choose j+1}j!\\ &=\sum_{j=0}^k{k\brace j}{n+1\choose j+1}j!\\ \end{aligned}$$

这其中关键的一步就是 $\sum_{j=0}^n\limits \frac{(e^x-1)^{j}}{j!}{n+1\choose j+1}j!\to \sum_{j=0}^k\limits \frac{(e^x-1)^{j}}{j!}{n+1\choose j+1}j!$ 的指标代换，这事实上也是斯特林反演的关键。它成立的原因是：$e^x-1$ 常数项是零，从而当 $j>k$ 时 $[x^k](e^x-1)^{j}=0$。

两个做法在第二步出现分歧，伯努利数直接对原式化简，而斯特林反演则是将关于 $e^x$ 的函数代换成 $e^x-1$。

第二个做法可以避免求逆元。

[NOI Online 2021 提高组] 愤怒的小 N

设集合 $S$ 对应 $F_S(x)=\sum_{u\in S}e^{ux}$ ，那么 $k![x^k]F_S(x)$ 就是 $S$ 中所有数的 $k$ 次幂之和。

我们想要求出长度为 $n$ 的前缀对应的集合 $S$。由于字符串的特殊性，可以考虑倍增。

分别维护长度为 $2^t$ 的前缀 $0/1$ 对应的集合的生成函数 $P_t(x),Q_t(x)$，有如下转移式成立：

$$P_{t+1}(x)=P_t(x)+e^{2^tx}Q_t(x),P_0(x)=1 \\ Q_{t+1}(x)=Q_t(x)+e^{2^tx}P_t(x),Q_0(x)=0$$

仔细瞪一下，作换元 $R_t=P_t+Q_t$ 与 $T_t=P_t-Q_t$ 能够给出更简洁的递推式：

$$R_{t+1}=(1+e^{2^tx})R_t,R_0=1 \\ T_{t+1}=(1-e^{2^tx})T_t,T_0=1$$

发现 $R_t=\sum_{i=0}^{2^k} e^x$ ，就是自然数幂和，同时 $(1-e^{2^tx})$ 常数项为 $0$ ，因此 $Q_t(x)≡0\pmod {x^k} \ (t≥k)$，我们只需要求前 $k$ 项。拿 $Q_t$ 凑 $n=(n_{len−1}⋯n_0)_2$ 的 $\text{bit}$，得到：

$$Q(x)=\sum_{i=0}^{l−1}[n_i=1](−1)^{count((n−(n_i⋯n_0)_2)}e^{(n−(n_i⋯n_0)_2)x}Q_i(x).$$

事实上指标 $i$ 的有效范围只到 $k$，直接 $O(k^3)$ 算出来。另一方面 $P(x)$ 的自然数幂和可以拉格朗日插值，最终复杂度是 $O(\log n+k^3)$ 。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n,k;
char s[500005];
const int md=1e9+7;
inline int pwr(int x,int y){
	int res=1;
	while(y){
		if(y&1)res=1ll*res*x%md;
		x=1ll*x*x%md;y>>=1;
	}
	return res;
}
int a[505],fac[505],ifac[505],inv[505],Pwr[500005];
int T[505][505],e[505],t[505],tmp[505];
inline void mul(int *a,int *b,int *c){
	for(int i=0;i<k;i++)c[i]=0;
	for(int i=0;i<k;i++){
		for(int j=0;i+j<k;j++)c[i+j]=(c[i+j]+1ll*a[i]*b[j])%md;
	}
}
int pre[505],suf[505];
inline int lagrange(int *y,int x){
	int res=0;pre[0]=x;suf[k+1]=1;
	for(int i=1;i<=k;i++)pre[i]=1ll*pre[i-1]*(x-i)%md;
	for(int i=k;i>=1;i--)suf[i]=1ll*suf[i+1]*(x-i)%md;
	for(int i=0;i<=k;i++){
		int tmp=1ll*ifac[i]*ifac[k-i]%md*((k-i)&1?md-1:1)%md;
		if(i)tmp=1ll*tmp*pre[i-1]%md;
		if(i!=k)tmp=1ll*tmp*suf[i+1]%md;
		res=(res+1ll*y[i]*tmp)%md;
	}
	return res;
}
inline void init(){
	Pwr[0]=fac[0]=fac[1]=inv[0]=inv[1]=ifac[0]=ifac[1]=1;
	for(int i=2;i<=k;i++)fac[i]=1ll*fac[i-1]*i%md;
	for(int i=2;i<=k;i++)inv[i]=1ll*(md-md/i)*inv[md%i]%md;
	for(int i=2;i<=k;i++)ifac[i]=1ll*ifac[i-1]*inv[i]%md;
	for(int i=1;i<=n;i++)Pwr[i]=2ll*Pwr[i-1]%md;
	T[0][0]=1;
	for(int i=1;i<k;i++){
		for(int j=1;j<k;j++)e[j]=1ll*ifac[j]*pwr(Pwr[i-1],j)%md;
		for(int j=1;j<k;j++)e[j]=(md-e[j])%md;mul(T[i-1],e,T[i]);
//		cout<<i<<endl;
//		for(int j=0;j<k;j++)cout<<T[i][j]<<" ";cout<<endl;
	}
	int sum=0,op=-1;
	for(int i=n-1;i>=k;i--){
		if(s[i])sum=(sum+Pwr[i])%md,op*=-1;
	}
	for(int i=k-1;~i;i--){
		if(s[i]){
			for(int j=0;j<k;j++)e[j]=1ll*ifac[j]*pwr(sum,j)%md;
			mul(T[i],e,tmp);
			for(int j=0;j<k;j++)t[j]=(t[j]+1ll*op*tmp[j]+md)%md;
			sum=(sum+Pwr[i])%md;op*=-1;
		}
	}
	for(int i=0;i<k;i++)t[i]=1ll*t[i]*fac[i]%md;
//	for(int i=0;i<k;i++)cout<<t[i]<<" ";cout<<endl;
	for(int i=0;i<k;i++)scanf("%d",&a[i]);
	for(int i=0;i<k;i++){
		tmp[0]=(i==0);
		for(int j=1;j<=k;j++)tmp[j]=(tmp[j-1]+pwr(j,i))%md;
		t[i]=1ll*(t[i]+lagrange(tmp,sum-1))%md*pwr(2,md-2)%md;
	}
//	for(int i=0;i<k;i++)cout<<t[i]<<" ";cout<<endl;
	int res=0;
	for(int i=0;i<k;i++)res=(res+1ll*a[i]*t[i])%md;
	printf("%d\n",res);
}
int main(){
	scanf("%s",s);n=strlen(s);
	reverse(s,s+n);for(int i=0;i<n;i++)s[i]-='0';
	scanf("%d",&k);init();


	return 0;
}

hdu6593 Coefficient

可以换元 $g(x)=\dfrac{1}{b}f(x{+}x_0)=\dfrac{1}{c+e^{ax}}$，因此我们就假定 $b=1,\;d=x_0=0$ 吧。

根据 $f'(x)=\dfrac{-ae^{ax}}{(c+e^{ax})^2}$，可以想到接下来的形式必然是 $\sum\dfrac{f_{j}e^{jax}}{(c+e^{ax})^{j+1}}$ 。系数的变化规律是：

$$\text{derivative's }f_j=ja\cdot(f_j-f_{j-1})$$

$a$ 是统一的，可忽略之，最后答案乘 $a^n$ 。初值 $f_0=1$ 。最后的答案就是 $\dfrac{1}{n!}\sum_{j\geqslant 0}\dfrac{f_j}{(c{+}1)^{j+1}}$ 。

考虑系数的 “移动”（因为是线性变换），可见路径权值是 $(-1)^k\prod_{i=1}^{k}i^{q_i}$，其中 $q_i\ne 0$ 且 $\sum q_i=n$，贡献到 $f_k$ 。这就可以用 $\rm OGF$ 轻易描述了：

$$f_k=[x^n](-1)^k\prod_{j=1}^{k}\frac{jx}{1-jx}$$

我不熟悉它，所以我没能指出：$\prod_{j=1}^{k}\dfrac{x}{1-jx}$ 是 第二类斯特林数 在第 $k$ 列上的 $\rm OGF$ 。其实想想 $\prod_{i=1}^{k}i^{q_i}$ 不就是每次要么放进新盒子、要么放进已有的盒子，第二类斯特林数吗！

“你看到的只是表象。” —— 太阳神 $\textsf{Tiw-Air-OAO}$

因此 $f_k=k!(-1)^k{n\brace k}$ 即得。用斯特林数在第 $n$ 行上的 $\rm OGF$ 去算即可。

再或者，我们可以换一种刺杀的方式。直接从生成函数的角度考虑它，答案就是

$$[x^n]\frac{1}{c+e^{ax}}$$

吸取 “$\text{Binomial Sum}$” 的经验，只需把 $F(x)=\dfrac{1}{1+c+x}$ 展开成

$$F(x)=(c{+}1)^{-1}\sum_{i=0}^{n}(-1)^i(c{+}1)^{-i}$$

我们先算出

$$[x^n](\text{e}^{ax}-1)^t=\sum_{i=0}^{t}{t\choose i}(-1)^{t-i}(ai)^n$$

Remark：把 $a$ 去掉后，其实 $\dfrac{1}{t!}(\text{e}^x{-}1)^t$ 就是第二类斯特林数在第 $t$ 列上的 $\rm EGF$ 。因此我们得到的也就是第二类斯特林数在第 $n$ 行上的 $\rm OGF$ 。

直接得到

$$ans=a^n\left[(c{+1})^{-1}\sum_{t=0}^{n}(c{+}1)^{-t}\sum_{i=0}^{t}{t\choose i}(-1)^{i}i^n\right]$$

系数 $\sum_{i=0}^{t}{t\choose i}(-1)^{i}i^n$ 是简单卷积，然后对 $q$ 个 $(c{+}1)^{-1}$ 多点求值。复杂度 $\mathcal O(n\log n+q\log^2 n)$ 。

点击查看代码

#4126. 国王奇遇记加强版之再加强版

求：

$$\left(\sum_{i=1}^ni^m\ \times m^i\right)\bmod (10^9+7)\ \ (n \leq 10^9,m\leq500000)$$

用 $m![x^m]e^{ix}=i^m$ 把 $m$ 放到指数外，有：

$$\begin{aligned} \textit{ans} &= \sum_{i=1}^ni^m\times m^i\\ &= m![x^m]\sum_{i=1}^ne^{ix}m^i\\ &= m![x^m]\frac{1-m^{n+1}e^{(n+1)x}}{1-me^x}. \end{aligned}$$

其中分母是常多项式，可以分析到答案是关于 $n+1$ 的多项式。令 $F(x)=\dfrac{1}{1−me^x}$，则：

$$\textit{ans}=m!(f_m-m^{n+1}\sum_{i=0}^m\frac{(n+1)^i}{i!}f_{m-i}).$$

可见，若令：

$$P(x)=\sum_{i=0}^m\frac{x^if_{m-i}}{i!}$$

则：

$$\textit{ans}=m!(P(0)-m^{n+1}P(n+1))$$

但是我们对 $P$，甚至对 $f$ 一无所知。怎么办？其实 $n$ 较小的时候，我们可以用答案最初的式子算出来，继而得到若干 $P(t)$ 与 $P(0)$ 的值的线性关系。再者，$P(x)$ 是 $m$ 次多项式，其 $m+1$ 阶差分为 $0$。所以：

$$\Delta^{m+1}P(0)=\sum_{i=0}^{m+1}(-1)^{m+1-i}\binom{m+1}{i}P(i).$$

利用已知线性关系解出 $P(0)$，得到点值，$\text{Lagrange}$ 差值。线性筛自然数幂、线性求一堆逆元，精细实现可以做到 $O(m)$ 。

#6718. 九个太阳「弱」化版

给定 $n,k$ ，满足 $k$ 是 $2$ 的幂，求

$$\sum_{k|i,0\le i\le n}{n\choose i}$$

对 $998244853$ 取模。

发现：

$$\sum_{0\le i\le n}{n\choose i}=(1+x)^n$$

进而 $\sum_{0\le i\le n}[k|i]{n\choose i}$ 就是 $(1+x)^n$ 中 $k$ 的倍数次项，即 $(1+x)^n\bmod (x^k-1)$ 的常数项，对 $(1+x)^n$ 做长度为 $k$ 的循环卷积即可。

「Karatsuba 算法」

设要计算 $2n−1$ 次多项式 $F(x)$ 和 $G(x)$ 的卷积，令 $F(x)=F_0(x)+x^nF_1(x)$，$G(x)=G_0(x)+x^nG_1(x)$，则：

$$F(x)G(x)=(F_0(x)+x^nF1(x))(G_0(x)+x^nG_1(x)) \\ =F_0(x)G_0(x)+x^n((F_0(x)+F_1(x))(G_0(x)+G_1(x))−F_0(x)G_0(x)−F_1(x)G_1(x))+x^{2n}F_1(x)G_1(x).$$

需要三次卷积，每次卷积递归计算。复杂度 $T(n)=3T(n/2)+O(n)=O(n^{\log_2 3})$。

不过还是过不了这题，还是要写 $\text{MTT}$ 或 $\text{long double}$ 的 $\text{FFT}$ 并每次取模。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const int md=998244853,wei=15;
int R[1000005];
int pw(int x,int y,int md){
	int ans=1;
	while(y){
		if(y&1)
			ans=1ll*ans*x%md;
		x=1ll*x*x%md;
		y>>=1;
	}
	return ans;
}
const long double pi=acos(-1.0);
struct CP{
	long double a,b;
	CP(){}
	CP(long double c,long double d){
		a=c;b=d;
	}
	friend CP operator + (const CP &x,const CP &y){
		return CP(x.a+y.a,x.b+y.b);
	}
	friend CP operator - (const CP &x,const CP &y){
		return CP(x.a-y.a,x.b-y.b);
	}
	friend CP operator * (const CP &x,const CP &y){
		return CP(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);
	}
	friend CP operator ~ (const CP &x){
		return CP(x.a,-x.b);
	}
};
void fft(CP a[],int len,int type){
	for(int i=0;i<len;i++)
		if(i<R[i])swap(a[i],a[R[i]]);
	for(int mid=1;mid<len;mid*=2){
		CP wn=CP(cos(pi/mid),type*sin(pi/mid));
		for(int j=0;j<len;j+=mid*2){
			CP w=CP(1,0);
			for(int k=0;k<mid;k++,w=w*wn){
				CP x=a[j+k],y=w*a[j+k+mid];
				a[j+k]=x+y;a[j+k+mid]=x-y;
			}
		}
	}
	if(type==1)return;
	for(int i=0;i<len;i++){
		a[i].a/=len;a[i].b/=len;
	}
}
CP A[1000005],B[1000005],C[1000005],D[1000005];
inline void mul(int a[],int b[],int n){
	int len=n;
	for(int i=0;i<len;i++){
		R[i]=(R[i>>1]>>1)|((i&1)*(len/2));
		A[i]=CP(a[i]>>wei,a[i]&((1<<wei)-1));
		B[i]=CP(b[i]>>wei,b[i]&((1<<wei)-1));
	}
	fft(A,len,1);fft(B,len,1);
	for(int i=0;i<len;i++){
		int j=((len-i)&(len-1));
		CP a0=(~A[j]+A[i])*CP(0.5,0),a1=(~A[j]-A[i])*CP(0,0.5),b0=(~B[j]+B[i])*CP(0.5,0),b1=(~B[j]-B[i])*CP(0,0.5);
		C[i]=a0*b0+a1*b1*CP(0.0,1.0);D[i]=a0*b1+a1*b0;
	}
	fft(C,len,-1);fft(D,len,-1);
	for(int i=0;i<n;i++){
		long long tmp1=(long long)(C[i].a+0.5)%md,tmp2=(long long)(D[i].a+0.5)%md,tmp3=(long long)(C[i].b+0.5)%md;
		a[i]=((tmp1<<(wei*2))+(tmp2<<wei)+tmp3)%md;
	}
}
int ans[1000005];
void solve(long long n,int k){
	if(n==1){
		ans[0]++;ans[1%k]++;
		return;
	}
	solve(n/2,k);
	mul(ans,ans,k);
	if(n&1){
		int tmp=ans[k-1];
		for(int i=k-1;i>=1;i--)ans[i]=(ans[i]+ans[i-1])%md;
		ans[0]=(ans[0]+tmp)%md;
	}
}
int main(){
	long long n;
	int k;
	scanf("%lld%d",&n,&k);
	solve(n,k);
	printf("%d",ans[0]);
	return 0;
}

hdu6355 fireflies

给定 $p_{1...n}$，以及 $\prod P_i$ 个点 $(x_1,x_2,...,x_n)$ ，$0\le x_i< p_i$ 。

一只萤火虫可以从 $(0,0,...,0)$ 开始起飞，每次飞翔使得某个坐标增加 $1$，并点亮路过的点，直到飞到终点 $(p_1-1,p_2-1,...,p_n-1)$。

请求出最少需要多少只萤火虫才能点亮所有点。$n\le 32$，$p_i\le 10^9$ 。

转最大反链。以下为方便起见，记 $a_i=p_i-1$ 。

上述问题可以形式化成：给定包含 $n$ 种元素的多重集 $S$，其中 $i$ 种元素出现 $a_i$ 次。定义子集之间的偏序为 $\subseteq$ ，求最大反链。

定理： 记 $M=\lfloor\dfrac{\sum a_i}{2}\rfloor$，则选择所有大小为 $M$ 的子集即达到最大反链。

我们称子集链 $T_1\subseteq T_2\subseteq\cdots\subseteq T_k$ 是对称链，当且仅当 $|T_i|+1=|T_{i+1}|$ 且 $|T_1|+|T_k|=\sum a_i$ 。注意 $k$ 可能为 $1$，此时 $|T_1|={a\over 2}\;(2\mid a)$ 。

引理： 存在对 $S$ 进行 “对称链分解” 的方式，使得 $S$ 的每个子集在恰好一条链中出现。

证明： 归纳构造。当 $n=0$ 时 $T_1=\varnothing$ 是一条链。然后每次加入 $a$ 个元素 $x$，对于原本的对称链 $T_1\subseteq T_2\subseteq\cdots\subseteq T_k$，可将其变为 $\min\{a{+}1,k\}$ 条对称链。

$$\newcommand\ss[0]{\subseteq} \begin{cases} T_1\ss\cdots\ss T_k\ss T_k{+}\{1{\times}x\}\ss\cdots\ss T_k{+}\{a{\times}x\}\\ \\ T_1{+}\{1{\times}x\}\ss\cdots\ss T_{k-1}{+}\{1{\times}x\}\ss T_{k-1}{+}\{2{\times}x\}\ss\cdots\ss T_{k-1}{+}\{a{\times}x\}\\ \\ T_1{+}\{2{\times}x\}\ss\cdots\ss T_{k-2}{+}\{2{\times}x\}\ss T_{k-2}{+}\{3{\times}x\}\ss\cdots\ss T_{k-2}{+}\{a{\times}x\}\\ \\ \cdots\\ \\ T_1{+}\{a{\times}x\}\ss T_2{+}\{a{\times}x\}\ss\cdots\ss T_{k-a}{+}\{a{\times}x\} \end{cases}$$

可以画图来直观感受它。原本的链就是 $T_i=(0,i)$，是 $y$ 轴上的点，而现在 $[0,a]$ 成为 $x$ 轴可选值。这个划分方式就是，先从 $(0,1)$ 开始，进行 “无脑搜索”——往上一直走，直到下一个点已走过，然后往右一直走。再从 $(1,1),(2,1),\dots,(\min\{a,k\},1)$ 分别进行 “无脑搜索”。

不难验证其仍为对称链。$\blacksquare$

注意到，大小为 $M=\lfloor\frac{\sum a_i}{2}\rfloor$ 的子集在每个对称链里都出现了恰好一次。因此，这样的集合的数量就是最长反链的长度（我们构造出了其对应的链覆盖）。

剩下的部分就是折半、容斥，求 $\sum_{i=1}^n x_i=M,x_i\le a_i$ 的整数解的个数。重点其实在于 $\rm Lemma$ 。

$$\begin{aligned} &\sum_{s\subseteq\{1,2,...,n\}}(-1)^{|s|}{M-\sum_{i\in s}p_i+n-1\choose n-1}\\ =&\sum_{i\in S}\sum_{j\in T}\alpha_i\beta_j{M-i-j+n-1\choose n-1}\\ =&\sum_{t=0}^{n-1}\sum_{i\in S}\alpha_i{M-i\choose t}\sum_{j\in T}\beta_j{n-1-j\choose n-1-t} \end{aligned}$$

点击查看代码

ioi2021 集训队互测 子集匹配

设 $\Bbb U=\{1,2,3,\dots,n\}$，记 $\mathcal F_k=\{S\;|\;S\subseteq\Bbb U,\;|S|=k\}$ 。若 $T\in\mathcal F_{k-1},\;S\in\mathcal F_k$ 满足 $T\subseteq S$ 则二者之间连边。请构造大小为 $|\mathcal F_k|$ 的匹配。保证 $2n\geqslant 2k>n$ 。

用对称链构造即可。

[WC2021] 斐波那契

定义 $F_0 = a$，$F_1 = b$，$F_i = (F_{i-1} + F_{i-2}) \bmod m$（$i \ge 2$）。

现在给定 $n$ 组询问，对于每组询问请找到一个最小的整数 $p$，使得 $F_p = 0$ 。

$1 \le n, m \le {10}^5$。

一个初等的思考方法是：转移矩阵的幂可以用斐波那契数填出。

如果互质性比较好，逆元存在，那么直接分离开来变成 $-\dfrac{a}{b}=\dfrac{f_n}{f_{n-1}}$ 。我们继续这个思路（分离 $a,b$ 与 $\{f_n\}$ ）：

考虑移项，并给 $a,b,m$ 同时除掉 $\gcd(a,b,m)$，得到：

$$a'f_{n-1}=(-b')f_n \pmod {m'}$$

此时 $\gcd(a',b',m')=1$，但 $\gcd(a',m'),\gcd(b',m')$ 仍可以 $>1$ 。

由于 $\operatorname{gcd}\left(f_{n-1}, f_{n}\right)=1$ ，此时应有:

$$\begin{array}{l} \operatorname{gcd}\left(f_{n}, m^{\prime}\right)=\operatorname{gcd}\left(a^{\prime}, m^{\prime}\right)=p \\ \\ \operatorname{gcd}\left(f_{n-1}, m^{\prime}\right)=\operatorname{gcd}\left(b^{\prime}, m^{\prime}\right)=q \end{array}$$

（符号问题请大家自觉忽略! 换个元 $b^{\prime \prime}=-b^{\prime}$ 就解决了! )

此时等式两边以及模数同除 $pq$ 即可实现互质，于是就分离开来了!

那么在 $m^{\prime}$ 处塞个三元组 $\left(p, q, \dfrac{f_{n}}{f_{n-1}} \bmod \dfrac{m^{\prime}}{p q}\right)$ ，之后直接查询即可。

斐波那契数的循环节是 $O(n)$ 的。那么该算法预处理的复杂度为 $O(\sigma(m) \times \log m)$ ，其 中 $\sigma(m)$ 是因子和，查询的复杂度是 $O(n \log m)$ 。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n,m;
void exgcd(int a,int b,int &x,int &y){
	if(!b){x=1;y=0;return ;}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
inline int Inv(int x,int md){
	int a,b;exgcd(x,md,a,b);
	return (a%md+md)%md;
}
struct node{
	int x,y,z;
	node(int _x,int _y,int _z){
		x=_x;y=_y;z=_z;
	}
	inline bool operator <(const node &b)const{
		if(x!=b.x)return x<b.x;
		if(y!=b.y)return y<b.y;
		return z<b.z;
	}
};
map<node,int> mp[100005];
int gcd(int x,int y){
	return y?gcd(y,x%y):x;
}
inline void init(){
	for(int i=2;i<=m;i++){
		if(m%i==0){
			int x=1,y=0;
			for(int j=0;;j++){
				if(x&&y){
					int p=gcd(x,i),q=gcd(y,i),m1=i/p/q;
					int k=1ll*(y/q)*Inv(x/p,m1)%m1;
					if(!mp[i].count(node(k,p,q)))mp[i][node(k,p,q)]=j;
				}
				swap(x,y);y=(x+y)%i;
				if(x==1&&y==0)break;
			}
		}
	}
}
int main(){
	scanf("%d%d",&n,&m);init();
	for(int i=1;i<=n;i++){
		int a,b;scanf("%d%d",&a,&b);
		b=(m-b)%m;
		if(a==0)puts("0");
		else if(b==0)puts("1");
		else {
			int d=gcd(gcd(a,b),m),m1=m/d;a/=d,b/=d;
			int p=gcd(a,m1),q=gcd(b,m1),m2=m1/p/q;a/=p,b/=q;
			int k=1ll*a*Inv(b,m2)%m2;
			if(mp[m1].count(node(k,q,p)))printf("%d\n",mp[m1][node(k,q,p)]);
			else puts("-1");
		}
	}


	return 0;
}

[AGC021F] Trinity

现有一个 $N$ 行 $M$ 列的、仅包含黑白格的表格，左上为 $(1,1)$ ，右下为 $(N, M)$ 。 对于一个表格，设长度为 $N$ 的数列 $A$ ，长度为 $M$ 的数列 $B$、$C$ 分别表示：

• $A_{i}$ 表示第 $i$ 行第一个黑格的列号。若不存在则为 $M+1$ 。
• $B_{i}$ 表示第 $i$ 列第一个黑格的行号。若不存在则为 $N+1$ 。
• $C_{i}$ 表示第 $i$ 列最后一个黑格的行号。若不存在则为 $0$ 。

现请你求出所有的 $2^{N M}$ 种表格中，不同的数列三元组 $(A, B, C)$ 的个数对 $998244353$ 取模的结果。

$1 \leq N \leq 8 \times 10^{3}, 1 \leq M \leq 200 \text { 。 }$

我们跳过 $dp$ 推理，得到 $dp$ 转移式如下

$$d p_{m, n}=\left(1+\left(\begin{array}{c} n+1 \\ 2 \end{array}\right)\right) d p_{m-1, n}+\sum_{i<n}\left(\begin{array}{c} n+2 \\ i \end{array}\right) d p_{m-1, i}$$

最后答案为 $\sum_{i=0}^{n}\left(\begin{array}{l}n \\ i\end{array}\right) d p_{m, i}$ 。

自然考虑写成 $\text{EGF}$: 记 $F_{m}(x)=\sum_{n \geq 0} \dfrac{d p_{m, n}}{n!} x^{n}$ 。

那么转移式可以改写成如下的微分方程:

$$F_{m}=e^{x} F_{m-1}+\left(2 e^{x}-x-2\right) F_{m-1}^{\prime}+\left(e^{x}-x-1\right) F_{m-1}^{\prime \prime}$$

之后可以考虑将 $F_{m}$ 写成 $\sum_{i=1}^{m} \sum_{j=1}^{m} f_{m, i, j} x^{i} e^{j x}$ 的形式，这样转移的代价就降到了 $O\left(m^{2}\right)$ 。

这样的时间复杂度其实是 $O\left(m^{3}\right)$ 与 $n$ 无关。

点击查看代码

#include <bits/stdc++.h>
#define mod 998244353
using namespace std;
inline int read(){
    int x=0,f=1;
    char c=getchar();
    while(c<'0'||c>'9')f=(c=='-')?-1:f,c=getchar();
    while(c>='0'&&c<='9')x=x*10+c-'0',c=getchar();
    return x*f;
}
inline int M(int x){return x>=mod?x-mod:x;}
int fsp(int bs,int p){
    int rt=1;
    while(p){
        if(p&1)rt=1ll*rt*bs%mod;
        bs=1ll*bs*bs%mod,p>>=1;
    }
    return rt;
}
inline int C2(int x){return x*(x-1)/2;}
int fac[10004],caf[10004];
inline int C(int x,int y){
	return 1ll*fac[x]*caf[y]%mod*caf[x-y]%mod;
}
inline void init(){
    int lim=10000;
    fac[0]=1;
    for(int i=1;i<=lim;i++)fac[i]=1ll*fac[i-1]*i%mod;
    caf[lim]=fsp(fac[lim],mod-2);
    for(int i=lim;i>=1;i--)caf[i-1]=1ll*caf[i]*i%mod;
}
int rtt[20004],O[20004],ny;
inline void getrtt(int w,int len){
    for(int i=1;i<len;i++)rtt[i]=rtt[i>>1]>>1|((i&1)<<w);
    for(int l=2;l<=len;l<<=1)O[l]=fsp(3,(mod-1)/l);
    ny=fsp(len,mod-2);
}
struct Poly{
    int x[20004];
    int &operator[](int p){return x[p];}
    inline void ntt(int len){
        for(int i=1;i<len;i++)if(rtt[i]>i)swap(x[rtt[i]],x[i]);
        for(int l=2;l<=len;l<<=1){
            for(int i=0,m=l>>1;i<len;i+=l){
                for(int j=i,tO=1,t;j<i+m;j++){
                    t=1ll*tO*x[j+m]%mod,tO=1ll*tO*O[l]%mod;
                    x[j+m]=M(x[j]-t+mod),x[j]=M(x[j]+t);
                }
            }
        }
    }
    inline void idft(int len){
        ntt(len),reverse(x+1,x+len);
        for(int i=0;i<len;i++)x[i]=1ll*x[i]*ny%mod;
    }
}f,g,tf;
int n,m;
int main(){
    n=read();m=read();
    init();
    int w=1,len=2;
    while(len<2*n+1)len<<=1,++w;
    for(int i=1;i<=n;i++)g[i]=caf[i+2];
    getrtt(w-1,len),g.ntt(len);
    tf[0]=f[0]=1;
    for(int i=1;i<=m;i++){
        f.ntt(len);
        for(int j=0;j<len;j++)f[j]=1ll*f[j]*g[j]%mod;
        f.idft(len);
        for(int j=0;j<=n;j++){
            tf[j]=M(1ll*f[j]*fac[j+2]%mod+1ll*tf[j]*(1+j+C2(j))%mod);
            f[j]=1ll*tf[j]*caf[j]%mod;
        }
        for(int j=n+1;j<len;j++)f[j]=0;
    }
    int res=0;
    for(int i=0;i<=n;i++)res=M(res+1ll*C(n,i)*tf[i]%mod);
    printf("%d\n",res);
    return 0;
}

[2019 江苏省队集训] Day1 T1 光影交错

两个人玩游戏，有 $p_{L}$ 的概率 $\text{Alice}$ 获胜，有 $p_{D}$ 的概率 $\text{Bob}$ 获胜，还有 $1-p_{L}-p_{D}$ 的概率两人平局。$0 \leq p_{L}, p_{D}$，$p_{L}+p_{D} \leq 1$ 。

当然两人可以一直玩下去，于是每局游戏后 $p$ 的概率两人停止游戏，$0<p \leq 1$ 。

问游戏停止时，有多少时刻 $\text{Alice}$ 的胜利局数严格多余 $\text{Bob}$ 的胜利局数。

原题容许较大的精度误差，但是我们可以给出其封闭形式!

首先特判 $p=1$ ，以下假设 $p \in(0,1)$ ，另还有 $p_{L}, p_{D}, p_{L}+p_{D} \in[0,1]$ 。

记 $p_{E}=1-p_{L}-p_{D} \in[0,1]$ ，则答案为:

$$\sum_{l>d \geq 0, e \geq 0}\left(\begin{array}{c} l+d+e \\ l, d, e \end{array}\right) p_{L}^{l} p_{D}^{d} p_{E}^{e}(1-p)^{l+d+e-1}$$

作换元 $q_{*} \leftarrow p_{*}(1-p)$ :

$$\begin{aligned} \frac{1}{1-p} \sum_{l>d \geq 0, e \geq 0} q_{L}^{l} q_{D}^{d} q_{E}^{e} &=\frac{1}{1-p} \sum_{l>d \geq 0}\left(\begin{array}{c} l+d \\ l \end{array}\right) q_{L}^{l} q_{D}^{d} \sum_{e \geq 0}\left(\begin{array}{c} l+d+e \\ e \end{array}\right) q_{E}^{e} \\ &=\frac{1}{1-p} \sum_{l>d \geq 0}\left(\begin{array}{c} l+d \\ l \end{array}\right) q_{L}^{l} q_{D}^{d} \frac{1}{\left(1-q_{E}\right)^{l+d+1}} \end{aligned}$$

注意 $p>0$ ，因此 $q_{E}=p_{E}(1-p)<1$ 在收敛域内，再作换元 $t_{*} \leftarrow \dfrac{q_{*}}{1-q_{E}}$ :

$$\frac{1}{(1-p)\left(1-q_{E}\right)} \sum_{l>d \geq 0}\left(\begin{array}{c} l+d \\ l \end{array}\right) t_{L}^{l} t_{D}^{d}$$

之前我的推法是直接处理 $S_{D}(l)=\sum_{d=0}^{l-1}\left(\begin{array}{c}l+d \\ l\end{array}\right) t_{D}^{d}$ ，但总感觉不是很自然。

现在想想，可以尝试作换元 $\delta=l-d-1 \geq 0$ 使得两元独立:

$$\frac{1}{(1-p)\left(1-q_{E}\right)} \sum_{\delta, d \geq 0}\left(\begin{array}{c} 2 d+\delta+1 \\ d \end{array}\right)\left(t_{L} t_{D}\right)^{d} t_{L}^{\delta+1}$$

记 $s=t_{L} t_{D}$ ，则 $\sum_{d \geq 0}\left(\begin{array}{c}2 d+\delta+1 \\ d\end{array}\right) s^{d}=\left(\dfrac{1-\sqrt{1-4 s}}{2 s}\right)^{\delta+1} /(\sqrt{1-4 s})$ ，可以通过拉格朗日反演推得，收敛域为 $[-1 / 4,1 / 4)$ 。

由于 $t_{L}+t_{D}=\dfrac{(1-p)\left(p_{L}+p_{D}\right)}{1-(1-p)\left(1-p_{L}-p_{D}\right)}=\dfrac{(1-p)\left(p_{L}+p_{D}\right)}{(1-p)\left(p_{L}+p_{D}\right)+p}<1$ ，由均值不等式可得 $s=t_{L} t_{D}<1 / 4$ ，又有 $s \geq 0$ 。

但是这个形式并无法兼容 $s=0$ 的情况，我的解决方案是，先假设 $s>0$ ，最后再代入 $s=0$ 来验证形式的统一。所以:

$$\frac{1}{(1-p)\left(1-q_{E}\right)} \sum_{\delta \geq 0} \frac{\left(\frac{1-\sqrt{1-4 s}}{2 s}\right)^{\delta+1}}{\sqrt{1-4 s}} t_{L}^{\delta+1}$$

然后一通分析 $\dfrac{1-\sqrt{1-4 s}}{2 s} t_{L}$ 收敛域 $(0,1)$ 内 (由于条件 $t_{L}+t_{D}<1$ 成立) :

$$\frac{1-\sqrt{1-4 s}}{2 s} t_{L} \times \frac{1}{\sqrt{1-4 s}(1-p)\left(1-q_{E}\right)} \times \frac{1}{1-\frac{1-\sqrt{1-4 s}}{2 s} t_{L}}$$

我们至少得到了一个封闭形式。但很遗恫的是，这个封闭形式对某些边界值末定义。尝试 继续化简:

$$\begin{aligned} & \frac{1}{\sqrt{1-4 s}(1-p)\left(1-q_{E}\right)} \times \frac{1-\sqrt{1-4 s}}{2 t_{D}-1+\sqrt{1-4 s}} \\ =& \frac{1}{\sqrt{1-4 s}(1-p)\left(1-q_{E}\right)} \times \frac{(1-\sqrt{1-4 s})\left(2 t_{D}-1-\sqrt{1-4 s}\right)}{\left(2 t_{D}-1\right)^{2}-(1-4 s)} \\ =& \frac{1}{\sqrt{1-4 s}(1-p)\left(1-q_{E}\right)} \times \frac{1-\sqrt{1-4 s}-2 t_{L}}{2 t_{D}+2 t_{L}-2} \\ =& \frac{1}{\sqrt{1-4 s}(1-p)} \times \frac{\sqrt{1-4 s}+2 t_{L}-1}{2\left(1-q_{E}\right)-2 q_{D}-2 q_{L}} \\ =& \frac{\sqrt{1-4 s}+2 t_{L}-1}{2 \sqrt{1-4 s}(1-p) p} \end{aligned}$$

最后这个形式就比较良好了。当然，也可以把换过的元给代回去，得到完全用题目中的变量表达的式子:

$$\frac{1}{2 p-2 p^{2}}\left(1+\frac{p_{L}-p_{D}-\frac{p}{1-p}}{\sqrt{\left(\frac{p}{1-p}+p_{L}+p_{D}\right)^{2}-4 p_{L} p_{D}}}\right)$$

点击查看代码