[BZOJ2688]Green Hackenbush

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

[BZOJ2688]Green Hackenbush

题意: 有 $n$ 棵随机的二叉树，每棵只知道大小为 $a_i$

博弈:每次选取一个子树删掉，只剩根不能操作，求先手获胜概率

考虑这个博弈，求出一棵树的 $\text{SG}$ 值

显然有:

1.只有一个点的树的 $\text{SG}$ 值为0

2.多个树组合的问题为 $\text{SG}$ 值异或

暴力 $dp$ ，对于树 $T$ 求答案，设 $T$ 所有可行的后继状态集合为 $N(T)$ ，则得到 $\text{SG}$ 值的表达式为

$\text{SG}(T)=\text{mex}_{R\in N(T)}\lbrace\text{SG(R)}\rbrace$

直接求解复杂度过高，考虑归纳性质

性质:

1.一棵根节点只有一个儿子的树，其 $\text{SG}$ 值为儿子的 $\text{SG}$ 值+1

考虑归纳证明:

设子树为 $T$ ，令 $T+u$ 表示 $T$ 子树上面接上自己作为根，问题变为求证 $\text{SG}(T+u)=\text{SG}(T)+1$

设已经归纳证明所有 $T$ 的子联通块成立

我们要求 $\text{SG}(T+u)$

$\text{SG}(T+u)=\text{mex}\{\text{SG}(u),\forall _{R\in N(T)}\text{SG}(R+u)\}$

由归纳的性质有

$\forall _{R\subsetneq T}\text{SG}(R+T)=\text{SG}(R)+1$

又因为 $\text{SG}(u)=0$ ，看做把所有儿子的情况平移了1，0的位置由自己占据，因而上式成立

2.多叉树的问题可以归纳为 根分别接上每个儿子得到的树 的问题的组合

因为儿子之间实际互不干扰，比较容易理解

由此得到，一棵树的 $\text{SG}$ 值为其所有儿子的 $\text{SG}$ 值+1的异或和

令 $dp_{n,i}$ 为一棵 $n$ 个节点的二叉树 $\text{SG}$ 值为 $i$ 的概率，为了便于转移，设空树的 $\text{SG}$ 值为-1

考虑直接枚举两棵子树的大小和 $\text{SG}$ 值

考虑对于 $n$ 个节点的二叉树，设其左儿子为 $i$ 时的总概率为 $F_i$

得到的 $\text{dp}$ 转移是

$dp_{n,(a+1)\oplus (b+1)}\leftarrow {dp_{i,a}\cdot dp_{n-i-1,b}\cdot F_i}$

我们知道 $n$ 个节点的二叉树方案数为 $Catalan(n)=\frac{(2n)!}{n!(n+1)!}$

由此得到 $\begin{aligned} F_i=\frac{Catalan(i)Catalan(n-i-1)}{Catalan(n)}\end{aligned}$

此题范围可以直接带入 $Catalan(i)$ 求解，但是依然要提一下递推的做法（似乎精度更有保障?）

$\begin{aligned} F_i=\frac{\frac{(2i)!}{i!(i+1)!}\cdot \frac{(2n-i-2)!}{(n-i-1)!(n-i)!}}{\frac{(2n)}{n!(n+1)!}}\end{aligned}$

递推求解 $F_i$ ，每次 $i$ 改变一阶乘只会改变1或者2，因此由 $F_{i-1}$ 得到 $F_i$ 的递推式为

$F_i=\left\{ \begin{aligned}\frac{n(n+1)}{2n(2n-1)}&& i=0\\ F_{i-1}\cdot \frac{2i(2i-1)}{(i+1)i}\frac{(n-i+1)(n-i)}{2(n-i)(2n-2i-1)} && i\in[1,n-1]\end{aligned}\right.$

化简之后应该是

$F_i=\left\{ \begin{aligned}\frac{(n+1)}{2(2n-1)}&& i=0\\ F_{i-1}\cdot \frac{(2i-1)}{(i+1)}\frac{(n-i+1)}{(2n-2i-1)} && i\in[1,n-1]\end{aligned}\right.$

至此得到一个朴素的 $O(n^4)$ 预处理，由于是异或，可以用 $\text{FWT}_{\oplus}$ 求解，复杂度为 $O(n^3)$

对于输入的每棵树，类似背包地叠加概率即可，复杂度为 $O(n^3)$

以下是朴素dp代码

#include<bits/stdc++.h>
using namespace std;
typedef double db;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

char IO;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=128;

int n;
db dp[N][N];

void FWT(db *a,int f){
	for(int i=1;i<N;i<<=1){
		for(int l=0;l<n;l+=i*2) {
			for(int j=l;j<l+i;j++){
				db t=a[j+i];
				a[j+i]=a[j]-t;
				a[j]+=t;
			}
		}
	}
	if(f==-1) rep(i,0,N-1) a[i]/=N;
}

db F[N],G[N];

int main(){
	dp[0][0]=1,dp[1][0]=1;
	rep(i,2,100) {
		F[0]=1.0/(2*i)/(2*i-1)*(i+1)*i;
		rep(j,1,i-1) {
			F[j]=F[j-1] *  (2*j)*(2*j-1)/(j+1)/j   * 1.0/(2*(i-j))/(2*(i-j)-1)*(i-j+1)*(i-j);
		}
		rep(a,0,i-1) rep(h1,0,N-1) if(dp[a][h1]>0) {
			rep(h2,0,N-1) if(dp[i-a-1][h2]) {
				int nxt=0;
				if(a>0) nxt^=h1+1;
				if(i-1-a>0) nxt^=h2+1;
				dp[i][nxt]+=dp[a][h1]*dp[i-a-1][h2]*F[a];
			}
		}
	}
	n=rd();
	rep(i,0,N-1) F[i]=0;
	F[0]=1;
	rep(i,1,n) {
		int x=rd();
		rep(j,0,N-1) G[j]=0;
		rep(j,0,N-1) if(F[j]) rep(k,0,N-1) G[j^k]+=F[j]*dp[x][k];
		rep(j,0,N-1) F[j]=G[j];
	}
	db ans=0;
	rep(i,1,N-1) ans+=F[i];
	printf("%.6lf\n",ans);
}

以下是FWT优化代码

#include<bits/stdc++.h>
using namespace std;
typedef double db;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

char IO;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=128;

int n;
db dp[N][N],T[N][N];

void FWT(db *a,int f){
	for(int i=1;i<N;i<<=1){
		for(int l=0;l<N;l+=i*2) {
			for(int j=l;j<l+i;j++){
				db t=a[j+i];
				a[j+i]=a[j]-t;
				a[j]+=t;
			}
		}
	}
	if(f==-1) rep(i,0,N-1) a[i]/=N;
}

db F[N],G[N];

int main(){
	dp[0][0]=1,dp[1][0]=1;
	T[0][0]=1,T[1][1]=1;
	FWT(T[0],1),FWT(T[1],1);

	rep(i,2,100) {
		F[0]=1.0/(2*i)/(2*i-1)*(i+1)*i;
		rep(j,1,i-1) {
			F[j]=F[j-1] *  (2*j)*(2*j-1)/(j+1)/j   * 1.0/(2*(i-j))/(2*(i-j)-1)*(i-j+1)*(i-j);
		}
		rep(j,0,i-1) rep(k,0,N-1) dp[i][k]+=T[j][k]*T[i-j-1][k]*F[j];
		FWT(dp[i],-1);
		rep(j,0,N-2) T[i][j+1]=dp[i][j];
		FWT(T[i],1);
	}
	n=rd();
	rep(i,0,N-1) F[i]=0;
	F[0]=1;
	rep(i,1,n) {
		int x=rd();
		rep(j,0,N-1) G[j]=0;
		rep(j,0,N-1) if(F[j]) rep(k,0,N-1) G[j^k]+=F[j]*dp[x][k];
		rep(j,0,N-1) F[j]=G[j];
	}
	db ans=0;
	rep(i,1,N-1) ans+=F[i];
	printf("%.6lf\n",ans);
}

「2019五校联考-雅礼」大凯的疑惑

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「2019五校联考-雅礼」大凯的疑惑

首先判断是否有无穷解，即判断 $\gcd(a,b)>1$ 时有无穷解

接下来我们由小凯的疑惑知道最大的无法表示的数是 $ab-a-b$ ，这能确定一个上界

考虑计算 $[1,R](R<ab)$ 中能用 $a,b$ 表示出来的数

因为 $\gcd(a,b)=1,R<ab$ ，所以每个数最多只有一种构成法，可以枚举其包含了几个 $b$ ，剩下的部分直接任意放 $a$

即得到计算能够被构成的个数的方法为 $\begin{aligned}\sum_{i=0}^{\lfloor \frac{R}{a}\rfloor} \lfloor \frac{R-ib}{a}\rfloor+1 \end{aligned}$

其中+1是计算了包含0个 $a$ 的情况

如果二分答案，复杂度为 $O(a\log (ab))$ ，恐怕难以通过

优化：

我的思路是是先确定了 $\lfloor \frac{R}{a}\rfloor$ ，那么此时确定了所有 $b$ 的个数的贡献

那么考虑枚举，找到答案所属的 $\lfloor \frac{R}{a}\rfloor$ 的区间，在这段区间里，判断一个数 $x$ 是否可以被构成即:

$x\equiv ib\pmod a(i\leq \lfloor \frac{R}{a}\rfloor)$ ，即考虑了不同 $b$ 的个数的贡献

用一个数组存下 $ib\bmod a$ ，那么可以 $O(1)$ 判断一个数是否合法，如果直接for过去是 $O(b)$ 的

显然这在一段中，构成情况每 $a$ 个一循环，那么先快速跳循环即可

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

const int N=2e7+10;

ll a,b,k;
ll gcd(ll a,ll b){ return b==0?a:gcd(b,a%b); }
bool mk[N];
int main() {
	freopen("math.in","r",stdin),freopen("math.out","w",stdout);
	scanf("%lld%lld%lld",&a,&b,&k);
	if(gcd(a,b)!=1 || a==1) return puts("-1"),0;
	if(a>b) swap(a,b);
	if(k==1) return printf("%lld\n",a*b-a-b),0;
	ll s=(a-1)*(b-1)/2,c=0; // s为总的个数
	if(s<k) return puts("-1"),0;
	k=s-k+1; // 改为求第k小
	int p=mk[0]=1;  // p为所属区间
	rep(i,1,a-1) {
		c+=i*b/a;
		ll t=i*(b-1)-c; // t为这段区间内无法被表示的个数
		if(t>=k){ p=i; break; }
		mk[i*b%a]=1;  // 把区间内的ib mod a 放进去
	}
	k-=(p-1)*(b-1)-(c-p*b/a); //还需要做的个数
	ll l=(p-1)*b+1,d=l%a; //l为区间开始位置
	ll i=l+(k-1)/(a-p)*a; // 每个长度为a的循环中已经有p个位置被标记，可以被表示，因此还有a-p个位置无法表示
	k-=(k-1)/(a-p)*(a-p); // 跳过循环
	for(;;++i,d++,d==a&&(d=0)) if(!mk[d] && --k==0) return printf("%lld\n",i),0; // 暴力for最后a个
}

「2020-2021 集训队作业」Yet Another Permutation Problem

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「2020-2021 集训队作业」Yet Another Permutation Problem

题目大意

对于一个初始为 $1,2,\ldots n$ 的排列，每次操作为选择一个数放到开头或者结尾，求 $k$ 次操作能够生成的排列数

对于 $k=0,1,\ldots ,n-1$ 求解

\[ \ \]

模型转化

容易发现，对于一个排列，生成它的最小次数取决于中间保留段的长度

而保留段实际上是任何一个上升子段

设一个排列的最长上升子段为 $l$ ，那么最少操作步骤就是 $n-l$

那么对于 $k$ ，合法的序列就是存在一个长度 $\ge n-k$ 的上升子段

存在不好算，改为计算任何一个上升子段 $<n-k$ 的数量

为了便于描述，令下文的 $k=n-k-1$

\[ \ \]

生成函数构造

考虑一个序列是由若干上升段构成的，设一个长度为 $l$ 的上升段的权值为 $[l\leq k]$

那么排列的权值就是上升段权值之积

容易想到用 $\text{EGF}$ 合并上升段，但是直接的统计，我们无法保证上升段之间无法拼接

假设我们确定了一个单位上升段的 $\text{EGF}$ 为 $G(x)$ ， $\text{OGF}$ 为 $F(x)$

那么按照上面 $\text{Naive}$ 的计算，上升段之间的合并为有序拼接，即 $\displaystyle \sum_{i=0}G^i(x)=\frac{1}{1-G(x)}$

容易发现，这样的计算，会导致一个长度为 $l$ 的极长上升段被分解成若干小段

也就是被计算了 $\displaystyle [x^l](\sum_{i=0}F^i(x))=[x^l]\frac{1}{1-F(x)}$ 次

在合法的计算中，我们希望， $[x^l]\frac{1}{1-F(x)}$ 恰好为权值 $[l\leq k]$

也就是说，我们希望 $\displaystyle \frac{1}{1-F(x)}=H(x)=\sum_{i=0}^kx^i=\frac{x^{k+1}-1}{x-1}$

那么可以反向由 $H(x)$ 构造出我们想要的 $F(x)$ ，从而得到 $G(x)$ ，再进行求解

\[ \ \]

答案计算

$\displaystyle F(x)=1-\frac{1}{H(x)}=1-\frac{x-1}{x^{k+1}-1}=\frac{x-x^{k+1}}{1-x^{k+1}}$

可以爆算得到 $F(x)$ ，从而得到 $G(x)$ ，然后暴力求逆就是 $O(n^2)$

优化：

$1-x^{k+1}$ 的逆，只包含 $\frac{n}{k+1}$ 项，所以 $G(x)$ 只含 $2\frac{n}{k+1}$ 项

即 $\displaystyle F(x)=\sum_{d=0}x^{d(k+1)+1}-\sum_{d=1}x^{d(k+1)}$ ， $G(x)$ 就是除一个阶乘

这样暴力求逆就是 $O(n^2\ln n)$

~~（不是你干嘛要真的求逆，直接进行 $G(x)$ 的叠加就可以了）~~

const int N=1010;

int n,P,I[N],J[N];
ll qpow(ll x,ll k=P-2){
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int F[N];

int main(){
	n=rd(),P=rd();
	rep(i,*J=1,n) J[i]=1ll*J[i-1]*i%P;
	I[n]=qpow(J[n]);
	drep(i,n,1) I[i-1]=1ll*I[i]*i%P;
	drep(k,n,1) {
		F[0]=1;
		rep(j,1,n) {
			F[j]=0;
			for(int d=1;d<=j;d+=k) F[j]=(F[j]+1ll*F[j-d]*I[d])%P;
			for(int d=k;d<=j;d+=k) F[j]=(F[j]-1ll*F[j-d]*I[d])%P;
		}
		printf("%d\n",int((1ll*(P+1-F[n])*J[n])%P));
	}
}

「2020noip模拟题-蒋凌宇」幂

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「2020noip模拟题-蒋凌宇」幂

Analysis

计算 $x$ 出现的次数，可以转化为枚举每一个 $x$ ，计算剩余 $n-1$ 个位置合法括号序列个数

因此我们只需要计算合法括号序列个数

定义一个辅助计数：不可分割 的合法括号序列

这样的括号序列，满足：其恰好为 $x$ ，或者序列两端是一对匹配的左右括号

而实际要求的括号序列的就是这样的 不可分割括号序列去掉两端的匹配括号

\[\ \]

Solution

设 $a_n$ 为长度为 $n$ 的不可分割 的合法括号序列数量

括号序列并非排列问题，因此我们用普通生成函数计算

设 $\text{OGF(a)}=A(x)$ ， $A(x)$ 容易发现 $A(x)$ 满足下面的递归式

$\displaystyle A(x)=x^2(\sum_{i=0}^{\infty} A^i(x))+x^1$

其中 $x^2$ 表示在外层加一对匹配括号， $\sum_{i=0}^{\infty} A^i(x)$ 枚举子括号中的分裂段， $x^1$ 表示单个 $x$

容易得到下面的化简过程

$\displaystyle A(x)=\frac{x^2}{1-A(x)}+x$

$A(x)-A^2(x)=x^2+xA(x)$

$A^2(x)-(x+1)A(x)+x^2+x=0$

带入求根公式，得到 $A(x)$ 的 收敛形式

$\displaystyle A_1(x)=\frac{x+1+\sqrt{-3x^2-2x+1}}{2},A_2(x)=\frac{x+1-\sqrt{-3x^2-2x+1}}{2}$

令 $\displaystyle F(x)=-3x^2-2x+1,G(x)=\sqrt{F(x)}$

容易手玩发现： $[x^0]G(x)=1$

而根据定义，我们知道 $[x^0]A(x)=0$ ，因此 $A(x)=A_2(x)$

接下来我们要求 $\displaystyle G(x)=F^\frac{1}{2}(x)$

\[ \ \]

下面介绍对于短多项式 $F(x)$ ，设 $\text{deg}(F(x))=m$ ，有理数 $k(k\ne 1)$

求解 $G(x)=F^k(x)$ 的前 $n$ 项的 $O(m^2+nm)$ 递推做法

变形

$G(x)=F^k(x)$

$\displaystyle G'(x)=kF^{k-1}(x)F'(x)$

$G'(x)F(x)=kF^k(x)F'(x)$

$G'(x)F(x)=kG(x)F'(x)$

\[ \ \]

求解递推式

对于等号两边，考虑 $[x^n]$ 一项的系数，容易求出 $F'(x)=-6x-2$

$\displaystyle \sum_{i=0}^m [x^{n-i}]G'(x)F_i=k\sum_{i=0}^{m-1}[x^{n-i}]G(x)F'_i$

$\displaystyle \sum_{i=0}^m [x^{n-i+1}](n-i+1)G(x)F_i=k\sum_{i=0}^{m-1}[x^{n-i}]G(x)F'_i$

$\displaystyle (n+1)[x^{n+1}]G(x)=k\sum_{i=0}^{m-1}[x^{n-i}]G(x)F'_i-\sum_{i=1}^m [x^{n-i+1}](n-i+1)G(x)F_i$

带入这题的 $k$ ，得到

$\displaystyle [x^n]=\frac{3(n-3)[x^{n-2}]+(2n-3)[x^{n-1}]}{n}$

递推边界 $[x^0]G(x)=1,[x^1]G(x)=-1$

然后由 $G(x)$ 得到 $A(x)$ 再得到最终答案即可

「2020联考北附1」命运歧途

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「2020联考北附1」命运歧途

排列 $dp$ 问题通常想到容斥，因为很难在 $dp$ 的同时保证排列元素不多次出现

对于这个问题，我们只需要考虑相邻关系

我们需要计算包含0个非法位置的排列个数，因此容斥容易定义为计算

包含至少 $i$ 个非法位置的排列个数，容斥系数可以简单设置为 $(-1)^i$

考虑一个序列的合法与非法情况，容易发现是类似下面的情况

序列会被分成若干段，每一段是形如 $x,x+k,x+2k,\cdots$ 或者 $x,x-k,x-2k$

显然 $x,y$ 会产生非法关系的必要条件是 $x\equiv t \pmod k$ ，因此按照 $x\mod k$ 分组

组内实际上是类似 $k=1$ 的子问题，不妨设一组包含 $m$ 个元素

显然不同组之间一定构成不同的段，而对于一组内的元素，我们可以强制某一些位置分段

最终会得到一个若干段的序列，段之间由于不确定相互关系，设最终分成 $i$ 段，则可以有 $i!$ 种段之间排列

(分成 $i$ 段的贡献在我们计算时实际上是至少包含了 $n-i$ 个非法位置)

现在问题就是落到了 $dp$ 序列分成 $i$ 段的方案数，显然对于每一组，是一个分组背包的问题

对于每一个大小为 $m$ 的组，考虑 $dp$ 其分段方案

令 $dp_{i,j}$ 表示当前已经确定的总大小为 $i$ 且已经分成了 $j$ 段

转移可以枚举下一段的大小 $k$ ，由于实际排列中的段是有递增和递减两种情况，因此对于任意 $k>1$ ，有两种段分配方法

转移式子是

$dp_{i,j}\rightarrow dp_{i+k,j+1} (k=1)$

$2dp_{i,j}\rightarrow dp_{i+k,j+1} (k>1)$

这个式子容易用前缀和优化，特殊的转移位置只有一个，处理一下即可

于是可以在 $O(n^2)$ 复杂度内计算任何一个大小的组的分段方案

暴力合并组之间的背包，对于每个查询，复杂度为 $O(n^2)$

因此总复杂度受限于查询，为 $O(qn^2)$

接下来考虑如何削减查询复杂度

由于 $q$ 的上限实际是 $n$ ，下文认为 $q=n$

优化1

上面的问题中，我们并没有具体地分析分组的情况

实际上对于 $n,k$ ，可能的组大小只有两种，即 $\lfloor \frac{n}{k}\rfloor ,\lceil \frac{n}{k}\rceil$

不同的 $\lfloor \frac{n}{k}\rfloor$ 只有 $O(\sqrt n)$ 种，不妨对于每种 $\lfloor \frac{n}{k}\rfloor$ 从小到大计算每一个 $k$ 的答案

简单观察即可发现：

在每个块内，按照 $k$ 递增，两种组的个数一个递增一个递减

如果在每个 $\lfloor \frac{n}{k}\rfloor$ 中，为最小的 $k$ 暴力 $O(n^2)$ 预处理出答案，然后不断增大

不妨维护两个单调的个数指针

如果能够支持背包回撤一个分组，增加一个分组，那么容易做到每个块内 $O(n^2)$ 递推答案

增加一个分组不必说，而对于去掉一个分组，不好求出模逆元

但是由于上面的 $dp_{i,j}$ 满足 $dp_{i,i}=1$ ，因此考虑从高到低递推每一位

模拟一个长除法即可

合理的常数优化也可以通过此题

const int N=2010;

int n,m,P;
int C[N][N],J[N];
int dp[N][N],S[N],F[N],G[N];

int len;
int pl=-1,nc1,nc2;
void clear() { 
	rep(i,0,len) F[i]=0;
	F[len=0]=1,nc1=nc2=0; 
}
int cnt=0;
void Mul(int m,int *A){ 
	drep(i,len+=m,0) {
		ull t=0;
		rep(j,1,min(i,m)) {
			t+=1ll*F[i-j]*A[j];
			((j&15)==0) && (t%=P);
		}
		F[i]=t%P;
	}
}
void Div(int m,int *A){ 
	static int G[N];
	rep(i,0,len) G[i]=F[i],F[i]=0;
	drep(i,len-=m,0) {
		ull t=G[i+m];
		rep(j,1,m) {
			t+=1ll*F[i+j]*(P-A[m-j]);
			((j&15)==0) && (t%=P);
		}
		F[i]=t%P;
	}
}

int Ans[N],Pow1[N],Pow2[N];
int clr;
void Solve(int m) {
	int l=(n-1)/m+1,c1=0,c2=0;
	rep(i,1,m) if((n-i)/m+1==l) c1++;
	else c2++;
	if(l==2) {
		rep(i,0,c1+c2) F[i]=0;
		rep(i,0,n-c1-c2) F[i+c1+c2]=1ll*C[c1][i]*Pow1[c1-i]%P*Pow2[i]%P;
		int x=1;
		rep(i,1,c2) x=1ll*x*dp[1][1]%P;
		rep(i,c1+c2,n) F[i]=1ll*F[i]*x%P;
	} else {
		if(l!=pl) {
			clr++;
			clear(); 
			pl=l;
		}
		cnt++;
		while(c1<nc1) Div(l,dp[l]),nc1--;
		while(c2<nc2) Div(l-1,dp[l-1]),nc2--;
		while(c1>nc1) Mul(l,dp[l]),nc1++;
		while(c2>nc2) Mul(l-1,dp[l-1]),nc2++;
	}
	rep(i,0,n) G[n-i]=1ll*F[i]*J[i]%P;
	int ans=0;
	rep(i,0,n) ans=(ans+1ll*((i&1)?-1:1)*G[i])%P;
	ans=(ans%P+P)%P;
	Ans[m]=ans;
}

int main() {
	freopen("fate.in","r",stdin),freopen("fate.out","w",stdout);
	n=rd(),m=rd(),P=rd();
	rep(i,0,n) rep(j,C[i][0]=1,i) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
	rep(i,J[0]=1,n) J[i]=1ll*J[i-1]*i%P;
	dp[0][0]=S[0]=1;
	rep(i,1,n) {
		drep(j,i,1) {
			dp[i][j]=(2ll*S[j-1]+P-dp[i-1][j-1])%P;
			S[j]+=dp[i][j],Mod1(S[j]);
		}
	}
	rep(i,Pow1[0]=1,n) Pow1[i]=1ll*Pow1[i-1]*dp[2][1]%P;
	rep(i,Pow2[0]=1,n) Pow2[i]=1ll*Pow2[i-1]*dp[2][2]%P;
	rep(i,1,n-1) Solve(i);
	Ans[n]=J[n];
	fprintf(stderr,"%d %d\n",cnt,clr);
	rep(i,1,m) printf("%d\n",Ans[rd()]);
}

\[ \ \]

优化2

从生成函数的角度，我们容易把所求的值归纳为

$H(x)=F^a(x)G^b(x)$

我们知道多项式快速幂的复杂度为 $\exp$ 的 $n\log n$

当然那是对于长度为 $n$ 的情况

而现在是长度之和为 $n$

由于EI在无数道题中介绍了这个东西，看到马上想到可以求导

$H'(x)=aF^{a-1}(x)F'(x)G(x)+bG^{b-1}(x)G'(x)F(x)$

$H'=H(a\cdot \frac{F'}{F}+b\cdot \frac{G'}{G})$

理想的情况是：对于这个式子，两边从低到高解方程确定每一项的值

然而首先遇到的多项式除法操作就难以解决

要除掉的 $F,G$ 没有常数项，而第一项系数为1或2，因此除法操作只需要计算 $2$ 的逆元

因此可以考虑对于模数中的 $2^t$ 提取出来，然后最后暴力exCRT合并

接下来就是分成两部分

计算 $\mod 2^t$

由上面知道，背包的每一个位置实际对于最后答案的贡献是 $i!(-1)^{n-i}$

因此对于 $i!\mod 2^t=0$ 的位置都不用考虑，只需要求出背包前 $O(\log P)$ 位，暴力即可

\[ \ \]

计算 $\mod P(2\not |P)$

接下来就是上面的递推过程

然而还有一个问题就是从 $[x^n]H'$ 得到 $[x^{n+1}]H$ ，显然这里需要一个逆元

EI为我们提供一个很好的思路，或许有助于解决任意模数的逆元问题：

因为答案计算的是 $k![x^k]H$ ，因此可以直接在计算过程中加入

这样求导的系数在阶乘中被省略，变成了单纯的平移

而乘法只需要额外添加一个权值 $\binom{i+j}{i}\cdot x^i\cdot x^j\rightarrow x^{i+j}$

接下来具体的方法是：

先将 $F(x),G(x)$ 平移一位去掉空余的项，处理过后 $H(x)$ 被平移了 $a+b$ 的位置，首项为 $x^0$

从低到高递推 $H(x)$ 的每一位，同步维护 $\frac{H}{F},\frac{H}{G}$

对于 $H'(x)$ 的第 $i$ 项累和得到，然后平移一位得到 $H(x)$ 的 $i+1$ 项

最后需要加上平移部分的贡献， $i!\rightarrow (i+a+b)!$ ，可以用一个组合数解决

const int N=2010;

int n,m,P;
int J[N];
int dp[N][N],S[N],F[N],G[N];

struct Solve2t{
	int F[N],P,U;
	void Init(int x) {
		P=x;
		int t=1;
		rep(i,1,n) {
			t=1ll*t*i%P;
			if(t==0) break;
			U=i;
		}
	}
	int Solve(int m) {
		rep(i,0,U) F[i]=0;
		F[0]=1;
		rep(x,1,m) {
			int c=(n-x)/m+1;
			rep(i,0,U) G[i]=F[i],F[i]=0;
			rep(i,0,U) if(G[i]) rep(j,1,min(c,U)) F[i+j]=(F[i+j]+1ll*G[i]*dp[c][j])%P;
		}
		int ans=0;
		rep(i,1,U) {
			F[i]=1ll*F[i]*J[i]%P;
			ans=(ans+(((n-i)&1)?-1:1)*F[i])%P;
		}
		ans=(ans%P+P)%P;
		return ans;
	}
} Sol1;

void Exgcd(ll a,ll b,ll &x,ll &y) {
	if(b==0) {
		x=1,y=0;
		return;
	}
	Exgcd(b,a%b,y,x),y-=a/b*x;
}
ll Inv(int a,int P) {
	ll x,y;
	Exgcd(a,P,x,y);
	return (x%P+P)%P;
}

struct SolveP{
	int C[N][N];
	int F[N],G[N];
	int A[N],B[N];
	int X[N],Y[N];
	int P,I2;
	void Init(int x) { 
		P=x; 
		rep(i,0,n) rep(j,C[i][0]=1,i) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
		I2=(P+1)/2;
	}
	int Solve(int m) {
		if(P==1) return 0;
		int l1=(n-1)/m+1,l2=l1-1;
		int c1=0,c2=0;
		rep(i,1,m) if((n-i)/m+1==l1) c1++;
		else c2++;
		int Inv=dp[l1][1]==1?1:I2;
		rep(i,1,l1) A[i-1]=1ll*J[i-1]*dp[l1][i]%P*Inv%P;
		Inv=dp[l2][1]==1?1:I2;
		rep(i,1,l2) B[i-1]=1ll*J[i-1]*dp[l2][i]%P*Inv%P;
		// 将A,B偏移补充常数项
		// 同时除掉常数项，为下面做除法铺路

		int d=c1+c2;
		// A,B存储偏移之后的值，d为偏移量
		// X,Y存储除法之后的值
		F[0]=X[0]=Y[0]=1;
		rep(i,1,n-d) {
			ull x=0,y=0;
			// 注意这里计算的是两边次数为x^{i-1}一项的值
			rep(j,1,min(l1-1,i)) {
				// 将除法的结果乘上求导以后的值，注意求导以后第一位消失了
				x+=1ll*X[i-j]*A[j]%P*C[i-1][j-1];
				(j&15)==0 && (x%=P);
				// 求导以后是j-1位，因此实际应该是x^{i-j}*x^{j-1}
			}
			rep(j,1,min(l2-1,i)) {
				y+=1ll*Y[i-j]*B[j]%P*C[i-1][j-1];
				(j&15)==0 && (y%=P);
			}
			// 由[x^{i-1}]H'(x)一项得到[x^i]H(x)
			F[i]=X[i]=Y[i]=(x%P*c1+y%P*c2)%P;
			// 接下来做二项除法(雾)
			x=0,y=0;
			rep(j,1,min(i,l1-1)) {
				x+=1ll*(P-A[j])*X[i-j]%P*C[i][j];
				(j&15)==0 && (x%=P);
			}
			rep(j,1,min(i,l2-1)) {
				y+=1ll*(P-B[j])*Y[i-j]%P*C[i][j];
				(j&15)==0 && (y%=P);
			}
			X[i]=(X[i]+x)%P,Y[i]=(Y[i]+y)%P;
		}
		int ans=0;
		rep(i,0,n-d) {
			F[i]=1ll*F[i]*C[i+d][d]%P*J[d]%P;
			ans=(ans+(((n-i-d)&1)?-1:1)*F[i])%P;
		}
        // 将除掉的系数乘回来
		rep(i,1,c1) ans=1ll*ans*dp[l1][1]%P;
		rep(i,1,c2) ans=1ll*ans*dp[l2][1]%P;
		ans=(ans%P+P)%P;
		return ans;
	}
} Sol2;

int t;
void Init() {
	t=1;
	for(t=1;P%2==0;) t*=2,P/=2;
	Sol1.Init(t),Sol2.Init(P);
}

void Solve(int m) {
	int x=Sol1.Solve(m),y=Sol2.Solve(m);
    // CRT
	ll res=1ll*((y-x)%P+P)*Inv(t,P)%P;
	ll mod=t*P;
	res=((res*t+x)%mod+mod)%mod;
	printf("%lld\n",res);
}

int main() {
	freopen("fate.in","r",stdin),freopen("fate.out","w",stdout);
	n=rd(),m=rd(),P=rd();
	rep(i,J[0]=1,n) J[i]=1ll*J[i-1]*i%P;
	dp[0][0]=S[0]=1;
	rep(i,1,n) {
		drep(j,i,1) {
			dp[i][j]=(S[j-1]*2ll+P-dp[i-1][j-1])%P;
			S[j]+=dp[i][j],Mod1(S[j]);
		}
	}
	Init();
	rep(i,1,m) Solve(rd());
}

「2019 集训队互测 Day 1」最短路径 (点分治+NTT/FFT+线段树)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「2019 集训队互测 Day 1」最短路径 (点分治+NTT/FFT+线段树)

题意：给定了一棵基环树，求所有的 $d(u,v)^k$ 的期望

当 $k$ 较小时，可以想到用斯特林数/二项式定理展开维护+1操作，对于树的可以从儿子合并上来，对于环上可以枚举每个块求得答案

复杂度为 $O(nk)$

当图为一棵树时，由于不好处理 $x^k$ ，考虑直接求出 $d(u,v)=i$ 的数量

比较容易想到用用点分治+ $\text{NTT}$ 求解，复杂度为 $O(n\log ^2n)$

环上的情况比较麻烦，不妨为每个块标号 $1,2,\cdots m$ ，每个块包含 $sz_i$ 个结点

显然 $(i,j)$ 的距离为 $\min\lbrace|i-j|,m-|i-j|\rbrace$

考虑计算所有块 $(i,j)(i<j)$ 之间的贡献，令 $d=\lfloor \frac{m}{2}\rfloor$ ，则对于 $j\in[i+1,i+d]$ 在环上的距离为 $j-i$ ，否则距离为 $m-(j-i)$

对于两种情况分类讨论，这里以计算 $j\in[i+1,i+d]$ 为例

因为是一段区间，考虑直接在线段树的 $[i+1,i+d]$ 加入 $i$ ，然后对于线段树上每个结点计算

推论1：能够被添加到线段树结点 $[l,r]$ 上的 $i$ 构成一段连续的区间

推论2：从区间 $[l,r]$ 的一端出发， $\text{dfs}$ 区间内的块得到的 $\max dis_u\leq \sum_{i=l}^r sz_i$

因此同样考虑用 $\text{NTT}$ 维护该答案，每次更新答案可以看做是区间 $[l1,r1],[l2,r2](r1<l2)$ 之间的贡献

分别从 $r1,l2$ 开始 $\text{dfs}$ 得到 $dis_u$ ，然后 $\text{NTT}$ 合并，不把 $[r1+1,l2-1]$ 这一部分在环上的加入 $\text{NTT}$ 大小

这样就能保证卷积大小 $\leq \sum_{i=l1}^{r1} sz_i+\sum_{i=l2}^{r2} sz_i$

同理可以类似处理 $j>i+d$ 的情况

分析复杂度：每个 $i$ 会出现在线段树上 $\log n$ 个位置，每个 $j$ 会在线段树上 $\log n$ 层被计算

因此每个点被加入卷积大小的次数为 $O(\log n)$ ，复杂度为 $O(n\log ^2 n)$ 与前面的点分治同阶

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

bool Mbe;
const int N=1<<18|10,P=998244353;

int n,m,k;
int A[N];
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int Pow[N];
struct Edge{
	int to,nxt;
}e[N];
int head[N],ecnt,deg[N];
void AddEdge(int u,int v) {
	e[++ecnt]=(Edge){v,head[u]};
	head[u]=ecnt,deg[v]++;
}
#define erep(u,i) for(int i=head[u];i;i=e[i].nxt)

int w[N];
void Init() {
	int R=1<<18;
	int t=qpow(3,(P-1)/R);
	w[R/2]=1;
	rep(i,R/2+1,R-1) w[i]=1ll*w[i-1]*t%P;
	drep(i,R/2-1,1) w[i]=w[i<<1];
}

int rev[N];
void NTT(int n,int *a,int f) {
	static int e[N>>1];
	rep(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=e[0]=1,t;i<n;i<<=1) {
		int *e=w+i;
		for(int l=0;l<n;l+=i*2) {
			for(int j=l;j<l+i;++j) {
				t=1ll*a[j+i]*e[j-l]%P;
				a[j+i]=a[j]-t,Mod2(a[j+i]);
				a[j]+=t,Mod1(a[j]);
			}
		}
	}
	if(f==-1) {
		reverse(a+1,a+n);
		ll base=qpow(n);
		rep(i,0,n-1) a[i]=a[i]*base%P;
	}
}
int Init(int n) {
	int R=1,c=-1;
	while(R<=n) R<<=1,c++;
	rep(i,0,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
	return R;
}

int Q[N],L,R,vis[N];

namespace pt1{ 
	const int N=1010;
	int dis[N];
	void Bfs(int u) {
		rep(i,1,n) dis[i]=-1;
		dis[Q[L=R=1]=u]=0;
		while(L<=R) {
			u=Q[L++];
			erep(u,i){
				int v=e[i].to;
				if(~dis[v]) continue;
				dis[v]=dis[u]+1,Q[++R]=v;
			}
		}
	}
	void Solve() {
		int ans=0;
		rep(i,2,n) {
			Bfs(i);
			rep(j,1,i-1) ans=(ans+Pow[dis[j]])%P;
		}
		ans=ans*qpow(n*(n-1)/2)%P;
		printf("%d\n",ans);
	}
}

int Ans[N],sz[N];
namespace pt2{ 
	int mi=1e9,rt;
	void FindRt(int n,int u,int f) {
		int ma=0; sz[u]=1;
		erep(u,i) {
			int v=e[i].to;
			if(v==u || v==f || vis[v]) continue;
			FindRt(n,v,u),sz[u]+=sz[v],cmax(ma,sz[v]);
		}
		cmax(ma,n-sz[u]);
		if(mi>ma) mi=ma,rt=u;
	}

	int F[N],A[N],B[N];
	void Solve(int n,int k) {
        // 容斥型 点分治
		int R=Init(n*2+1);
		rep(i,0,R) F[i]=0;
		rep(i,0,n) F[i]=A[i];
		NTT(R,F,1);
		rep(i,0,R-1) F[i]=1ll*F[i]*F[i]%P;
		NTT(R,F,-1);
		if(k==1) rep(i,0,n*2) Ans[i]+=F[i],Mod1(Ans[i]);
		else rep(i,0,n*2) Ans[i]-=F[i],Mod2(Ans[i]);
	}
	int maxd;
	void dfs(int u,int f,int d=0) {
		A[d]++,sz[u]=1,cmax(maxd,d);
		erep(u,i) {
			int v=e[i].to;
			if(v==u || v==f || vis[v]) continue;
			dfs(v,u,d+1),sz[u]+=sz[v];
		}
	}
	void Divide(int n,int u) {
		mi=1e9,FindRt(n,u,0),u=rt;
		vis[u]=1;
		int D=0;B[0]=1;
		erep(u,i) {
			int v=e[i].to;
			if(vis[v]) continue;
			maxd=0,dfs(v,u,1);
			Solve(maxd,-1);
			rep(j,0,maxd) B[j]+=A[j],A[j]=0;
			cmax(D,maxd);
		}
		rep(i,0,D) A[i]=B[i],B[i]=0;
		Solve(D,1);
		rep(i,0,D) A[i]=0;
		erep(u,i) {
			int v=e[i].to;
			if(vis[v]) continue;
			Divide(sz[v],v);
		}
	}
	void Solve() {
		rep(i,1,n) vis[i]=0;
		Divide(n,1);
		int ans=0;
		rep(i,1,n) ans=(ans+1ll*Ans[i]*Pow[i])%P;
		ans=ans*qpow(1ll*n*(n-1)%P)%P;
		printf("%d\n",ans);
	}
}

int QL[N<<2],QR[N<<2];
void Add(int p,int l,int r,int ql,int qr,int x) {
    // 在线段树上加入结点
	if(ql<=l && r<=qr) {
		if(!QL[p]) QL[p]=x;
		QR[p]=x;
		return;
	}
	int mid=(l+r)>>1;
	if(ql<=mid) Add(p<<1,l,mid,ql,qr,x);
	if(qr>mid) Add(p<<1|1,mid+1,r,ql,qr,x);
}

int typ;
int X[N],Y[N],D;

void dfs(int *C,int u,int f,int d) {
	cmax(D,d),C[d]++;
	for(int i=head[u];i;i=e[i].nxt) {
		int v=e[i].to;
		if(v==f || vis[v])  continue;
		dfs(C,v,u,d+1);
	}
}

void Mark(int i,int k) {
	int l=A[i==1?m:i-1],r=A[i==m?1:i+1];
	vis[l]=vis[r]=k;
}

void Get(int p,int l,int r) { 
	if(QL[p]) {
        // 计算区间QL,QR到l,r的贡献
		if(typ==0) {
			int qr=QR[p];
			rep(x,QL[p],QR[p]) Mark(x,1),dfs(X,A[x],0,qr-x),Mark(x,0);
			int T=D; D=0;
			rep(x,l,r) Mark(x,1),dfs(Y,A[x],0,x-l),Mark(x,0);
			int R=Init(T+D+1);
			NTT(R,X,1),NTT(R,Y,1);
			rep(i,0,R-1) X[i]=1ll*X[i]*Y[i]%P;
			NTT(R,X,-1);
			rep(i,0,T+D) Ans[i+l-qr]+=X[i],Mod1(Ans[i+l-qr]);
			rep(i,0,R) X[i]=Y[i]=0;
		} else {
			int ql=QL[p];
			rep(x,QL[p],QR[p]) Mark(x,1),dfs(X,A[x],0,x-ql),Mark(x,0); 
			int T=D; D=0;
			rep(x,l,r) Mark(x,1),dfs(Y,A[x],0,r-x),Mark(x,0);
			int R=Init(T+D+1);
			NTT(R,X,1),NTT(R,Y,1);
			rep(i,0,R-1) X[i]=1ll*X[i]*Y[i]%P;
			NTT(R,X,-1);
			int d=ql+m-r;
			rep(i,0,T+D) Ans[i+d]+=X[i],Mod1(Ans[i+d]);
			rep(i,0,R) X[i]=Y[i]=0;
		}
		QL[p]=QR[p]=0;
	}
	if(l==r) return;
	int mid=(l+r)>>1;
	Get(p<<1,l,mid),Get(p<<1|1,mid+1,r);
}

int main() {
	freopen("path.in","r",stdin),freopen("path.out","w",stdout);
	n=rd(),k=rd();
	rep(i,1,n) Pow[i]=qpow(i,k);
	rep(i,1,n) {
		int u=rd(),v=rd();
		AddEdge(u,v),AddEdge(v,u);
	}
	if(n<=1000) return pt1::Solve(),0;
	Init(),L=1;
    // 拓扑求环
	rep(i,1,n) if(deg[i]==1) sz[Q[++R]=i]=1;
	while(L<=R) {
		int u=Q[L++]; vis[u]=1;
		for(int i=head[u];i;i=e[i].nxt) {
			int v=e[i].to;
			if(deg[v]<=1) sz[u]+=sz[v];
			if(--deg[v]==1) Q[++R]=v;
		}
	}
	for(int u=1;u<=n;++u) if(!vis[u]) {
		while(1) {
			vis[u]=1,A[++m]=u;
			int nxt=-1;
			for(int i=head[u];i;i=e[i].nxt) {
				int v=e[i].to;
				if(!vis[v]) nxt=v;
			}
			if(nxt==-1) break;
			u=nxt;
		}
		break;
	}
	if(m==1) return pt2::Solve(),0;
	fprintf(stderr,"Circle Length =%d\n",m);
	rep(i,1,n) vis[i]=0;

	k=m/2;
	rep(i,1,m) {
		Mark(i,1);
		pt2::Divide(sz[A[i]],A[i]);
		Mark(i,0);
	}
	rep(i,1,n) Ans[i]=1ll*Ans[i]*(P+1)/2%P;
	rep(i,1,n) vis[i]=0;
	rep(i,1,m-1) Add(1,1,m,i+1,min(i+k,m),i);
	typ=0,Get(1,1,m);
	rep(i,1,m-k-1) Add(1,1,m,i+k+1,m,i);
	typ=1,Get(1,1,m);
	int ans=0;
	rep(i,1,n) ans=(ans+1ll*Ans[i]*Pow[i])%P;
	ans=ans*qpow(1ll*n*(n-1)/2%P)%P;
	printf("%d\n",ans);
}

「CodePlus 2017 11 月赛」Yazid 的新生舞会

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「CodePlus 2017 11 月赛」Yazid 的新生舞会

最基本的分析这里只保留： $cnt>\frac{len}{2}\Rightarrow 2cnt>len$

对于每一个合法的区间，合法的众数显然只有一个

考虑对于每一个众数计算答案，把 $x$ 出现的位置拿出来成一个序列 $A_i$

如果选择的区间恰好包含 $A_i,A_{i+1},\cdots ,A_j$ ，那么合法的情况就是 $2(j-i+1)>R-L+1,L\in[A_{i-1}+1,A_i],R\in[A_{j},A_{j},A_{j+1}-1]$

参数分离得到 $2j-R>2i-L$

如果对于每一个 $L$ 更新答案，那么更新的是一段区间，不妨设其为 $UL,UR$

对于每个 $R$ 查询则是一段前缀和的区间

我们知道树状数组维护区间修改区间查询需要做一次差分，而这次是区间前缀和

也就是说是再高一维。。

不妨在 $UL$ 上加, $UR$ 上减，那么在 $p$ 处的更新对于在 $k$ 处的查询的贡献是

$\cfrac{(k-p+1)(k-p+2)}{2}=\cfrac{k^2+p^2-2pk+3k-3p+2}{2}$

那么直接处理这个式子即可，需要维护 $updval,p\cdot updval,p^2\cdot updval$

查询时加入 $k$ 的贡献

#include<bits/stdc++.h>
using namespace std;
using ll=long long;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
int rd(){
	int s=0; char c;
	while(c=getchar(),c<48);
	do s=s*10+c-'0';
	while(c=getchar(),c>47);
	return s;
}
enum{N=1000010};
int n;
vector <int> A[N];
ll ans,s1[N],s2[N],s3[N];
void Add(int p,ll x) {
	p--; ll a=x,b=x*p,c=x*p*p;
	for(p+=n+1;p<=n*2;p+=p&-p) s1[p]+=a,s2[p]+=b,s3[p]+=c;
}
void Add(int l,int r,int x) { Add(l,x),Add(r+1,-x); }
ll Que(int p) {
	ll r1=0,r2=0,r3=0,tp=p;
	for(p+=n;p;p-=p&-p) r1+=s1[p],r2+=s2[p],r3+=s3[p];
	return ((tp*tp+tp)*r1-(2*tp+1)*r2+r3)/2;
}
ll Que(int l,int r) { return Que(r)-Que(l-1); }
int main(){ 
	freopen("party.in","r",stdin),n=rd(),rd();
	rep(i,1,n) A[rd()].push_back(i);
	rep(k,0,n-1) if(A[k].size()) {
		rep(i,0,A[k].size()-1) {
			int p=A[k][i],l=i?A[k][i-1]:0,r=i<(int)A[k].size()-1?A[k][i+1]:n+1;
			Add(2*i-p+1,2*i-l,1);
			ans+=Que(2*(i+1)-r,2*(i+1)-p-1);
		}
		rep(i,0,A[k].size()-1) {
			int p=A[k][i],l=i?A[k][i-1]:0;
			Add(2*i-p+1,2*i-l,-1);
		}
	}
	fprintf(fopen("party.out","w"),"%lld\n",ans);
}

Luogu P7445「EZEC-7」线段树

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

Luogu P7445「EZEC-7」线段树

显然一个点是否被 $\text{push_down}$ 仅取决于所有完全包含它的操作区间权值之和

那么可以考虑对于每个节点计算概率，然后累加

反向计算一个节点不被 $\text{push_down}$ 的概率，即权值之和为 $0$ 的概率

而每个节点有自己被覆盖的概率，即 $p_i=\cfrac{l\cdot (n-r+1)}{n(n+1)/2}$

而覆盖的次数 $c$ 决定了这个概率贡献的权值，即 $p_i^c(1-p_i)^{m-c}$

由此得到一个思路：

先通过计算得到 $k$ 次覆盖权值为0的函数 $A(x)$

容易发现这样得到每个点的概率为： $A(\cfrac{p_i}{1-p_i})(1-p_i)^m$

Naive地带入多点求值，就能暴力得到

计算 $k$ 次被覆盖权值和为0的方案数

容易发现就是

$\displaystyle [x^0](\sum_{i=-1}^V x^i)^k=[x^0](x^{-1}\frac{1-x^{V+2}}{1-x})^k$

暴力展开这个式子会需要对于 $(x^{V+2}-1)^k$ 有用的项只有 $\frac{k}{V+2}$ 项

即原式 $\displaystyle =[x^k](\frac{1-x^{V+2}}{1-x})^k=\sum_{i=0}^{k}\binom{2k-i-1}{k-1}[x^i](1-x^{V+2})^k$

$\displaystyle =\sum_{i=0}^{\frac{k}{V+2}} \binom{2k-i(V+2)-1)}{k-1} (-1)^i \binom{k}{i}$

(第一个组合数是组合意义插板，第二个是二项展开)

求 $k$ 的权值需要 $O(\frac{k}{V})$ ，并不好直接卷积优化

由于涉及了类似 $[x^n]G^k(x)$ 的形式，考虑用 "另类拉格朗日反演" 求解

如果想参考一下，但是EI的课件我是真的看不懂

处理一下 $x$ 的负指数，设 $\displaystyle F(x)=\sum _{i=0}^{V+1}x^i$ ，转化为 $[x^0](\frac{F(x)}{x})^k$

然而不管是 $F(x)$ 还是 $\frac{F(x)}{x}$ 都不存在复合逆，但是 $\frac{x}{F(x)}$ 有

设 $G(x)$ 为 $\frac{x}{F(x)}$ 的复合逆

$\displaystyle [x^0](\frac{F(x)}{x})^k=[x^0](\frac{x}{F(x)})^{-k}=[x^{k}]\frac{xG'(x)}{G(x)}$

求解 $G(x)$ 即可得到所有 $k$ 的值

$G(x)$ 为 $\cfrac{x}{F(x)}=\cfrac{x (x-1)}{x^{V+2}-1}$ 的复合逆

即满足 $\cfrac{G(G-1)}{G^{V+2}-1}=x$

$xG^{V+2}-G^2+G-x=0$

这个形式还是比较易于进行牛顿迭代的

$f(z)=xz^{V+2}-z^2+z-x$

$f'(z)=(V+2)xz^{V+1}-2z+1$

$\displaystyle A=B-\frac{f(B)}{f'(B)}=B-\frac{xB^{V+2}-B^2+B-x}{(V+2)xB^{V+1}-2B+1}$

边界条件为 $[x^0]G(x)=0$ ， $[x^1]G(x)=1$

~~结果牛顿迭代需要多项式快速幂~~

有关多点求值的优化

由于我们并不需要知道每个点被操作的概率，只需要一个求和，因此可以对于每一项求出

设 $a_i=\cfrac{p_i}{1-p_i},b_i=(1-p_i)^m$ ，容易发现实际上求出的式子是

$A(\cfrac{p_i}{1-p_i})(1-p_i)^m=\sum A_j\sum_i a_i^j b_i$

对于每个 $j$ 求解，就是

$\displaystyle [x^j]\sum _i \frac{b_i}{1-a_ix}$

可以分治 $\text{NTT}$ 通分，也就是写成下式

$\displaystyle \frac{1}{\prod (1-a_ix)} \sum b_i \prod_{i!=j} (1-a_jx)$

右边就是一个经典的分治 $\text{NTT}$ 问题，再加上一次求逆即可

~~好像也不一定快吧~~

接下来就是套板板时间

const int N=1<<19,P=998244353;

int n,m,k;

ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

typedef vector <int> V;
int rev[N],w[N];
int Inv[N+1],I[N+1],J[N+1];
void Init_w() {
	int N=1; while(N<=max(n,m+1)*2+4) N<<=1;
	int t=qpow(3,(P-1)/N);
	w[N/2]=1;
	rep(i,N/2+1,N-1) w[i]=1ll*w[i-1]*t%P;
	drep(i,N/2-1,1) w[i]=w[i<<1];
	Inv[0]=Inv[1]=1;
	rep(i,2,N) Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
	rep(i,*I=*J=1,N) {
		I[i]=1ll*I[i-1]*Inv[i]%P;
		J[i]=1ll*J[i-1]*i%P;
	}
}
int Init(int n){
	int R=1,c=-1;
	while(R<=n) R<<=1,c++;
	rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
	return R;
}

/*
被隐藏的部分:!!
NTT
operator *
operator +
operator -
*/

V operator ~ (V a) {
	int n=a.size();
	if(n==1) return V{(int)qpow(a[0],P-2)};
	V b=a; b.resize((n+1)/2); b=~b;
	int R=Init(n*2);
	NTT(R,a,1),NTT(R,b,1);
	rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
	NTT(R,a,-1),a.resize(n);
	return a;
}
void Exp_Solve(V &A,V &B,int l,int r){
	static int X[N],Y[N];
	if(l==r) {
		B[l]=1ll*B[l]*Inv[l]%P;
		return;
	}
	int mid=(l+r)>>1;
	Exp_Solve(A,B,l,mid);
	int R=Init(r-l+2);
	rep(i,0,R) X[i]=Y[i]=0;
	rep(i,l,mid) X[i-l]=B[i];
	rep(i,0,r-l-1) Y[i+1]=A[i];
	NTT(R,X,1),NTT(R,Y,1);
	rep(i,0,R-1) X[i]=1ll*X[i]*Y[i]%P;
	NTT(R,X,-1);
	rep(i,mid+1,r) B[i]+=X[i-l],Mod1(B[i]);
	Exp_Solve(A,B,mid+1,r);
}
V Deri(V a){
	rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
	a.pop_back();
	return a;
}
V Integ(V a) {
	a.pb(0);
	drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Inv[i]%P;
	return a[0]=0,a;
}

V operator << (V A,const int &x){
	A.resize(A.size()+x);
	drep(i,A.size()-1,x) A[i]=A[i-x];
	rep(i,0,x-1) A[i]=0;
	return A;
}
V operator >> (V A,const int &x){
	rep(i,x,A.size()-1) A[i-x]=A[i];
	A.resize(A.size()-x);
	return A;
}

V Ln(V a){
	int n=a.size();
	a=Deri(a)*~a,a.resize(n-1);
	return Integ(a);
}
V Exp(V F){
	int n=F.size(); F=Deri(F);
	V A(n); A[0]=1;
	Exp_Solve(F,A,0,n-1);
	return A;
}
V Pow(V x,int k) {
	int d=0,n=x.size();
	while(d<n && !x[d]) d++;
	if(1ll*d*k>=n){
		rep(i,0,x.size()-1) x[i]=0;
		return x;
	}
	x=x>>d,x.resize(n),x=Ln(x);
	rep(i,0,n-1) x[i]=1ll*x[i]*k%P;
	x=Exp(x)<<(d*k),x.resize(n);
	return x;
}

V Evaluate(V F,V X){
	static int ls[N<<1],rs[N<<1],cnt;
	static V T[N<<1];
	static auto TMul=[&] (V F,V G){
		int n=F.size(),m=G.size();
		reverse(G.begin(),G.end());
		int R=Init(n);
		NTT(R,F,1),NTT(R,G,1);
		rep(i,0,R-1) F[i]=1ll*F[i]*G[i]%P;
		NTT(R,F,-1); V T(n-m+1);
		rep(i,0,n-m) T[i]=F[i+m-1];
		return T;
	};
	static function <int(int,int)> Build=[&](int l,int r) {
		int u=++cnt; ls[u]=rs[u]=0;
		if(l==r) {
			T[u]=V{1,P-X[l]};
			return u;
		}
		int mid=(l+r)>>1;
		ls[u]=Build(l,mid),rs[u]=Build(mid+1,r);
		T[u]=T[ls[u]]*T[rs[u]];
		return u;
	};
	int n=F.size(),m=X.size();
	cmax(n,m),F.resize(n),X.resize(n);
	cnt=0,Build(0,n-1);
	F.resize(n*2+1),T[1]=TMul(F,~T[1]);
	int p=0;
	rep(i,1,cnt) if(ls[i]) {
		swap(T[ls[i]],T[rs[i]]);

		int R=Init(T[i].size()),n=T[i].size(),m1=T[ls[i]].size(),m2=T[rs[i]].size();
		NTT(R,T[i],1);
		reverse(T[ls[i]].begin(),T[ls[i]].end()); reverse(T[rs[i]].begin(),T[rs[i]].end());
		NTT(R,T[ls[i]],1); NTT(R,T[rs[i]],1);
		rep(j,0,R-1) {
			T[ls[i]][j]=1ll*T[ls[i]][j]*T[i][j]%P;
			T[rs[i]][j]=1ll*T[rs[i]][j]*T[i][j]%P;
		}
		NTT(R,T[ls[i]],-1); NTT(R,T[rs[i]],-1);
		rep(j,0,n-m1) T[ls[i]][j]=T[ls[i]][j+m1-1];
		T[ls[i]].resize(n-m1+1);
		rep(j,0,n-m2) T[rs[i]][j]=T[rs[i]][j+m2-1];
		T[rs[i]].resize(n-m2+1);
			
	} else X[p++]=T[i][0];
	X.resize(m);
	return X;
}

V operator * (V A,const int &x){
	rep(i,0,A.size()-1) A[i]=1ll*A[i]*x%P;
	return A;
}

V Newton(int n){
	if(n==1) return V{0,1}; 
	V G=Newton((n+1)/2); G.resize(n);
	V T=Pow(G,k+1);
	V A=((G*T)<<1)-G*G+G-V{0,1},B=(T<<1)*(k+2)-G*2+V{1};
	A.resize(n+1),B.resize(n+1),A=A*~B,A.resize(n+1);
	return G-A;
}

V X,Y;
void Build(int l,int r){ 
	if(l==r) return;
	int prob=1ll*l*(n-r+1)%P*Inv[n]%P*Inv[n+1]%P*2%P;
	Y.pb(P+1-prob);
	X.pb(prob*qpow(P+1-prob)%P);
	int mid=(l+r)>>1;
	Build(l,mid),Build(mid+1,r);
}

int main(){
	n=rd(),m=rd(),k=rd(),Init_w();
	V F=Newton(m+1); 
	F=Deri(F)*~(F>>1),F.resize(m+1);
	int t=1,inv=qpow(k+2);
	rep(i,0,m) {
		F[i]=1ll*F[i]*t%P*J[m]%P*I[i]%P*I[m-i]%P;
		t=1ll*t*inv%P;
	}
	Build(1,n);
	X=Evaluate(F,X);
	int ans=0;
	rep(i,0,X.size()-1) ans=(ans+P+1-X[i]*qpow(Y[i],m))%P;
	Mod2(ans),printf("%d\n",ans);
}

$\displaystyle F(x,z)=\frac{1}{\displaystyle 1-z\sum_{i=-1}^{V} x^i}$

我们希望知道 $[x^0]F(x,z)$ ，然后根据 $[z^k]$ 就能得到 $k$ 次操作权值和为 $0$ 的方案数

考虑拉格朗日反演解二元函数

设 $\displaystyle G(z)=z \sum_{i=-1}^V x^i$ ，转化为求 $\displaystyle [z^1]\frac{z}{1-G(z)}$

设 $H(z)$ 为 $G(z)$ 的复合逆，带入扩展拉格朗日反演

$\displaystyle [z^1]\frac{z}{1-G(z)}=[z^0]\frac{1}{(1-z)^2}\frac{z}{H(z)}$

$H(z)$ 满足 $=z$

「FJWC2020Day5-zzq」lg

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「FJWC2020Day5-zzq」lg

设模数为 $P$

考虑对于每一个 $\gcd$ 计算 $\text{lcm}$ 之积 $F(m)$

那么可以想到强制每个数都是 $\gcd$ 的倍数，问题转化为求 $\lfloor \frac{m}{gcd}\rfloor$ 以内所有 $\text{lcm}$ 的积 $G(m)$

那么对于每个质因数依次考虑,则得到一个简单的式子

\[G(m)=\begin{aligned}\prod p_i^{\sum_{j=1}m^n-(m-\lfloor \frac{m}{p_i^j}\rfloor )^n} \end{aligned}\]

其中枚举的 $j$ 是 $p_i$ 至少出现 $j$ 次的方案数，枚举的 $j$ 是 $\log m$ 级别的

肯定是先求出指数 $\mod \varphi(P)$ ，可以线性预处理出所有的 $i^n \mod \varphi(P)$

对于每个 $p_i$ 求出指数后还要快速幂，复杂度就是 $O(|p_i|\log P)=O(m)$ ，实际带有一些常数

那么求 $G(m)$ 的复杂度上界是 $O(m\log m)$ ，实际上 $\lfloor \frac{m}{i}\rfloor$ 有很多重复，复杂度要低很多

得到的每个 $\gcd$ 的答案还要把强制取出的 $\gcd$ 补上，是 $^

「GDOI2020模拟赛day2」我的朋友们

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

「GDOI2020模拟赛day2」我的朋友们

默认翻转了 $a_i$ 顺序，下文多项式省略 $(x)$

设当前区间为 $(i-L,i]$ 到最终结束的期望步数为 $dp_i$

令 $A_i=a_ix+1-a_i$

令 $P_i=\prod_{j\in(i-L,i]}A_j$

令 $F_i=\sum_{j=1}^i dp_jx^j$

令 $\displaystyle G_i=F_{i-1}P_i,dp_i=\frac{[x^i]G_i}{coef_i}$

其中常数 $coef_i=1-[x^0]P_i=1-\prod_{j\in(i-L,i]}(1-a_j)$ ，容易预处理得到

用分治 $\text{NTT}$ ，过程中维护

$P_{l,r}=\prod_{j\in(r-L,l]}A_j$

$\displaystyle G_{l,r}=F_{l-1}P_{l,r}$

显然 $G_{i,i}=G_i$

分治转移如下：

$\displaystyle G_{l,mid}=G_{l,r}\prod_{j\in (mid-L,min(l,r-L)]} A_j$

$P_{l,mid}=P_{l,r}\prod_{j\in (mid-L,min(l,r-L)]} A_j$

先求出 $G_{l,mid}$

$P_{mid+1,r}=P_{l,r}\prod_{j\in (max(r-L,l),mid+1]} A_j$

$\displaystyle G_{mid+1,r}=(G_{l,r}+(F_{mid}-F_{l-1})P_{l,r})\prod_{j\in (max(r-L,l),mid+1]} A_j$

归纳上述过程，发现实际上同步维护的部分只需要维护

$G_{l,r}$ 的 $[l-(r-l),r]$ 项， $P_{l,r}$ 的 $[0,r-l]$ 项

只需要实现

通过 $G_{l,r}$ 的 $[l-(r-l),r]$ 项得到 $G_{l,mid}$ 的 $[l-(mid-l),mid]$ ，以及 $G_{mid+1,r}$ 的 $[l,r]$

通过 $P_{l,r}$ 的 $[0,r-l]$ 得到 $P_{l,mid}$ 的 $[0,mid-l]$ ， $P_{mid+1,r}$ 的 $[0,r-mid-1]$

初始状态

$P_{L,n}=\prod_{j\in(n-L,L]}A_j$

$G_{L,n}=0$

在分治过程中有非常多的 $P_{l,r}$ 都是空的。。

可以分治 $\text{NTT}$ 预处理出转移系数 $\displaystyle \prod_{i=l}^r A_i$ ，或者用记忆化暴力求出，复杂度为 $O(n\log ^2n)$

在转移过程中维护的部分长度与 $r-l$ 同阶，因此复杂度为 $O(n\log^2 n)$

orangejuice's blog