Tiw_Air_OAO 数学杂题选讲
伯努利数求自然数幂前缀和
\(\text{EGF}\)（指数型生成函数）形式更为常用，不如说 \(\text{OGF}\)（普通型生成函数）几乎没用。
例：
后面的分式可以拆成 \(\dfrac{x}{e^{x}-1}\)（只有常数项不为零的形式幂级数才可以求逆）与 \(\dfrac{e^{(n+1)x}-1}{x}\) 的乘积。
前者与 \(n\) 无关，可以预先求出。事实上，这就是伯努利数的 \(\text{EGF}\) 。
与斯特林数的异同
斯特林数也可以用于求自然数幂和。
但其实它们是差不多的！回忆斯特林数的 \(\text{EGF}\) ：
那么：
这其中关键的一步就是 \(\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)\)，有如下转移式成立：
仔细瞪一下，作换元 \(R_t=P_t+Q_t\) 与 \(T_t=P_t-Q_t\) 能够给出更简洁的递推式：
发现 \(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}\)，得到：
事实上指标 \(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}}\) 。系数的变化规律是：
\(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\) 轻易描述了：
我不熟悉它，所以我没能指出：\(\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\) 去算即可。
再或者，我们可以换一种刺杀的方式。直接从生成函数的角度考虑它，答案就是
吸取 “\(\text{Binomial Sum}\)” 的经验，只需把 \(F(x)=\dfrac{1}{1+c+x}\) 展开成
我们先算出
Remark：把 \(a\) 去掉后，其实 \(\dfrac{1}{t!}(\text{e}^x{-}1)^t\) 就是第二类斯特林数在第 \(t\) 列上的 \(\rm EGF\) 。因此我们得到的也就是第二类斯特林数在第 \(n\) 行上的 \(\rm OGF\) 。
直接得到
系数 \(\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. 国王奇遇记加强版之再加强版
求：
用 \(m![x^m]e^{ix}=i^m\) 把 \(m\) 放到指数外，有：
其中分母是常多项式，可以分析到答案是关于 \(n+1\) 的多项式。令 \(F(x)=\dfrac{1}{1−me^x}\)，则：
可见，若令：
则：
但是我们对 \(P\)，甚至对 \(f\) 一无所知。怎么办？其实 \(n\) 较小的时候，我们可以用答案最初的式子算出来，继而得到若干 \(P(t)\) 与 \(P(0)\) 的值的线性关系。再者，\(P(x)\) 是 \(m\) 次多项式，其 \(m+1\) 阶差分为 \(0\)。所以：
利用已知线性关系解出 \(P(0)\)，得到点值，\(\text{Lagrange}\) 差值。线性筛自然数幂、线性求一堆逆元，精细实现可以做到 \(O(m)\) 。
#6718. 九个太阳「弱」化版
给定 \(n,k\) ，满足 \(k\) 是 \(2\) 的幂，求
对 \(998244853\) 取模。
发现：
进而 \(\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)\)，则：
需要三次卷积，每次卷积递归计算。复杂度 \(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\}\) 条对称链。
可以画图来直观感受它。原本的链就是 \(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\) 。
点击查看代码
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)\)，得到：
此时 \(\gcd(a',b',m')=1\)，但 \(\gcd(a',m'),\gcd(b',m')\) 仍可以 \(>1\) 。
由于 \(\operatorname{gcd}\left(f_{n-1}, f_{n}\right)=1\) ，此时应有:
（符号问题请大家自觉忽略! 换个元 \(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\) 转移式如下
最后答案为 \(\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}\) 写成 \(\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]\) ，则答案为:
作换元 \(q_{*} \leftarrow p_{*}(1-p)\) :
注意 \(p>0\) ，因此 \(q_{E}=p_{E}(1-p)<1\) 在收敛域内，再作换元 \(t_{*} \leftarrow \dfrac{q_{*}}{1-q_{E}}\) :
之前我的推法是直接处理 \(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\) 使得两元独立:
记 \(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\) 来验证形式的统一。所以:
然后一通分析 \(\dfrac{1-\sqrt{1-4 s}}{2 s} t_{L}\) 收敛域 \((0,1)\) 内 (由于条件 \(t_{L}+t_{D}<1\) 成立) :
我们至少得到了一个封闭形式。但很遗恫的是，这个封闭形式对某些边界值末定义。尝试 继续化简:
最后这个形式就比较良好了。当然，也可以把换过的元给代回去，得到完全用题目中的变量表达的式子:
点击查看代码
