做题总结

CF671E Organizing a Race

题意

有一条无限长的公路，公路上依次有 \(n​\) 个城市。每经过一个城市可以加 \(g_i​\) 升的油，一升油可以行使 \(1​\) 单位距离，两个相邻城市之间的距离是 \(w_i​\)。出发时油量为所在城市的\(g_i​\)。

你需要在 \(l,r\) 之间往返，即需要从 \(l\) 出发向右开到 \(r\)，并且从 \(r\) 出发向左开到 \(l\)。

你可以让任意个城市的加油量 \(+x(x>0)\)，要求 \(\sum x\le k\)。求最多能在多少个城市之间往返，即最大的 \(r-l+1\)。

数据范围：\(n\le 10^5,w_i,g_i,k\le 10^9\)

solution

记 \(pre_i\) 表示从 \(1\) 出发走到 \(i\)，至少要再加多少油。\(suf_i\) 表示从 \(i\) 出发走到 \(1\)，至少要再加多少油。

\[\begin{aligned} pre_i=pre_{i-1}-g_{i-1}+w_{i-1}\\ suf_i=suf_{i-1}-g_i+w_{i-1} \end{aligned} \]

这样我们做差分就可以分别求出 \(l\to r\) 至少要再加多少油和 \(r\to l\) 至少要再加多少油。

现在我们要做的就是尽量使这两部分共用的尽量多。先满足\(l\to r\) ，所需的加油量至少是 \(cost(l,r)=\max_{l<i\le r}\{pre_i\}-pre_l\)，贪心地，让 \(l\to r\) 加油的地方一定要越靠右越好。也就是，如果存在一个 \(j\) 满足 \(pre_j>pre_i\)，那么最优的方案一定是在 \(g_{j-1}\) 的地方加。这个关系可以使用单调栈维护，并且会形成一个森林。

我们在森林上 DFS，二分出最靠右的一个 \(r\)，满足 \(\max_{i<j\le r}\{pre_j\}-pre_i\le k\) ，这个 \(r\) 一定是某个祖先的左边的一个。

现在我们需要只满足 \(r\to l\) ，所以把剩下的部分全部加在 \(g_r\) 上即可。这一部分的代价是 \(suf'_r-\min_{l\le i\le r}\{suf'_i\}\)。这里的 \(suf'\) 指的是在上面的修改 \(g\) 的影响下，重新计算的 \(suf\)。

也就是，在满足了上面二分出的 \(r\) 的条件下，只需要满足 \(cost(l,r)+suf'_r-\min_{l\le i\le r}\{suf'_i\}\le k\) 即可。注意到 \(cost(l,r)\) 造成的修改与 \(suf'_{r}\) 的关系。\(cost(l,r)\) 的修改对 \(suf\) 造成了后缀减，也就是 \(cost(l,r)\) 每 \(+1\)，\(suf_r\) 都会 \(-1\)。

\[cost(l,r)+suf'_r=suf'r \]

需要满足的条件变成

\[suf_r-\min_{l\le i< r}\{suf'_i\}\le k \]

考虑线段树上二分，\(i\le l\) 的限制可以让 \([1,l-1]\) 的部分的 \(suf'\) 变成 \(+\infin\)，\(i< r\) 的限制可以让 \([r,n]\) 的部分的 \(suf'\) 变成 \(-\infin\)

现在我们需要的就是维护 \(a_i-\min_{1\le j<i} b_j\) 这个东西，需要支持区间加。

线段树的一个节点维护 \(3\) 个东西

• \(\min suf\)，记为 \(mn1\)
• \(\min suf'\)，记为 \(mn2\)
• \(\min_{mid<x\le r}\{suf_x-\min_{l\le i< x}\{suf'_i\}\}\)，即考虑整个区间的 \(suf'\)，\(r\) 只在右子树里取的最小代价。记为 \(val\)

\(update\) 的时候需要 \(O(\log n)\) 在子树里二分。假设当前正在考虑的区间的左侧的 \(\min suf'\) 是 \(mn\)，那么

• 如果 \(mn2_{ls}<mn\)，那么当前区间左侧的 \(mn2\) 不需要考虑了，答案在 \(val_{id}\) 和 左子树里取 \(\min\)
• 否则，左区间的 \(\min suf'\) 被 \(mn\) 覆盖，答案在 \(mn1_{ls}-mn\) 和右子树里取 \(\min\)

复杂度 \(O(n\log ^2n)\)

线段树上二分的时候，假设当前正在考虑的区间的左侧的 \(\min suf'\) 是 \(mn\)，那么

• 如果 \(mn<mn2_{ls}\)，\(ls\) 内的条件变为 \(suf_r-mn\le k\)，直接二分即可，然后继续考虑右子树（注意每一个线段树上的节点最多会向左二分一次，复杂度也是 \(O(n\log^2n)\)）
• 否则与 \(mn\) 无关，根据 \(val_{id}\) 与 \(k\) 的大小关系判断是左子树还是右子树

总时间复杂度 \(O(n\log^2n)\)

view code

#include <bits/stdc++.h>
using namespace std;
const int N=1e5+5,inf=1e9;
#define ll long long
const ll Inf=0x3f3f3f3f3f3f3f3f;
int n,g[N],w[N],fa[N];
ll pre[N],suf[N],k;
vector<int> e[N];
#define pb push_back
int stac[N],top;
namespace SGT{
	#define lid (id<<1)
	#define rid (id<<1|1)
	#define mid ((l+r)>>1)
	ll mn1[N<<2],mn2[N<<2],val[N<<2],tag[N<<2];
	inline void add(int id,ll a){
		mn2[id]+=a;val[id]-=a;tag[id]+=a;
	}
	inline void pushdown(int id){
		add(lid,tag[id]);
		add(rid,tag[id]);
		tag[id]=0;
	}
	ll calc(int id,int l,int r,ll mn){
		if(l==r)return mn1[id]-mn;
		pushdown(id);
		if(mn>=mn2[lid])return min(val[id],calc(lid,l,mid,mn));
		else return min(mn1[lid]-mn,calc(rid,mid+1,r,mn));
	}
	inline void update(int id,int l,int r){
		mn1[id]=min(mn1[lid],mn1[rid]);
		mn2[id]=min(mn2[lid],mn2[rid]);
		val[id]=calc(rid,mid+1,r,mn2[lid]);
	}
	void build(int id,int l,int r){
		if(l==r){
			mn1[id]=mn2[id]=suf[l];
			val[id]=0;
			return;
		}
		build(lid,l,mid);build(rid,mid+1,r);
		update(id,l,r);
	}
	void modify(int id,int l,int r,int L,int R,ll val){
		if(L<=l&&r<=R){add(id,val);return;}
		pushdown(id);
		if(L<=mid)modify(lid,l,mid,L,R,val);
		if(R>mid)modify(rid,mid+1,r,L,R,val);
		update(id,l,r);
	}
	int query1(int id,int l,int r,ll lim){
		if(l==r)return mn1[id]<=lim?l:0;
		pushdown(id);
		if(mn1[rid]<=lim)return query1(rid,mid+1,r,lim);
		else return query1(lid,l,mid,lim);
	}
	int query(int id,int l,int r,ll mn){
		if(l==r)return mn1[id]-mn<=k?l:0;
		pushdown(id);
		if(mn<=mn2[lid])
			return max(query1(lid,l,mid,k+mn),query(rid,mid+1,r,mn));
		else{
			if(val[id]<=k)return query(rid,mid+1,r,mn2[lid]);
			else return query(lid,l,mid,mn);
		}
	}
}
int ans=0;
void dfs(int u){
	stac[++top]=u;
	if(u!=n+1&&fa[u]<=n)
		SGT::modify(1,1,n,fa[u]-1,n,pre[u]-pre[fa[u]]);
	if(u!=n+1){
		int l=1,r=top-1,ret=1;
		while(l<=r){
			if(pre[stac[mid]]-pre[u]>k)ret=mid,l=mid+1;
			else r=mid-1;
		}
		ret=stac[ret]-1;
		if(u>1)SGT::modify(1,1,n,1,u-1,Inf);
		SGT::modify(1,1,n,ret,n,-Inf);
		ans=max(ans,SGT::query(1,1,n,Inf)-u+1);
		if(u>1)SGT::modify(1,1,n,1,u-1,-Inf);
		SGT::modify(1,1,n,ret,n,Inf);
	}
	for(int v:e[u])dfs(v);
	if(u!=n+1&&fa[u]<=n)
		SGT::modify(1,1,n,fa[u]-1,n,-pre[u]+pre[fa[u]]);
	--top;
}
int main(){
	n=read();k=read();
	for(int i=1;i<n;++i)w[i]=read();
	for(int i=1;i<=n;++i)g[i]=read();
	for(int i=2;i<=n;++i){
		pre[i]=pre[i-1]+w[i-1]-g[i-1];
		suf[i]=suf[i-1]+w[i-1]-g[i];
		
	}
	pre[n+1]=suf[n+1]=Inf;stac[top=1]=n+1;
	for(int i=n;i>=1;--i){
		while(pre[stac[top]]<=pre[i])--top;
		fa[i]=stac[top];e[stac[top]].pb(i);
		stac[++top]=i;
	}
	top=0;
	SGT::build(1,1,n);
	dfs(n+1);
	printf("%d\n",ans);
	return 0;
}

P5441 【XR-2】伤痕

题意

给一个 \(n\) 个点的混合完全图的边定向，\(n\) 为奇数，最多有 \(n\) 条双向边。选择四个点作为四元组，要求这四个点的导出子图强联通。最大化这样的四元组个数，并构造方案。

数据范围：\(1\le n\le 99\)（其实可以出 \(5000\) 的，这个范围只是方便 spj）

solution

四个点不强连通一定是以下三种情况之一

• 一个点向剩下 \(3\) 个点连有向边，自己没有入边
• 剩下 \(3\) 个点向一个点连有向边，自己没有出边
• 有两条没有公共点的双向边，其中一条的两个端点分别向另一条的两个端点连有向边

首先考虑计算答案。双向边一定是越多越好，剩下 \(\frac{n(n-3)}{2}\) 条有向边。设每个点的（有向边的）出度为 \(d_i\)，满足 \(\sum d_i=\frac{n(n-3)}{2}\)。那么第一种情况的个数就是

\[\sum_{i=1}^n\binom{d_i}{3} \]

注意到

\[\forall x<y,\binom{x}{3}+\binom{y}{3}\le \binom{x-1}{3}+\binom{y+1}{3} \]

所以每个点的 \(d_i\) 都相等，\(=\frac{n-3}{2}\) 时，第一种情况的数量最少。大胆猜想一定可以构造出方案使得第一种情况的个数最少，并且不存在第 \(2,3\) 种情况。方案如下：

• 对于每个点，向 \(i+1\sim i+\frac{n-3}{2}\) 连有向边，向 \(i+\frac{n-1}{2}\) 和 \(i+\frac{n+1}{2}\) 连双向边（ \(i+\frac{n-1}{2}\) 和 \(i+\frac{n+1}{2}\) 向 \(i\) 同样会向 \(i\) 连双向边，所以双向边的数量是 \(n\)）

第二种情况存在，当且仅当 \(a<b<c<d<a+\frac{n-1}{2}\)，\(a,b,c\) 都会向 \(d\) 连边，但是因为 \(a\) 也会向 \(b,c,d\) 连边，所以属于第一种。

第三种情况一定不存在。假设存在 \((a,b)\) 和 \((c,d)\) 之间存在双向边，那么 \(a<c<b<d\)，而 \(a\) 会向 \(c\) 连边，\(c\) 向 \(b\) 连边，与第三种情况不符。

综上，答案为 \(\binom{n}{4}-n\times \binom{(n-3)/2}{3}\)，构造方案为对于每个点，向 \(i+1\sim i+\frac{n-3}{2}\) 连有向边，向 \(i+\frac{n-1}{2}\) 和 \(i+\frac{n+1}{2}\) 连双向边

view code

#include <bits/stdc++.h>
using namespace std;
const int N=1005;
#define ll long long
ll C(int x,int y){
	ll ret=1;
	for(int i=0;i<y;++i)ret*=x-i;
	for(int i=1;i<=y;++i)ret/=i;
	return ret;
}
int n,g[N][N];
inline void adde(int u,int v){
	u=(u-1)%n+1;v=(v-1)%n+1;
	g[u][v]=1;
}
int main(){
	scanf("%d",&n);
	if(n==1){puts("0\n0");return 0;}
	else if(n==3){
		puts("0\n0 1 1\n0 0 1\n0 1 0");
		return 0;
	}
	printf("%lld\n",C(n,4)-n*C((n-3)/2,3));
	for(int i=1;i<=n;++i)
		for(int j=1;j<=(n+1)/2;++j)
			adde(i,i+j);
	for(int i=1;i<=n;++i,puts(""))for(int j=1;j<=n;++j)printf("%d ",g[i][j]);
	return 0;
}

P3642 [APIO2016]烟火表演

题意

给定一棵有根树，边有边权，将一条边的边权修改 \(1\) 的代价是 \(1\)。求最小的代价使得所有叶子到根的距离相等。

数据范围：\(n\le 3\times 10^5\)

solution

有一个显然的DP，记 \(f_{u,x}\) 表示以 \(u\) 为根的子树，所有叶子到 \(u\) 的距离为 \(x\) 的最小代价。转移有

\[f_{u,x}=\sum \min_{y\le x}\{f_{v,y}+|w-x+y|\} \]

记 \(f_u(x)\) 为 \(f_{u,x}\)，容易发现 \(f_u(x)\) 为连续的下凸函数，且斜率均为整数。

注意到原函数斜率均为整数，\(|w-x+y|\) 的斜率 \(\in \{-1,0,1\}\)，只需要考虑 \(y\) 在 \(f_v\) 斜率 \(=0\) 的时候特别考虑，记 \(f_v\) 斜率 \(=0\) 的横坐标为 \([L,R]\)，对 \(y\) 在 \([0,L],L,[L,R],R\) 分别考虑。

考虑 \(f_v(y)\) 对 \(f_u(x)\) 的影响：

• 对于 \(x\le L\)，那么 \(y\) 一定取最大值 \(x\) 最优，也就是 \(f_{u}(x)\gets f_v(x)+w\)。因为如果 \(y\) 每减 \(a\)，\(f_v(y)\) 至少增加 \(y\)（斜率\(\ge -1\)），\(|w-x+y|\) 最多减少 \(y\)。
• 对于 \(x\ge R+w\)，那么 \(y\) 一定取 \(R\) 最优，也就是 \(f_u(x)\gets f_v(R)+(x-w-R)\)。原因与上面类似。
• 对于 \(L\le x\le L+w\)，那么 \(y\) 一定取 \(L\) 最优，也就是 \(f_u(x)\gets f_v(L)+(L+w-x)\)。原因与上面类似。
• 对于 \(L+w\le x\le R+w\)，都存在一个 \(y\in [L,R]\) 满足 \(|w-x+y|=0\)，所以 \(f_u(x)\gets f_v(L)=f_v(x-w)\)

所以，新的函数是：\(\le L\) 的部分向上平移，\([L+w,R+w]\) 的部分由原来的 \([L,R]\) 平移得到，中间断开的部分是 \([L,L+w]\)，用一条斜率为 \(-1\) 的线连接，最后一段斜率为 \(1\) 的直线。

考虑使用堆来维护分段函数：恒坐标从大往小，每一个点表示斜率在这个点之后 \(-1\)，有多个重复的点表示斜率 \(-k\)。

那么，进行一次合并的时候，就把两个堆合并起来（使用左偏树）即可。问题是如何找到斜率 \(=0\) 的横坐标区间 \([L,R]\)，我们发现一个节点的儿子数量一定和 \(f_u\) 最右侧的斜率一致，因为每一个儿子合并上来的时候在末尾都有且仅有一段斜率为 \(1\) 的直线。点 \(u\) 所表示的堆里，弹出 \(|son_u|\) 个，即可知道 \(R\)，再弹一个，就是 \(L\)。

那么最后怎么求答案呢。根据上面的做法我们可以求出整个堆里的点所表示出的分段函数的斜率，斜率即 \(f_1(0)\)，应该是所有边的边权之和，那么还原出整个分段函数就是很简单的了。求出斜率为 \(0\) 的区间的对应 \(f_1(L)\) 即可。

时间复杂度 \(O(n\log n)\)。

view code

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=3e5+5;
namespace Heap{
	int ls[N<<1],val[N<<1],rs[N<<1],inde,dep[N<<1];
	inline int newnode(int x){
		val[++inde]=x;
		return inde;
	} 
	inline int merge(int x,int y){
		if(!x||!y)return x|y;
		if(val[x]<val[y])swap(x,y);
		rs[x]=merge(rs[x],y);
		if(dep[ls[x]]<dep[rs[x]])swap(ls[x],rs[x]);
		dep[x]=dep[rs[x]]+1;
		return x;
	}
	inline void pop(int &x){
		x=merge(ls[x],rs[x]);
	}
}
using namespace Heap;
int n,m,fa[N],c[N],deg[N],rt[N];
#define ll long long
ll sum;
signed main(){
	n=read();m=read();
	for(int i=2;i<=n+m;++i){
		fa[i]=read();c[i]=read();
		++deg[fa[i]];
		sum+=c[i];
	}
	cerr<<sum<<endl;
	for(int i=n+m;i>1;--i){
		for(;deg[i]>1;--deg[i])pop(rt[i]);
		int w=c[i];
		int R=val[rt[i]];pop(rt[i]);
		int L=val[rt[i]];pop(rt[i]);
		rt[fa[i]]=merge(merge(rt[fa[i]],rt[i]),merge(newnode(L+w),newnode(R+w)));
	}
	for(;deg[1]--;)pop(rt[1]);
	for(;rt[1];)sum-=val[rt[1]],pop(rt[1]);
	printf("%lld\n",sum);
	return 0;
}

CF1583F Defender of Childhood Dreams

题意

一个 \(n\) 个点的有向完全图，若 \(x<y\) ，则有一条 \(x\to y\) 的有向边。你需要给边染色，要求满足不存在一条长度 \(\ge k\) 的路径上的边颜色全部相同。求最小颜色数并构造方案。

数据范围：\(2\le k<n\le 1000\)

solution

颜色数为 \(\lceil\log_kn\rceil\) 一定可行，下面给出构造方案：

我们将 \(1\sim k,k+1\sim 2k,\cdots,(m-1)k+1\sim mk\) 之间的点连颜色为 \(1\) 的边，因为这些点之间的路径最长为 \(k-1\)。然后考虑跨过这些分组之间的连边，容易发现，这等价于一个 \(1\sim m\) 的子问题，被连边的组在原图上对应的点一一连颜色为 \(2\) 的边。不断做下去即可。

时间复杂度 \(O(n^2\log_k n)\)。

view code

#include <bits/stdc++.h>
using namespace std;
int n,k,cnt=0,s=1;
int main(){
	scanf("%d %d",&n,&k);
	for(;s<n;)s*=k,++cnt;
	printf("%d\n",cnt);
	for(int i=0;i<n;++i)for(int j=i+1;j<n;++j){
		cnt=0;
		for(int x=i,y=j;x!=y;x/=k,y/=k)++cnt;
		printf("%d ",cnt);
	}
	return 0;
}

CF1583G Omkar and Time Travel

题意

有 \(n\) 个事件，当你在时间 \(b_i\) 时，你需要从时间 \(b_i\) 传送到 \(a_i(a_i<b_i)\)，完成这个事件，所有 \(a_i,b_i\) 互不相同。传送之后，所有 \(a_j\) 在 \([a_i,b_i]\) 内的事件会变成未完成。给定一个事件集合 \(S\)，求完成这个集合内所有事件需要传送的次数\(\bmod10^9+7\) 。

数据范围：\(n\le 2\times 10^5\)

solution

注意到如果 \(S\) 内的两个事件满足 \(a_i<a_j\land b_i<b_j\)，那么我们只需要满足完成了 \(j\) 即可，因为完成 \(j\) 的时候一定完成了 \(i\)。现在剩下的事件就是 \(a_1<a_2<\cdots <a_n<b_n<\cdots b_2<b_1\)，称之为包含关系。最后的应该是形如：走到 \(b_1\)，然后回到 \(a_1\)，然后走到 \(b_2\)，然后回到 \(a_2,\cdots\)

记 \(f_i\) 表示从 \(b_i\) 回到 \(a_i\) 再走到 \(b_i\)，这之间需要的传送次数。这个可以 \(O(n\log n)\) 使用树状数组求出。设 \(g_i\)（ \(i\) 是一个必须完成的时间 ）表示完成所有被 \(i\) 包含的，必须完成的事件的传送次数。设被 \(i\) 包含的最大的事件区间为 \(j\)。那么 \(g_i=g_j+\sum_{a_i<a_p<a_j,b_p<b_j}f_{p}\)，同样可以树状数组维护。

在最外层加一个事件 \(n+1\)，\(a_{n+1}=1,b_{n+1}=\max_{i\in S}b_S\) 也就是要求在所有存在包含关系的必选事件之前的事件也必须经历一遍，答案就是 \(g_{n+1}\)。

view code

#include <bits/stdc++.h>
using namespace std;
inline int read(){
	int s=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9')s=(s<<3)+(s<<1)+ch-'0',ch=getchar();
	return s*f;
}
const int mod=1e9+7,N=4e5+5;
inline int add(int a,int b){return (a+b>=mod)?a+b-mod:a+b;}
inline void Add(int &a,int b){a+=(a+b>=mod)?b-mod:b;}
inline void Max(int &a,int b){if(a<b)a=b;}
inline int lowbit(int x){return x&(-x);}
namespace T1{
	int mx[N],len;
	inline void modify(int x,int v){
		for(;x;x-=lowbit(x))Max(mx[x],v);
	}
	inline int query(int x){
		int ret=0;
		for(;x<=len;x+=lowbit(x))Max(ret,mx[x]);
		return ret;
	}
}
namespace T2{
	int sum[N],len;
	inline void modify(int x,int v){
		for(;x;x-=lowbit(x))Add(sum[x],v);
	}
	inline int query(int x){
		int ret=0;
		for(;x<=len;x+=lowbit(x))Add(ret,sum[x]);
		return ret;
	}
}
int a[N],b[N],id[N],n,m,q,f[N],g[N];
vector<int> e[N];
#define pb push_back
bool tag[N];
int main(){
	n=read();m=n<<1;
	for(int i=1;i<=n;++i){
		a[i]=read();b[i]=read();
		id[a[i]]=i;id[b[i]]=i;
	}
	T1::len=T2::len=m;
	q=read();
	for(int i=1;i<=q;++i)tag[read()]=1;
	for(int i=1;i<=m;++i){
		if(tag[id[i]]&&i==b[id[i]]){
			int x=T1::query(a[id[i]]);
			x=id[x];
			e[x].pb(id[i]);
			T1::modify(a[id[i]],b[id[i]]);
		}
	}
	a[n+1]=1;b[n+1]=m;
	e[id[T1::query(1)]].pb(n+1);
	for(int i=1;i<=m;++i){
		if(i!=b[id[i]])continue;
		f[i]=add(T2::query(a[id[i]]),1);
		if(tag[id[i]]){
			Add(g[id[i]],1);
			for(int v:e[id[i]])
				g[v]=add(g[id[i]],T2::query(a[v]));
		}
		T2::modify(a[id[i]],f[i]);
	}
	printf("%d\n",g[n+1]);
	return 0;
}

CF1603D

题意

记 \(c(l,r)=\sum\limits_{i=l}^r\sum\limits_{j=i}^r[\gcd(i,j)\ge l]\)。你可以将\(1\sim n\) 的\(n\) 个数分成 \(k\) 段 \(0=x_1<x_2<x_3<\cdots<x_k<x_{k+1}=n\)，最小化

\[\sum_{i=1}^{k}c(x_i+1,x_{i+1}) \]

数据范围：\(n,k\le 10^5,T\) 组询问， \(T\le 3\times 10^5\)。

solution

一个暴力的 DP 做法是：设 \(f_{i,k}\) 表示前 \(i\) 个数分成 \(k\) 段的最小值。转移是

\[f_{i,k}=\min_{j=1}^i\{f_{j-1,k-1}+c(l,r)\} \]

边界：\(f_{i,0}=+\infin\)，\(f_{0,k}=0\)。

询问数很多，可以预处理出 DP 数组 \(O(1)\) 查询。

DP状态数是 \(O(n^2)\) 的，无法承受。观察 \(c(l,r)\) 的性质，发现 \(c(l,r)\ge r-l+1\)，也就是说 \(f_{n,k}\ge n\)。其次，\(c(l,2l-1)=(2l-1)-l+1=l\)，也就是说，若 \(n<2^k\)，一定可以把 \(n\) 分成若干段，满足 \(c(l,r)=r-l+1\)，即\(f_{n,k}=n(2^k>n)\)。剩下的需要考虑的 \(k\) 就是 \(O(\log n)\) 级别的了。

状态数没有太大的优化空间了，考虑优化 \(c(l,r)\) 的计算。大力推式子：

\[\begin{aligned} c(l,r)&=\sum_{i=l}^r\sum_{j=i}^r[\gcd(i,j)\ge l]\\ &=\sum_{k=l}^r\sum_{i=l}^r\sum_{j=i}^r[\gcd(i,j)= k]\\ &=\sum_{k=l}^r\sum_{i=\lceil\frac{l}{k}\rceil}^{\lfloor\frac{r}{k}\rfloor}\sum_{j=i}^{\lfloor\frac{r}{k}\rfloor}[\gcd(i,j)=1]\\ &=\sum_{k=l}^r\sum_{j=1}^{\lfloor\frac{r}{k}\rfloor}\sum_{i=1}^j[\gcd(i,j)=1]\\ &=\sum_{k=l}^r\sum_{i=1}^{\lfloor\frac{r}{k}\rfloor}\varphi(i) \end{aligned} \]

记 \(s(n)=\sum_{i=1}^n\varphi(i)\)，即 \(\varphi\) 的前缀和。则

\[c(l,r)=\sum_{k=l}^rs(\lfloor\frac{r}{k}\rfloor) \]

预处理出 \(s(n)\) 之后即可 \(O(\sqrt{r})\) 查询。实际上，\(c(l,r)\) 的计算复杂度大概是 \(O(\sqrt{r-l})\) 的。或者可以 \(O(n\sqrt{n})\) 预处理 \(O(1)\) 查询。

考虑优化 DP 转移。根据这个DP分段的形式，大胆猜测 \(c(l,r)\) 满足四边形不等式，那么这样DP转移就可以使用决策单调性优化转移了。\(c(l,r)\) 满足四边形不等式的证明如下：

设 \(l_1<l_2<r_1<r_2\)，我们需要证明 \(c(l_1,r_1)+c(l_2,r_2)\le c(l_1,r_2)+c(l_2,r_1)\)。

设 \(f(l,r,p)=\sum_{k=l}^p(\lfloor\frac{r}{k}\rfloor)\)，那么

\[\begin{aligned} c(l_1,r_2)+c(l_2,r_1)&=f(l_1,r_2,r_2)+f(l_2,r_1,r_1)\\ &=(f(l_1,r_2,l_2-1)+f(l_2,r_2,r_2))+(f(l_1,r_1,r_1)-f(l_1,r_1,l_2-1))\\ &=f(l_1,r_2,l_2-1)-f(l_1,r_1,l_2-1)+c(l_1,r_1)+c(l_2,r_2)\\ \end{aligned} \]

显然 \(f(l_1,r_2,l_2-1)\ge f(l_1,r_1,l_2-1)\)，得证。

剩下的部分可以使用1D1D决策单调性优化DP的常见套路（分治/单调队列）进行，可以参考我的博客。因为有 \(O(\log n)\) 层，每层的复杂度是 \(O(n\log n)\) 的，转移的时间复杂度是 \(O(n\log^2n)\) 的，加上预处理的复杂度，总时间复杂度 \(O(n\log^2n+n\sqrt{n})\) ，空间复杂度 \(O(n\sqrt{n})\)。

实际上，采用分治做法，假设当前的转移区间是 \([L,R]\)，当前需要寻找转移点的是 \(mid\)。首先\(O(\sqrt{r-l})\) 求出 \(c(R+1,mid)\)，然后从 \(\min(mid,R)\) 到 \(L\) 枚举转移点 \(i\)，\(c(i,mid)=c(i+1,mid)+s(\lfloor\frac{mid}{i}\rfloor)\)。这个做法的时间复杂度不太好证，不过预处理跑的很快，而且代码复杂度很低，下面的代码正是这个写法。

view code

#include <bits/stdc++.h>
using namespace std;
const int N=1e5+5,LOG=18;
inline int read(){
	int s=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9')s=(s<<3)+(s<<1)+ch-'0',ch=getchar();
	return s*f;
}
#define ll long long
const ll inf=1e18;
int prime[N],pcnt;
ll phi[N];
bool vis[N];
void init(int n){
	phi[1]=1;
	for(int i=2;i<=n;++i){
		if(!vis[i]){
			prime[++pcnt]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=pcnt&&prime[j]*i<=n;++j){
			vis[prime[j]*i]=1;
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
	for(int i=1;i<=n;++i)phi[i]+=phi[i-1];
}
inline ll calc(int l,int r){
	ll ret=0;
	for(int i;l<=r;l=i+1){
		i=min(r/(r/l),r);
		ret+=1ll*phi[r/l]*(i-l+1);
	}
	return ret;
}
inline ll c(int l,int r){return phi[r/l];}
int T,n,k;
ll f[N][LOG];
void solve(int l,int r,int L,int R){
	if(l>r)return;
	int mid=(l+r)>>1;ll sum=calc(R+1,mid);
	int pos=0;ll mn=inf; 
	for(int i=min(mid,R);i>=L;--i){
		sum+=phi[mid/i];
		if(sum+f[i-1][k-1]<=mn){
			mn=sum+f[i-1][k-1];
			pos=i;
		}
	}
	f[mid][k]=mn;
	solve(l,mid-1,L,pos);
	solve(mid+1,r,pos,R); 
}
int main(){
	n=100000;init(n);
	for(int i=1;i<=n;++i)f[i][1]=(1ll*i*(i+1))>>1;
	for(k=2;k<=17;++k)
		solve(1,n,1,n);
	T=read();
	while(T--){
		n=read();k=read();
		if(k>=20||n<(1<<k)){
			printf("%d\n",n);
			continue;
		}
		printf("%lld\n",f[n][k]);
	}
	return 0;
}

CF1603E A Perfect Problem

题意

求长为 \(n\)，值域在 \([1,n+1]\) 的序列，满足其所有子序列 \(a\) 满足 \(\min\{a_i\}\times \max\{a_i\}\ge \sum a_i\) 的序列数 \(\bmod p\)。

数据范围：\(n\le 200,p\) 为质数。

solution

容易发现，一个序列满足条件当且仅当其从小到大排完序之后所有前缀都满足条件。设排完序后依次是 $a_1,a_2,\cdots,a_n $。

那么原来的限制就是 \(a_1a_n\ge \sum a_i​\)，逐步发现性质

• \(a_i\ge i\)

否则 \(\sum_{j=1}^ia_j\ge a_1i>a_1a_i\)

• 若存在 \(a_i=i\)，则 \(\forall j\in [1,i],a_i=a_j\)

因为 \(a_1i=a_1a_i\ge \sum_{j=1}^ia_j\ge a_1i\)，所以 \(\forall j\in [1,i],a_i=a_j\)

• 若不存在 \(a_i=i\)，那么 \(a_n=n+1\)，则有 \(a_1(n+1)\ge \sum a_i\)，移项，得 \(a_1\ge \sum a_i-a_1\)

设 \(f_{i,j,d}\) 表示考虑填 \([i,n+1]\) 内的数，从大到小填了 \(j\) 个数，这 \(j\) 个数 \(-i\) 的差之和为 \(d\) ，不存在一个 \(a_i=i\) 的方案数。因为是序列数而不是有序序列数，所以需要乘上 \(n!\)。

转移形如

\[f_{i,j+cnt,d+j}\gets \frac{f_{i+1,j,d}}{cnt!} \]

对答案的贡献可以枚举最后一段填什么，这样有了一个状态数 \(O(n^3)\)，转移 \(O(n)\) 的 \(O(n^4)\) 的做法，已经足以通过。

• 若 \(a_1\le n-\sqrt{2n}\)，那么不存在一个合法的序列

因为 \(\ge a_1\) 的数有 \(\sqrt{2n}\) 个，这些数 \(-a_1\) 的差之和是 \(\sum_{i=1}^{\sqrt{2n}} i=n+\frac{\sqrt{2n}}{2}\)，已经 \(> a_1\)了。

所以状态数只剩下 \(O(n^2\sqrt{n})\)，总时间复杂度 \(O(n^3\sqrt{n})\) 且常数较小。

view code

#include <bits/stdc++.h>
using namespace std;
const int N=305;
int n,mod,fac[N],invfac[N];
inline int quick_pow(int a,int b){
	int ret=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ret=1ll*ret*a%mod;
	return ret;
}
inline int Inv(int a){return quick_pow(a,mod-2);}
inline void init(){
	fac[0]=1;for(int i=1;i<=n;++i)fac[i]=1ll*fac[i-1]*i%mod;
	invfac[n]=Inv(fac[n]);
	for(int i=n-1;~i;--i)invfac[i]=1ll*invfac[i+1]*(i+1)%mod;
}
int f[N][N][N],ans;
int main(){
	scanf("%d%d",&n,&mod);
	init(); 
	int lim=max(1,(int)(n-sqrt(2*n)-2));
	f[n+2][0][0]=1;
	for(int i=n+1;i>=lim;--i)for(int j=max(0,n-i);j<n;++j){
		for(int d=0;d+j<=n;++d)if(f[i+1][j][d]){
			int v=f[i+1][j][d];
			if(i>n-j){
				for(int cnt=0;cnt+j<=n&&d+j<=n;++cnt){
					f[i][j+cnt][d+j]=(f[i][j+cnt][d+j]+1ll*v*invfac[cnt])%mod;
					if(cnt&&j+cnt==n&&i>=d+j)ans=(ans+1ll*v*invfac[cnt])%mod;
				}
			}else if(i==n-j&&j+d<=i)
				ans=(ans+1ll*v*invfac[i])%mod;
		}
	}
	printf("%d\n",1ll*ans*fac[n]%mod);
	return 0;
} 

CF1606F Tree Queries&P7897 [Ynoi2006] spxmcq

题意

一棵以 \(1\) 为根的有根树，点有点权。一次询问会将所有点的点权 \(+x\)，求 \(u\) 的子树内，包含 \(u\) 的连通点集中，点权的权值之和最大是多少。

数据范围：\(n,q\le 10^6,|x|,|a_i|\le 10^8\)

solution

首先把 \(u​\) 去掉，考虑 \(u​\) 的所有儿子。有显然的DP：设 \(f_p\) 表示 \(p\) 的子树内（必须选 \(p\)）的答案，那么有：

\[f_p=a_p+x+\sum_{v\in son_p}\max(0,f_v) \]

\(\max\) 非常麻烦，如果所有的 \(f_v\) 都 \(\ge 0\) 的话，记 \(sz_p\) 表示 \(p\) 的子树大小，\(sum_u\) 表示 \(p​\) 的子树内的权值和，那么：

\[f_p=sz_p\times x+sum_p \]

也就是说，\(f_p< 0\) 的 \(p\) 不对父亲做出贡献，可以视作这条边不存在，也就是成为了若干有根树森林。在 \(f_p\ge 0\) 时 \(p\) 对父亲做出贡献，需要把这条边连上。

而 \(f_p\ge 0\)，则意味着 \(sz_p\times x+sum_p\ge 0\)，即

\[x\ge \left\lceil\frac{-sum_p}{sz_p}\right\rceil \]

也就是 \(x\ge \left\lceil\frac{-sum_p}{sz_p}\right\rceil\) 时，\(p\) 与父亲合并。将询问的 \(x\) 从小到大排序，找出需要与父亲合并的最小的阈值（使用 priority_queue 维护）。然后需要做的是将一段祖先进行链加。改成单点加，dfn 区间求和，可以使用树状数组维护。

时间复杂度 \(O(n\log n)\)。

view code

const int N=1e6+5;
struct Edge{int to,next;}e[N];
int head[N],ecnt;
inline void adde(int u,int v){e[++ecnt]=(Edge){v,head[u]};head[u]=ecnt;}
int n;
#define pb push_back
int q,dep[N],pa[N],inde,dfn[N],ed[N];
void dfs1(int u,int fat){
	dfn[u]=++inde;
	dep[u]=dep[fat]+1;
	pa[u]=fat;
	for(int i=head[u],v;i;i=e[i].next){
		v=e[i].to;
		if(v==fat)continue;
		dfs1(v,u);
	}
	ed[u]=++inde;
}
#define ll long long 
namespace dsu{
	int fa[N];
	inline int find_fa(int x){return x==fa[x]?x:fa[x]=find_fa(fa[x]);}
	inline void merge(int x,int y){
		x=find_fa(x);y=find_fa(y);
		if(x==y)return;
		if(dep[x]>dep[y])swap(x,y);
		fa[y]=x;
	}
	inline void init(){for(int i=1;i<=n;++i)fa[i]=i;}
}
struct BIT{
	ll tr[N<<1];
	inline int lowbit(int x){return x&(-x);}
	inline void add(int x,ll v){
		for(;x<=inde;x+=lowbit(x))tr[x]+=v;
	}
	inline ll query(int x){
		ll ret=0;
		for(;x;x-=lowbit(x))ret+=tr[x];
		return ret;
	}
	inline ll query(ll l,ll r){return query(r)-query(l-1);}
}T1,T2;
struct Query{int u,k,id;}qry[N];
inline bool operator <(Query a,Query b){
	return a.k==b.k?dep[a.u]>dep[b.u]:a.k<b.k;
}
inline int Ceil(ll a,ll b){return (a<0)?a/b:(a-1)/b+1;}
inline ll calc(int u){
	return Ceil(-T1.query(dfn[u],ed[u]-1),T2.query(dfn[u],ed[u]-1));
}
inline void merge(int x,int y){
	x=dsu::find_fa(x);
	ll sum1=T1.query(dfn[y],ed[y]-1),sz1=T2.query(dfn[y],ed[y]-1);
	T1.add(ed[y],sum1);T1.add(ed[x],-sum1);
	T2.add(ed[y],sz1);T2.add(ed[x],-sz1);
	dsu::merge(x,y);
}
#define pr pair<ll,int>
#define se second
#define mp make_pair
priority_queue<pr,vector<pr>,greater<pr> > que;
ll ans[N],a[N];
bool vis[N];
int main(){
	read(n);read(q);
	for(int i=2,u;i<=n;++i){read(u);adde(u,i);}
	for(int i=1;i<=n;++i)read(a[i]);
	dsu::init();dfs1(1,0);
	for(int i=1,x,k;i<=q;++i){
		read(x);read(k);
		qry[i]=(Query){x,k,i};
	}
	sort(qry+1,qry+1+q);
	for(int u=2;u<=n;++u)
		T2.add(dfn[u],1),T2.add(ed[u],-1),
		T1.add(dfn[u],a[u]),T1.add(ed[u],-a[u]);
	for(int i=2;i<=n;++i)que.push(mp(calc(i),i));
	for(int i=1;i<=q;++i){
		while(!que.empty()&&qry[i].k>=que.top().fi){
			pr cur=que.top();que.pop();
			int u=cur.se;
			if(vis[u])continue;
			vis[u]=1;
			merge(pa[u],u);
			u=dsu::find_fa(u);
			if(u>1)que.push(mp(calc(u),u));
		}
		int u=qry[i].u,k=qry[i].k;
		ans[qry[i].id]=T1.query(dfn[u]+1,ed[u]-1)+k*T2.query(dfn[u]+1,ed[u]-1)+a[u]+k;
	}
	for(int i=1;i<=q;++i)putint(ans[i]);
	flush();
	return 0;
}

【UR #2】树上GCD

题意

给定一个以 \(1\) 为根的有根树，记 \(p=LCA(u,v)\)，\(f(u,v)=\gcd(dis(u,p),dis(v,p))\)。特别的，\(\gcd(0,x)=\gcd(x,0)=x\)。

对于所有 \(i\in[1,n)\)，求 \(f(u,v)=i\) 的无序对 \((u,v),u\neq v\) 的数量。

数据范围：\(n\le 2\times 10^5\)

solution

首先把答案转化为求 \(f(u,v)\) 是 \(i\) 的倍数的点对数量，最后莫比乌斯反演回来即可。

考虑点分治，对于一个分治中心 \(u\)，求出经过它的所有点对 对答案的贡献。

因为是有根树，所以原树上 \(u\) 的儿子之间的贡献是方便计算的，可以 \(O(n\log n)\) 调和级数预处理深度是 \(i\) 的倍数的点的个数，同一子树内的简单容斥掉即可。

然后考虑 \(u\) 在当前连通块内，原树上的所有父亲的其他子树内的点，与 \(u\) 子树内的点的贡献。

对于一个 \(i\)，计算 \(\gcd\ge i\) 的点对数量。也就是要求 \(dis+len(u,x)\equiv 0\pmod i\) 。考虑根号分治，对于 \(i\ge \sqrt{dep}\) 的部分，直接枚举，每次 \(+i\)，复杂度 \(O(n\sqrt{dep})\)。剩下的部分，对 \(u\) 子树内的点预处理出 \(cnt_{i,j}\) 表示 \(len(u,x)\bmod i=j\) 的 \(x\) 的数量，预处理复杂度 \(O(n\sqrt{dep})\)。

那么一层复杂度是 \(O(n\sqrt{n})\)，根据主定理，总复杂度也是 \(O(n\sqrt{n})\)。

view code

#define ll long long
#define pb push_back
const int N=200005,B=505,inf=1e9;
inline void Max(int &a,int b){if(a<b)a=b;}
inline int MOD(int x,int mod){return x<0?x+mod:x;}
namespace Mu{
	int mu[N],prime[N],pcnt;
	bool vis[N];
	inline void init(int n){
		mu[1]=1;
		for(int i=2;i<=n;++i){
			if(!vis[i])prime[++pcnt]=i,mu[i]=-1;
			for(int j=1;j<=pcnt&&prime[j]*i<=n;++j){
				vis[i*prime[j]]=1;
				if(i%prime[j]==0)break;
				mu[i*prime[j]]=-mu[i];
			}
		}
	}
}
using Mu::mu;
vector<int> e[N];
int cnt[B][B],pa[N],tot1[N],tot2[N],sum[N];
int mxdep;
inline void build(int *tot){
	for(int i=1;i<=mxdep;++i){
		sum[i]=0;
		for(int j=i;j<=mxdep;j+=i)sum[i]+=tot[j];
	}
}
bool vis[N];
int sz[N],rt,mn;
void getsz(int u,int fa){
	sz[u]=1;
	for(int v:e[u])
		if(v!=fa&&!vis[v])getsz(v,u),sz[u]+=sz[v];
}
void dfs1(int u,int fa,int tot){
	int mx=tot-sz[u];
	for(int v:e[u]){
		if(v==fa||vis[v])continue;
		dfs1(v,u,tot);
		Max(mx,sz[v]);
	}
	if(mx<mn)mn=mx,rt=u;
}
inline int findroot(int x){mn=inf,rt=0;getsz(x,0);dfs1(x,0,sz[x]);return rt;}
ll ans[N],ans2[N];
void calc(int u,int fa,int dep,int *tot){
	++tot[dep];
	Max(mxdep,dep);
	for(int v:e[u]){
		if(vis[v]||v==fa)continue;
		calc(v,u,dep+1,tot);
	}
}
int blk;
void dfs2(int u,int fa,int dep){
	for(int j=1;j<=blk;++j)++cnt[j][dep%j];
	Max(mxdep,dep);++tot2[dep];
	for(int v:e[u]){
		if(vis[v]||v==fa)continue;
		dfs2(v,u,dep+1);
	}
}
void dfs3(int u,int fa,int dep){
	if(u!=1){++ans2[1];--ans2[dep+1];}
	for(int v:e[u])if(v!=fa)dfs3(v,u,dep+1);
}
void dfs(int u){
	vis[u]=1;mxdep=0;tot1[0]=1;
	for(int v:e[u]){
		if(vis[v]||v==pa[u])continue;
		calc(v,u,1,tot1);
	}
	build(tot1);
	blk=max(1,(int)sqrt(mxdep));
	for(int i=1;i<=mxdep;++i)ans[i]+=1ll*sum[i]*(sum[i]-1)/2;
	for(int i=1;i<=blk;++i)++cnt[i][0];
	int dep1=mxdep;
	for(int v:e[u]){
		if(vis[v]||v==pa[u])continue;
		mxdep=0;
		dfs2(v,u,1);
		build(tot2);
		for(int i=1;i<=mxdep;++i)ans[i]-=1ll*sum[i]*(sum[i]-1)/2;
		memset(tot2,0,(mxdep+1)<<2);
	}
	int top=u,dis=1;
	for(int x=pa[u];x&&!vis[x];vis[x]=1,x=pa[x],++dis){
		top=x;
		mxdep=0;
		calc(x,pa[x],0,tot2);
		build(tot2);
		for(int i=1;i<=mxdep;++i){
			if(i<=blk){
				ans[i]+=1ll*sum[i]*cnt[i][MOD((-dis)%i,i)];
			}else{
				int cur=0;
				for(int j=MOD((-dis)%i,i);j<=dep1;j+=i)
					cur+=tot1[j];
				ans[i]+=1ll*sum[i]*cur;
			}
		}
		memset(tot2,0,(mxdep+1)<<2);
		vis[x]=1;
	}
	memset(tot1,0,(dep1+1)<<2);
	for(int x=u;;x=pa[x]){vis[x]=0;if(x==top)break;}
	for(int i=1;i<=blk;++i)memset(cnt[i],0,i<<2);
	vis[u]=1;
	for(int v:e[u])if(!vis[v])dfs(findroot(v));
}
int n;
int main(){
	n=read();
	Mu::init(n);
	for(int i=2;i<=n;++i){
		pa[i]=read();
		e[pa[i]].pb(i),e[i].pb(pa[i]);
	}
	dfs(findroot(1));
	dfs3(1,0,0);
	for(int i=1;i<=n;++i)ans2[i]+=ans2[i-1];
	for(int i=1;i<n;++i){
		ll cur=0;
		for(int j=1;j<=n&&j*i<=n;++j)cur+=ans[i*j]*mu[j];
		printf("%lld\n",cur+ans2[i]);
	}
	return 0;
}

loj2743. 「JOI Open 2016」摩天大楼

题意

一个互不相同的整数序列 \(a\)，求满足 \(\sum_{i=1}^{n-1}|a_{p_i}-a_{p_{i+1}}|\le L\) 的排列 \(p\) 的数量 \(\bmod 10^9+7\)。

数据范围：\(n\le 100,a_i,L\le 1000\)

solution

考虑一对相邻的数 \((a_i,a_j)\) 的贡献，假设 \(a_j>a_i\)，那么 \(a_{j}-a_{i}=\sum_{p=i}^{j-1}a_{p+1}-a_p\)。

也就是说，\(a_{p+1}-a_{p}\) 的代价会被所有相邻的数中，一个 \(\ge a_{p+1}\)，另一个 \(\le a_p\) 的位置计算。也就是说，考虑 \(a_{p+1}-a_p\) 的贡献时，只需要考虑 \(\le a_p\) 的连续段和 \(\ge a_{p+1}\) 的连续段交替的次数。

考虑一个 连续段DP，设 \(f_{i,j,k,d}\) 表示：考虑了前 \(i\) 个数，有 \(j\) 个连续段，之前的 \(a_{i}-a_{i-1}\) 的代价之和是 \(k\)，有 \(d\) 个线段占据了两端（\(d\le 2\)）。

因为连续段交替的次数与 \(d\) 有关，也就是 \(2\times j-d\)。那么新的贡献是 \((2\times j-d)\times (a_{i+1}-a_{i})\)。

考虑转移：

• \(f_{i+1,j+1,k',d}\gets f_{i,j,k,d}\times (j+1-d)\)，自己作为新的一段，插入到不为边界的空隙中
• \(f_{i+1,j+1,k',d+1}\gets f_{i,j,k,d}\times (2-d)\)，自己作为新的一段，插入边界中（\(d<2\)）
• \(f_{i+1,j.k',d}\gets f_{i,j,k,d}\times (2\times j-d)\)，自己插入到某一段的旁边
• \(f_{i+1,j-1,k',d}\gets f_{i,j,k,d}\times (j-1)\)，合并两个段（\(j>1\)）
• \(f_{i+1,j,k',d+1}\gets f_{i,j,k,d}\times (2-d)\)，原本不是端点的一段填到端点

最后的答案是 \(\sum f_{n,1,i,2}\)

时间复杂度是 \(O(n^2L)\)。注意特判 \(n=1\)。

view code

#include <bits/stdc++.h>
using namespace std;
const int N=105,M=1005,mod=1e9+7;
inline void Add(int &a,int b){a+=(a+b>=mod)?b-mod:b;}
int n,l,a[N];
int f[2][N][M][3];
int main(){
	scanf("%d %d",&n,&l);
	if(n==1){
		puts("1");
		return 0;
	}
	for(int i=1;i<=n;++i)scanf("%d",&a[i]);
	sort(a+1,a+1+n);
	int pre=0,cur=1;
	f[pre][0][0][0]=1;
	for(int i=1;i<=n;++i,pre^=1,cur^=1){
		memset(f[cur],0,sizeof(f[cur]));
		for(int k=0;k<=l;++k)for(int j=0;j<=i;++j){
			for(int d=0;d<=2;++d){
				int nxt=k+(2*j-d)*(a[i]-a[i-1]),v=f[pre][j][k][d];
				if(nxt>l||!v)continue;
				Add(f[cur][j+1][nxt][d],1ll*v*(j+1-d)%mod);
				if(j>1)Add(f[cur][j-1][nxt][d],1ll*v*(j-1)%mod);
				if(j)Add(f[cur][j][nxt][d],1ll*v*(2*j-d)%mod);
				if(d<2)Add(f[cur][j+1][nxt][d+1],1ll*v*(2-d)%mod);
				if(d<2&&j)Add(f[cur][j][nxt][d+1],1ll*v*(2-d)%mod);
			}
		}
	}
	int ans=0;
	for(int i=0;i<=l;++i)Add(ans,f[pre][1][i][2]);
	printf("%d\n",ans);
	return 0;
}