orangejuice's blog

挂了请务必叫我

Posted on Edited on In

TopCoder - 12349 SRM579 Round1 Div1 RockPaperScissors (概率dp)

题目大意：

有 \(n\) 个骰子，每个骰子有300个面，其中有 \(a_i,b_i,c_i\) 分别为石头/布/剪刀

每轮你选择出石头/剪刀/布，然后会从剩下的骰子中随机取一个再随机结果，但是你不知道取的是什么骰子

赢一局的权值为3，平局为1

求最优情况下最大的权值期望

\[ \ \]

题目分析：

显然当前的局面只和已经抽出的石头/布/剪刀的数量有关（因为这是你唯一的决策依据，也是唯一影响局面的）

乍一看非常抽象的决策过程，实在无法通过分析得到每一种局面下应该作出的决策

于是考虑能否直接把每一种局面出现的概率求出，决策后将期望线性相加

不妨把随机取出的骰子放在一个序列上，一个局面是由已经出现的骰子和当前第 \(i\) 次决策的骰子组成的

令 \(dp_{i,j,k,typ}\) 表示已经出现的骰子中出现 \(i,j,k\) 个石头/布/剪刀，当前决策时骰子结果为 \(typ\) 的概率

实际在 \(dp\) 时， \(typ\) 一维需要额外加入一个值表示未确定下一个是什么

考虑对于每个骰子，枚举它是已经出现/正在决策/在第 \(i\) 次决策之后出现，将其概率累和

\(dp\) 状态为 \(O(n^3)\) ，转移次数为 \(O(n)\) ,复杂度为 \(O(n^4)\)

最后对于每种局面计算最优的决策即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include<bits/stdc++.h>
using namespace std;

#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)); }

class RockPaperScissors {
private:
	static const int N=53;
	static const int eps=1e-9;
	double F[N][N][N][4],G[N][N][N][4],w[3];
	double C[N][N];

public:
	
	double bestScore(vector <int> w1, vector <int> w2, vector <int> w3) {
		int n=w1.size();
		rep(i,0,n) rep(j,C[i][0]=1,i) C[i][j]=C[i-1][j-1]+C[i-1][j];
		memset(F,0,sizeof F),F[0][0][0][3]=1;
		rep(i,1,n) {
			w[0]=w1[i-1]/300.0; w[1]=w2[i-1]/300.0; w[2]=w3[i-1]/300.0;
			rep(a,0,i) rep(b,0,i-a) rep(c,0,i-a-b) rep(d,0,3) G[a][b][c][d]=F[a][b][c][d];
			rep(a,0,i) rep(b,0,i-a) rep(c,0,i-a-b) rep(d,0,3) if(G[a][b][c][d]>eps) {
				// 枚举这个骰子在之前出现过了
				F[a+1][b][c][d]+=G[a][b][c][d]*w[0];
				F[a][b+1][c][d]+=G[a][b][c][d]*w[1];
				F[a][b][c+1][d]+=G[a][b][c][d]*w[2];
				// 枚举这个骰子为下一个出现的
				if(d==3) rep(e,0,2) F[a][b][c][e]+=G[a][b][c][d]*w[e];
			}
		}
		double ans=0;
		rep(a,0,n-1) rep(b,0,n-1-a) rep(c,0,n-1-a-b) {
			// 决策这一轮出什么
			double ma=0;
			rep(d,0,2) cmax(ma,F[a][b][c][d]+F[a][b][c][(d+2)%3]*3);
			ans+=ma/C[n][a+b+c]/(n-a-b-c);
		}
		return ans;
	}
};

Posted on Edited on In

SweetFruits TopCoder - 12141(Matrix-Tree)

问题看起来很复杂，不可写，所以先考虑分解一下

假设最后生效的点集为 \(V\) ，那么答案只和 \(\sum sweetness[V_i]\) 和 \(|V|\) 有关

所以可以考虑对于每一种 \(|V|\) ，先预处理出方案数

得知每一种 \(|V|\) 的方案数之后，可以用 \(\text{meet in the middle}\) 法枚举得到 \(\sum sweetness[V_i]\leq maxsweetness\) 的方案数

方案数是有限制的生成树个数，所以考虑用 \(\text{Matrix-Tree}\) 求

限定有 \(a\) 个点生效， \(b\) 个点不生效， \(c\) 个点是 \(-1\)

那么可能出现的边是 \(a-a,a-c,b-c\) 三种

但是我们无法保证 \(a\) 中的点一定生效，所以可以用容斥/二项式反演得到

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#include<bits/stdc++.h>
using namespace std;

#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))

#define reg register
typedef long long ll;
#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)

#define pb push_back
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;
	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=50,P=1e9+7;

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 n,m;

struct Mat{
	int a[N][N];
	Mat(){ memset(a,0,sizeof a); }
	int Det(int n){  // 用于Matrix-Tree的
		static int G[N][N];
		rep(i,1,n) rep(j,1,n) G[i][j]=a[i][j];
		int res=1;
		rep(i,1,n) {
			if(!G[i][i]) rep(j,i+1,n) if(G[j][i]){ res=P-res; swap(G[i],G[j]); break; }
			if(!G[i][i]) continue;
			ll Inv=qpow(a[i][i]);
			rep(j,i+1,n) {
				ll t=a[j][i]*Inv%P;
				rep(k,i,n) a[j][k]=(a[j][k]-1ll*a[i][k]*t%P+P)%P;
			}
		}
		rep(i,1,n) res=1ll*res*a[i][i]%P;
		return res;
	}
};

int w[N],G[N][N],C[N][N];
void Init(){
	rep(i,0,N-1) rep(j,C[i][0]=1,i) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
	rep(i,0,n) {
		Mat T;
		rep(j,1,m) { // -1
			rep(k,1,n+m) T.a[j][k]=-1;
			T.a[j][j]=n+m-1;
		}
		rep(j,m+1,m+i) { // 生效
			rep(k,1,m+i) if(j!=k) T.a[j][k]=-1;
			T.a[j][j]=m+i-1;
		}
		rep(j,m+i+1,m+n) { // 不生效
			rep(k,1,m) T.a[j][k]=-1;
			T.a[j][j]=m;
		}
		w[i]=T.Det(n+m-1);
	}
	rep(i,0,n) rep(j,0,i-1) w[i]=(w[i]-1ll*C[i][j]*w[j]%P+P)%P; // 容斥/二项式反演
}

int val[N];
vector <int> st[N/2];


int Meet_In_The_Middle(int Max){
	int ans=0,mid=n/2;
	rep(i,0,mid) st[i].clear();
	rep(S,0,(1<<mid)-1) {
		int c=0,s=0;
		rep(i,0,mid-1) if(S&(1<<i)) s+=val[i],c++;
		if(s<=Max) st[c].pb(s);
	}
	rep(i,0,mid) sort(st[i].begin(),st[i].end());
	rep(S,0,(1<<(n-mid))-1) {
		int c=0,s=0;
		rep(i,0,n-mid-1) if(S&(1<<i)) s+=val[i+mid],c++;
		if(s>Max) continue;
		rep(j,0,mid) {
			int x=upper_bound(st[j].begin(),st[j].end(),Max-s)-st[j].begin();
			if(x) ans=(ans+1ll*x*w[c+j])%P;
		}
	}
	return ans;
}

class SweetFruits {
public:
	int countTrees(vector <int> Val, int Max) {
		sort(Val.begin(),Val.end(),greater <int> ());
		m=0; while(Val.size() && *Val.rbegin()==-1) m++,Val.pop_back();
		n=Val.size();
		rep(i,0,n-1) val[i]=Val[i];
		Init();
		return Meet_In_The_Middle(Max);
	}
};


Posted on Edited on In

[TopCoder - 12244 SRM 559 Round1 Div1] CircusTents

小而精的计算几何题

题目大意：有 \(n\) 个实心圆（不能从内部经过）

在第一个圆上选出一个点，使得从其他任意圆上到达它的最小距离 最大

分析：要最小值最大，显然可以想到二分答案

不能穿过其他圆这一条件让计算答案变得十分困难，但是可以发现，如果路径经过了一号圆以外的圆

那么从路径与该圆的交点直接过去的距离一定更近

也就是说，可以把 距离 看做 可以穿过一号圆以外的圆 的距离

考虑从一个圆 \(O\) 到达一号圆上的某一点 \(X\) 的最小距离，大致可以成两种情况

  1. \(OX\) 连线不穿过一号圆，那么可以直接走 \(OX\) 连线
QQ截图20201102162335.png

最优路径就是绿色线

2.先走一条切线，然后绕着圆周走一段

QQ截图20201102162651.png

其中 \(Y\) 是 \(O\) 点对于一号圆的切线的切点，绿色线+圆弧是最优路径

那么二分答案 \(mid\) 之后，可以发现满足距离 \(\ge mid\) 的选点位置是一段圆弧，可以从角度取交集判断是否有解

实现上，可以先把一号圆平移到远点

然后对于其他的圆，按照角度范围是否在切线内部可以分类讨论

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include<bits/stdc++.h>
using namespace std;
typedef double db;
#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,const T &b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,const T &b){ ((a<b)&&(a=b)); }

const int N=110;
const db eps=1e-10,Pi=acos(-1),D=2*Pi;

int n;
struct Cir{
    int x,y,r;
}A[N];

db X[N],Y[N],H[N];
int C,S[N];
db Norm(db x){ return x<0?x+D:x; }
int I(db x){ return lower_bound(H+1,H+C+1,x)-H; }

int Check(db L){
    memset(S,0,sizeof S);
    H[1]=0,H[C=2]=D;
    rep(i,1,n) {
        db dis=A[i].x*A[i].x+A[i].y*A[i].y;
        db t=sqrt(dis-A[0].r*A[0].r)-A[i].r;
        //t为走过切线的距离
        dis=sqrt(dis);
        if(dis-A[0].r-A[i].r>=L){ X[i]=0,Y[i]=D; continue; }
        db l,r;
        if(t>=L) {
            // 说明范围在切线位置以内
            // 此时，满足d=dis((x0,y0),A[i].O)-A[i].r>=L
            db a=dis,b=A[0].r,c=L+A[i].r;
            db co=(a*a+b*b-c*c)/(2*a*b);
            //余弦定理
            db x=acos(co),y=atan2(A[i].y,A[i].x);
            l=y+x,r=y-x+D;
        } else {
            db d=acos(A[0].r/dis);
            db x=atan2(A[i].y,A[i].x);
            db y=(L-t)/A[0].r+d;
            // 圆弧长度/半径得到圆弧弧度
            if(y>Pi) return 0; // 特判一下全部覆盖的情况
            l=x+y,r=x-y+D;
        }
        // 求 [l,r] 的交
        if(r>D) l-=D,r-=D;
        if(r<=0) l+=D,r+=D;
        H[++C]=Norm(X[i]=l),H[++C]=Y[i]=r;
    }

    sort(H+1,H+C+1);
    rep(i,1,n) {
        if(X[i]>=0) S[I(X[i])]++,S[I(Y[i])]--;
        else S[I(0)]++,S[I(Y[i])]--,S[I(X[i]+D)]++,S[I(D)]--;
    }
    // 暴力求交
    rep(i,1,C) if((S[i]+=S[i-1])==n) return 1;
    return 0;
}

class CircusTents {
    public:
        double findMaximumDistance(vector <int> _x, vector <int> _y, vector <int> _r) {
            n=_x.size()-1;
            rep(i,0,n) A[i]=(Cir){_x[i]-_x[0],_y[i]-_y[0],_r[i]};

            db l=0,r=1e5;
            while(r-l>eps) {
                db mid=(l+r)/2;
                if(Check(mid)) l=mid;
                else r=mid;
            }
            return l;
        }
};

Posted on Edited on In

【UER #9】知识网络

bitset写错没调出来。。。

\(O(n(n+m))\)

暴力枚举起点，建立转移虚点，得到一个边权为 \(0/1\) ，点数 \(n+k\) ，边数为 \(O(n+m)\) 的图

然后广搜双端队列维护即可

\[ \ \]

\(O(k(n+m)+\frac{n(n+m)}{w})\)

考虑枚举颜色 \(k\) ，对于所有这种颜色的点，假设一开始 \(dis_u\) 均为 \(2\)

并由此广搜预处理一个最短路图，复杂度为 \(O(k(n+m))\)

由于图上有0边无0环，因此最短路图是拓扑图

那么选定其中一个点为起点，会将这个点的 \(dis\rightarrow dis-1\) ，其他同色点 \(dis\) 不变

同时，这个点在最短路图上的所有后记节点 \(dis\rightarrow dis-1\)

那么对于最短路图上所有点，统计有多少个起点能够到达它，就能够知道以这种颜色点为起点时，这个点不同的 \(dis\) 出现次数

这是一个拓扑图 \(dp\) 问题，不好处理，因此考虑用 \(\text{bitset}\) 暴力存储所有能够转移的状态

由于需要维护的起点总数为 \(n\) ，单次只维护一个集合，因此无法直接用 \(\text{std::bitset}\)

复杂度为 \(O(\frac{nm}{w})\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef pair <int,int> Pii;
#define reg register
#define mp make_pair
#define pb push_back
#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;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) f|=IO=='-';
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=51000,INF=1e9+10;

int n,m,k;
struct Edge{
	int to,nxt,w;
} e[N*4];
int head[N],ecnt;
void AddEdge(int u,int v,int w=1){
	e[++ecnt]=(Edge){v,head[u],w};
	head[u]=ecnt;
}

int Q[N*4],D[N*4],L,R,dis[N];
unsigned Ans[N];
vector <int> G[200];
int len;
struct BITSET{
	ull a[N/128];
	int count(){
		int c=0;
		rep(i,0,len) c+=__builtin_popcountll(a[i]);
		return c;
	}
	void set(int x){ a[x>>6]|=1ull<<x; }
	void clear(){ memset(a,0,(len+1)<<3); }
	void copy(const BITSET &t){ memcpy(a,t.a,(len+1)<<3); }
	void operator |= (const BITSET &t) { rep(i,0,len) a[i]|=t.a[i]; }
} BS[N];

int cnt,vis[N],ind[N];
void Topo(int sz){
	rep(i,1,n+k) ind[i]=0;
	rep(u,1,n+k) for(int i=head[u];i;i=e[i].nxt) if(dis[e[i].to]==dis[u]+e[i].w) ind[e[i].to]++;
	L=1,R=0;
	rep(i,1,n+k) if(dis[i]<1e6 && !ind[i]) Q[++R]=i;
	while(L<=R){
		int u=Q[L++];
		if(u<=n) {
			int t=BS[u].count();
			Ans[dis[u]]+=sz-t;
			Ans[dis[u]-1]+=t;
		}
		for(int i=head[u];i;i=e[i].nxt) if(dis[e[i].to]==dis[u]+e[i].w) {
			BS[e[i].to]|=BS[u];
			if(--ind[e[i].to]==0) Q[++R]=e[i].to;
		}
	}
}

void Bfs(int c){
	if(!G[c].size()) return;
	memset(dis,63,(n+k+2)<<2);
	L=2*(n+k),R=2*(n+k)-1;
	rep(i,0,G[c].size()-1) Q[++R]=G[c][i],D[R]=dis[G[c][i]]=2;
	while(L<=R){
		int u=Q[L++];

		for(reg int i=head[u];i;i=e[i].nxt){
			reg int v=e[i].to,w=e[i].w;
			if(dis[v]<=dis[u]+w) continue;
			dis[v]=dis[u]+w;
			if(!w) Q[--L]=v,D[L]=dis[v];
			else Q[++R]=v,D[R]=dis[v];
		}
	}
	len=G[c].size()/64;
	if(len>=N/128){
		len=N/128-1;
		rep(i,1,n+k) BS[i].clear();
		rep(i,0,G[c].size()-1) if(i<len*64) BS[G[c][i]].set(i);
		Topo(len*64);

		int tmp=len*64;
		len=(G[c].size()-tmp)/64;
		rep(i,1,n+k) BS[i].clear();
		rep(i,0,G[c].size()-1) if(i>=tmp) BS[G[c][i]].set(i-tmp);
		Topo(G[c].size()-tmp);
	} else {
		rep(i,1,n+k) BS[i].clear();
		rep(i,0,G[c].size()-1) BS[G[c][i]].set(i);
		Topo(G[c].size());
	}
}

int col[N];

int main(){
	n=rd(),m=rd(),k=rd();
	rep(i,1,n) {
		int x=col[i]=rd();
		AddEdge(n+x,i,0),AddEdge(i,n+x);
		G[x].pb(i);
	}
	rep(i,1,m) {
		int u=rd(),v=rd();
		if(col[u]==col[v]) continue;
		ind[u]++,ind[v]++;
		AddEdge(u,v),AddEdge(v,u);
	}
	rep(i,1,k) Bfs(i);
	Ans[1]=0;
	rep(i,2,k*2) Ans[i]/=2;
	Ans[k*2+1]=1ll*n*(n-1)/2;
	rep(i,1,k*2+1) printf("%u ",Ans[i]),Ans[k*2+1]-=Ans[i];
}

Posted on Edited on In

【UER #9】赶路

-----一定有解。。

\(x_1\leq x_i\leq x_n\)

将中间的点按照 \((x_i,y_i)\) 排序，然后依次连过去即可

\[ \ \]

\(x_1=y_1=0\) ，四个象限均存在点

将所有点极角排序，然后走一圈即可

\[ \ \]

\(O(n\log n)\)

不妨设 \(x_1<x_n\)

将 \(1,n\) 以外所有点分成三部分，即左边| 1 | 中间 | n | 右边

左边右边考虑极角排序转圈走，中间按照 \((x_i,y_i)\) 走

发现两边极角排序之后转圈走不一定能够走到中间去，可能会与转圈时的路径相交

但是实际上画图就会发现，如果顺时针走的路径会相交，逆时针走一定不相交

因此对于左右枚举顺时针还是逆时针即可，4中情况，每种 \(O(n)\) 检查线段相交

总能构造一组合法解

预处理需要排序，因此复杂度为 \(O(n\log n)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef pair <int,int> Pii;
#define reg register
#define mp make_pair
#define pb push_back
#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;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) f|=IO=='-';
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=510,INF=1e9+10;
const db eps=1e-7;

int n,m;
struct Node{
	db x,y;
	Node(){}
	Node(db x,db y):x(x),y(y){ }
	bool operator < (const Node __) const { return x!=__.x?x<__.x:y<__.y; }
	Node operator - (const Node &t) const { return Node(x-t.x,y-t.y); }
	db operator * (const Node &t) const { return x*t.x+y*t.y; }
	db operator ^ (const Node &t) const { return x*t.y-y*t.x; }
	void turn(db t){
		db a=x*cos(t)-y*sin(t),b=x*sin(t)+y*cos(t);
		x=a,y=b;
	}
	db tan() const { return atan2(x,y); }
} A[N];
int Cross(int x,int y,int a,int b){
	if(((A[x]-A[y])^(A[a]-A[b]))==0) return 0;
	db l=(A[a]-A[x])^(A[y]-A[x]),r=(A[b]-A[x])^(A[y]-A[x]);
	if((l<0)^(r>0)) return 0;
	swap(a,x),swap(b,y);
	l=(A[a]-A[x])^(A[y]-A[x]),r=(A[b]-A[x])^(A[y]-A[x]);
	if((l<0)^(r>0)) return 0;
	return 1;
}

int L[N],LC,R[N],RC,T[N],TC;
int P[N];

int Work(){
	int C=0;
	P[++C]=1;
	rep(i,1,LC) P[++C]=L[i];
	rep(i,1,TC) P[++C]=T[i];
	rep(i,1,RC) P[++C]=R[i];
	P[++C]=n;
	if(LC) rep(i,2,LC) if(Cross(P[LC+1],P[LC+2],P[i],P[i-1])) return 0;
	if(RC) drep(i,n,n-RC+2) if(Cross(P[n-RC],P[n-RC-1],P[i],P[i-1])) return 0;
	return 1;
}
void Solve(){
	rep(i,0,1) {
		rep(j,0,1) {
			if(Work()) {
				rep(i,1,n) printf("%d ",P[i]);
				puts("");
				return;
			}
			reverse(R+1,R+RC+1);
		}
		reverse(L+1,L+LC+1);
	}
}

int main(){
	rep(kase,1,rd()){
		rep(i,1,n=rd()) A[i].x=rd(),A[i].y=rd();
		db t=1.0*(rand()+2)/(rand()+2);
		rep(i,1,n) A[i].turn(t);
		if(A[1].x>A[n].x) rep(i,1,n) A[i].x=-A[i].x;
		LC=RC=TC=0;
		rep(i,2,n-1) if(A[i].x<A[1].x-eps) L[++LC]=i;
		else if(A[i].x-eps>A[n].x) R[++RC]=i;
		else T[++TC]=i;
		sort(T+1,T+TC+1,[&](int x,int y){ return A[x].x<A[y].x; });
		sort(L+1,L+LC+1,[&](int x,int y){ return (A[x]-A[1]).tan()<(A[y]-A[1]).tan(); });
		sort(R+1,R+RC+1,[&](int x,int y){ return (A[x]-A[n]).tan()<(A[y]-A[n]).tan(); });
		Solve();
	}
}

Posted on Edited on In

TopCoder - 12584 SRM 582 Div 1 SemiPerfectPower (莫比乌斯反演)

题目大意：

给定 \(L,R\) ，求 \([L,R]\) 中能够表示为 \(a\cdot b^c(1\leq a<b,c>1)\) 的数(SemiPerfect数)的个数

\(R\leq 8\cdot 10^{16}\)

解题思路

首先显然可以通过作差转化为求 \([1,n]\) 内的个数

接下来考虑简化 \(c\) 的情形

推论：任何一个SemiPerfect数可以表示为 \(c\leq 3\) 的形式

证明：若 \(c>3\) ，对于 \(n=a\cdot b^c(c>3)\)

当 \(2|c\) 时，显然存在形如 \(n=a\cdot (b^{\frac{c}{2}})^{2}\) 的表示

当 \(2\not |c\) 时，可以表示为 \(n=(a\cdot c)\cdot (b^{\frac{c-1}{2}})^2\) 同样合法

接下来考虑对于两种情况分类讨论

为了便于叙述，令 \(F(n)=\max \lbrace k\in \N^+|\ k^2|n\rbrace\)

\(G(n)=[\nexists k>1,k^3|n]\)

Part1 \(c=2\)

为了避免重复，强制每一个数 \(n\) 的唯一表示为

\(n=a\cdot b^2(F(a)=1,a<b)\)

由于 \(a<b\) ，所以显然 \(n>a^3\) ，即 \(a<n^{\frac{1}{3}}\)

暴力枚举 \(a\) ，预处理 \(n^{\frac{1}{3}}\) 中所有的 \(G(i)\) 即可

\[ \ \]

Part \(c=3\)

同样的，限制条件为 \(n=a\cdot b^3,(G(a)=1,a<b)\) ，得到 \(a<n^{\frac{1}{4}}\)

但是由于 \(c=2,3\) 两部分有重复，还需额外考虑强制 \(n\) 不存在形如 \(n=a'\cdot b'^2\) 的表示

假设已知 \(n=a\cdot b^3\) ，不存在 \(n=a'\cdot b'^2\) 的判定条件是

\(a'=\frac{n}{F^2(n)}\ge F(n)\) ，即 \(F(n)\leq n^{\frac{1}{3}}\)

同时由于 \(F(n)=F(a\cdot b)b\)

得到 \(F(a,b)\leq a^{\frac{1}{3}}\)

由于等号右边包含 \(a\) ，考虑枚举 \(a\) ，易求出 \(L=F(a),d=\frac{a}{L^2}\) ，得到 \(F(a\cdot b)\) 的另一种表达形式

\(F(a\cdot b)=L \cdot \gcd (d,b)\cdot F(\frac{b}{\gcd(d,b)})\leq a^{\frac{1}{3}}\)

上面的转化意为： \(L\) 为 \(a\) 中已经成对的部分自然取出，然后优先考虑为 \(D\) 匹配 \(b\) 中的因数成对，对于剩下的部分再重新计算答案

\[ \ \]

化简该式，得到 \(L\cdot \gcd(d,b)F(\frac{b}{\gcd(d,b)})\leq a^{\frac{1}{3}}\)

式子包含 \(\gcd\) ，似乎具有莫比乌斯反演的性质

考虑计算 \(b\in [a+1,(\frac{n}{a})^{\frac{1}{3}}]\) 的数量

观察到 \(a^{\frac{1}{3}}\leq n^{\frac{1}{12}}\) ，上限只有 \(25\) 左右，可以考虑直接枚举 \(F(\frac{b}{\gcd(b,d)})\)

令枚举的 \(g=\gcd(b,d),F(\frac{b}{g})=f\) ，计算 \(\gcd(\frac{b}{g},\frac{d}{g})=1,g\cdot f\cdot L\leq a^{\frac{1}{3}}\) 的 \(b\) 的数量

考虑直接从 \(\frac{b}{g}\) 中把 \(f^2\) 的因数提取出来，令 \(\begin{aligned} L'=\lceil \frac{a+1}{gf^2}\rceil,R'=\lfloor \frac{(\frac{n}{a})^{\frac{1}{3}}}{gf^2}\rfloor \end{aligned}\) ，令 \(i=\frac{b}{gf^2},x=\frac{d}{g}\) ，得到新的限制条件式子为

\(\gcd(x,f)=1,\gcd(i,x)=1,F(i)=1,i\in[L',R']\)

在确定了 \(g,f\) 之后，需要考虑的限制就是 \(\gcd(i,x)=1,F(i)=1,i\in[L',R']\)

由于包含 \(\gcd\) ，不妨用莫比乌斯反演计算该式，得到表达式为

\(Sum=\sum_{k|x}\mu(k)\sum_{i\in [L',R']} [k|i\ \text{and}\ F(i)=1]\)

对于 \(k\) ，计算 \(\sum_{i\in[L',R']}[k|i\ \text{and}\ F(i)=1]\) 可以归纳为计算

\(\sum_{i=\frac{n}{k}} [F(ik)=1]\) ，一共有 \(n^{\frac{1}{3}}\ln n\) 种不同的权值，可以暴力预处理得到

枚举 \(d\) 的因数对于所有的上面的式子计算，可能的 \(g,f\) 并不多，可以直接枚举

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define pb push_back
#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)

class SemiPerfectPower {
public:
	static const int N=450000,M=20000;
	int w[N],notpri[N],pri[N],pc,F[N],G[N];
	vector <int> S[N],Fac[N];
	ll gcd(ll a,ll b){ return b==0?a:gcd(b,a%b); }
	SemiPerfectPower(){
		// 预处理F,G ,w[i]=\mu(i) S[k][i]=\sum_{j \in [1,i]}  F(i,k)=1
		rep(i,1,sqrt(N)) rep(j,1,(N-1)/i/i) F[i*i*j]=i;
		rep(i,2,pow(N,1/3.0)) rep(j,1,(N-1)/i/i/i) G[i*i*i*j]=1;
		rep(i,1,M-1) for(int j=i;j<M;j+=i) Fac[j].pb(i);
		w[1]=1;
		rep(i,2,N-1) {
			if(!notpri[i]) pri[++pc]=i,w[i]=-1;
			for(int j=1;i*pri[j]<N && j<=pc;++j) {
				notpri[i*pri[j]]=1;
				if(i%pri[j]==0) {
					w[i*pri[j]]=0;
					break;
				} 
				w[i*pri[j]]=-w[i];
			}
		}
		rep(i,1,N-1) if(w[i]) {
			S[i].resize(N/i+1);
			rep(j,1,(N-1)/i) S[i][j]=S[i][j-1]+(F[i*j]==1);
		}
	}

	ll Solve2(ll n){
		ll ans=0,UX=pow(n,1/3.0);
        // 防止浮点误差
		if((UX+1)*(UX+1)*(UX+1)<=n) UX++;
		if(UX*UX*UX>n) UX--;
		rep(i,1,UX) if(F[i]==1) {
			ll UY=sqrt(n/i);
            // 防止浮点误差
			if(i*(UY+1)*(UY+1)<=n) UY++;
			if(i*UY*UY>n) UY--;
			ans+=max(0ll,UY-i);
		}
		return ans;
	}

	ll Solve3(ll n){
		ll UX=pow(n,0.25); 
		// 枚举c的上界
		ll ans=0;
		if((UX+1)*(UX+1)*(UX+1)*(UX+1)<=n) UX++;
		rep(x,1,UX) if(!G[x]) {
			ll UY=pow(n/x,1.0/3.0),U=pow(x,1/3.0);
			
			// 防止浮点误差
			if((U+1)*(U+1)*(U+1)<=x) U++;
			if(U*U*U>x) U--;
			if(x*(UY+1)*(UY+1)*(UY+1)<=n) UY++;
			if(x*UY*UY*UY>n) UY--;

			int L=F[x],D=x/L/L;
			for(int G:Fac[D]) {
				for(int g:Fac[G]) {
					if(g*L>U) break;
					rep(f,1,U/g/L) if(gcd(f,D/g)==1) {
						int L=x/G/f/f,R=UY/G/f/f;
						ans+=w[G/g]*(S[G/g][R]-S[G/g][L]);
					}
				}
			}
		}
		return ans;
	}
	ll Solve(ll n) {
		return Solve2(n)+Solve3(n);
	}
	ll count(ll L,ll R) {
		return Solve(R)-Solve(L-1);
	}
};

Posted on Edited on In

Topcoder SRM 569 Div1 - MegaFactorial (矩阵)

首先是对于末尾0个数的处理，设最后得到的数中包含 \(i\) 的指数为 \(F(i)\)

对于 \(B=2,3,5,7\) 的情况，可以直接计算答案 \(\sum_{i=1}\sum_{j=1}F(j\cdot B^i)\)

对于 \(B\) 为质因子组合的情况，即 \(B=6(2\times 3),10(2\times 5)\) ，因为 \(F(i)\) 实际有单调性，可以直接取较大的质因子

对于 \(B\) 为质因子次方的情况，即 \(B=2^2,2^3,3^2\) 的情况，设 \(B=p^k\) 则答案可以表示为

\(\begin{aligned} \lfloor \frac{\sum_{i=1}\sum_{j=1}F(j\cdot p^i)}{k}\rfloor \end{aligned}\)

由于要取模，实际上要做一点魔改，设模数为 \(m\) ，答案可以表示为 \(c=ak+b\) 的形式，这个式子求出的是 \(a\)

则 \(\lfloor \frac{(ak+b)\mod km }{k}\rfloor =a+\lfloor \frac{(b\mod km)}{k}\rfloor =a\)

由于 \(k\leq 3\) ，扩大模数后可以用unsigned int 存

下面考虑用矩阵求解上式

\(nk!\) (下面用 \(f(n,k)\) 表示)这个东西可以看作从 \(n\) 向下的一个递推式

因此考虑以 \(k\) 为矩阵元素，求出每个 \(f(n,k)\) 被调用的次数

注意这样递推就是反向的了

递推的转移式子是 \(f(n,k)\rightarrow f(n,k-1),f(n-1,k)\) ，其中 \(f(n,k-1)\) 的转移需要在层内完成

据此构造矩阵即可，注意 \(f(n,0)\) 不能向 \(f(n-1,0)\) 转移

考虑对于 \(\sum_{i=1}\sum_{j=1}F(j\cdot B^i)\) 的每个 \(i\) 求解，一共有 \(\frac{n}{B^i}\) 个 \(j\) ，每个 \(j\) 出现的递推层数为等差数列

即 \(n\mod B^i,n\mod B^i+B^i,n\mod B^i +2 \cdot B^i\cdots\)

我们要求的其实是每一层的 \(f(i,0)\) ，所以考虑求出每次 \(B^i\) 层的转移矩阵

然后是依次累和，把矩阵的转移中 \(0\rightarrow 0\) 的转移赋为1即可做到

tips:首项是 \(n\mod B^i\)

一共有 \(\log _Bn\) 种不同的 \(i\) ，因此复杂度为 \(O(\log_Bn\log_2 n\cdot k^3)\)

当然更优的做法，咕咕咕

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include<bits/stdc++.h>
using namespace std;

#define reg register
typedef long long ll;
typedef unsigned int U;
#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)

#define pb push_back
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;
	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=20;


int n,d;
U P=1e9+9;
struct Mat{
	U a[N][N];
	Mat(){ memset(a,0,sizeof a);} 
	void One(){ rep(i,0,d) a[i][i]=1; }
	U* operator [] (int x){ return a[x]; }
	Mat operator * (const Mat &x) const {
		Mat res;
		rep(i,0,d) rep(j,0,d) rep(k,0,d) res.a[i][k]=(res.a[i][k]+1ll*a[i][j]*x.a[j][k])%P;
		return res;
	}
	void Show(){
		rep(i,0,d) { rep(j,0,d) printf("%d ",a[i][j]); puts(""); }
	}
} A,B,C;

Mat qpow(Mat x,int k){
	Mat res; res.One();
	for(;k;k>>=1,x=x*x) if(k&1) res=res*x;
	return res;
}
int Factor(int &x) {
	int p=-1,c=0;
	rep(i,2,x) if(x%i==0) {
		while(x%i==0) c++,x/=i;
		p=i;
		break;
	}
	return x=p,c;
}

class MegaFactorial {
	public:
		int countTrailingZeros(int N, int K, int b) {
			A=Mat(),n=N,d=K;

			if(b==10) b=5;
			if(b==6) b=3;
			int t=Factor(b);
			P*=t;

			drep(i,d,0) {
				A[i][i]=1;
				rep(j,0,d) A[j][i]+=A[j][i+1];
			}
			A[0][0]=0;
			ll ans=0;
			for(ll i=b;i<=n;i*=b) {
				B=qpow(A,i); B[0][0]=1;
				Mat res=qpow(A,n%i)*qpow(B,n/i-1);
				rep(i,0,d) (ans+=res[i][0])%=P;
			}
			P/=t,ans/=t;

			return ans;
		}
};


Posted on Edited on In

[BJ United Round #3] 押韵

先%%%%%%%%%%%%%%%%% EI

\[ \ \]

\[ \ \]

下文默认模数为 \(P\)

简要题意：求：用 \(k\) 种元素，每种元素使用 \(d\) 的倍数次，排成一个长度为 \(nd\) 的序列 的方案数

这个题目的设定就让人想到两个离不开的元素 ： (模数暗示了？)

指数型生成函数 + 单位根反演

显然可以得到每一种元素的指数型生成函数为

\(\begin{aligned} \text{EGF(Element)}=F(x)=\sum_{d|i} \frac{x^i}{i!}\end{aligned}\)

带入单位根反演 \(\begin{aligned}\ [d|n]=\frac{\sum_0^{d-1} \omega_d^{in}}{d}\end{aligned}\)

即 \(F(x)=\begin{aligned}\frac{1}{d}\cdot \sum_{i=0}^{d-1}e^{\omega_d^ix}\end{aligned}\)

而总的生成函数就是 \(G(x)=F^k(x)\)

即 \(\begin{aligned} G(x)=\frac{1}{d^k}\cdot (\sum_{i=0}^{d-1}e^{\omega_d^ix})^k\end{aligned}\)

其中的和式幂次展开会得到一个 \(k^d\) 项的多项式，我们要求 \([x^n]G(x)\) ，就需要展开得到每一项的幂系数

所以显然我们需要先合并同类项一下。。。

而幂系数是一个单位根之和的形式，这就需要我们寻找单位根之间的关系

这里得到一个思路：用 \(d\) 次单位根中的 \(\varphi(d)\) 个作为基底，以简单的 有理数/整系数 表示出所有的 \(\omega_d^i\)

\[ \ \]

对于 \(d=4\) 的情况比较简单， \(\varphi(d)=2\) ，可以得到四个单位根分别为 \(1,\omega,-1,-\omega\)

可以枚举得到的和为 \(x+y\omega\) ，然后求系数

优先考虑组合意义，可以发现就是在平面上每次可以走四个方向， \(k\) 步之后最终到达 \((x,y)\) 的方案数

两个维度分立的情况，还需要枚举每个维度走了几步，所以用一种巧妙的转化两个维度联系在一起

将平面旋转 \(\frac{\pi}{8}\) ，并且扩大 \(\sqrt 2\) 倍，得到新的坐标为 \((x-y,x+y)\) ，新的行走方向是 \((+1,+1),(-1,-1),(-1,+1),(+1,-1)\)

这样以来，每次每个维度都有行走，可以确定每个维度 \(+1\) 和 \(-1\) 的次数，直接组合数排列即可得到答案

\[ \ \]

\[ \ \]

对于 \(d=6\) ，甚至是更一般的情况的情况

只在代数层面来看单位根似乎十分抽象，不如从复平面单位根上找一找灵感

下面是 \(d=6\) 的情形

planeomega.png

\(\varphi(6)=2\) ，假设以 \(\overrightarrow{OA},\overrightarrow{OB}\) 作为基底，可以直观地得到基底表达

\(\overrightarrow{OA}=1\) \(\overrightarrow{OB}=\omega\)
\(\overrightarrow{OA}=\omega_6^0\) 1 0
\(\overrightarrow{OB}=\omega_6^1\) 0 1
\(\overrightarrow{OC}=\omega_6^2\) -1 1
\(\overrightarrow{OD}=\omega_6^3\) -1 0
\(\overrightarrow{OE}=\omega_6^4\) 0 -1
\(\overrightarrow{OF}=\omega_6^5\) 1 -1

由此我们得到了一个 \(\varphi(d)\) 维数的表达方法

把每一维看做不同元，也就是说，得到了一个 \(\varphi(d)\) 维， \(O(1)\) 次的多项式，需要我们求高维多项式幂次

令 \(N=k^{\varphi(d)}\)

直接压位暴力多项式复杂度为 \(O(N\log N-N\log^2N)\) ，而且面临着模数难以处理，常数大的问题

所以 \(\text{EI}\) 又用出了一个巧妙的暴力方法解决这个问题，以 \(d=6\) 为例，先做一下处理，得到要求的多项式

似乎每次 \(k\) 次幂总是求导+递推？

\(f(x,y)=x^2y+xy^2+y^2+y+x+x^2,g(x,y)=f^k(x,y)\)

\(g(x,y)\) 对于 \(x\) 求偏导，得到 \(g'(x,y)=kf^{k-1}(x,y)f'(x,y)\)

即 \(g'(x,y)f(x,y)=kg(x,y)f'(x,y)\)

\(f'(x,y)=2xy+2x+y^2+1\)

然后我们要解这个方程，考虑乘积为 \([x^ny^m]\) 一项两边的系数

左边 \(=[x^{n-2}y^{m-1}]+[x^{n-1}y^{m-2}]+[x^{n}y^{m-2}]+[x^{n}y^{m-1}]+[x^{n-1}y^{m}]+[x^{n-2}y^{m}]\)

换成 \(g(x,y)\) 的系数应该是

左边 \(=(n-1)[x^{n-1}y^{m-1}]+n[x^{n}y^{m-2}]+(n+1)[x^{n+1}y^{m-2}]+(n+1)[x^{n+1}y^{m-1}]+n[x^{n}y^{m}]+(n-1)[x^{n-1}y^{m}]\)

右边 \(=2k[x^{n-1}y^{m-1}]+2k[x^{n-1}y^m]+k[x^ny^{m-2}]+k[x^ny^m]\)

其中 \([x^{n+1}y^{m-1}]\) 只出现了一次，按照先 \(n\) 递增再 \(m\) 递增的顺序进行递推，即

\(\begin{aligned}\ [x^ny^m]=\frac{2k[x^{n-2}y^{m}]+2k[x^{n-2}y^{m+1}]+k[x^{n-1}y^{m-1}]+k[x^{n-1}y^{m+1}]}{n}\end{aligned}\)

\(\begin{aligned}-\frac{(n-2)[x^{n-2}y^{m}]+(n-1)[x^{n-1}y^{m-1}]+n[x^{n}y^{m-1}]+(n-1)[x^{n-1}y^{m+1}]+(n-2)[x^{n-2}y^{m+1}]}{n}\end{aligned}\)

边界条件是 \(\begin{aligned}\ [x^i]=[y^i](i\ge k)=\binom{k}{i-k}\end{aligned}\) (由系数 \(x,x^2\) 或 \(y,y^2\) 得到)

由此带入递推即可

综上，得到的每项的系数的复杂度为 \(O(d\cdot k^{\varphi(d)})\) ，其中 \(d\) 为递推每项需要的时间

由系数得到答案仍然需要一次快速幂，因此依然带一个 \(\log P\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#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=min(a,b); }
template <class T> inline void cmax(T &a,T b){ a=max(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;
}

const int N=2e3+10,P=1049874433,G=7;

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

int main() {
    n=rd(),k=rd(),d=rd();
    w[0]=1,w[1]=qpow(G,(P-1)/d);
    rep(i,2,90) w[i]=1ll*w[i-1]*w[1]%P;
    rep(i,0,k) rep(j,C[i][0]=1,i) C[i][j]=C[i-1][j]+C[i-1][j-1],Mod1(C[i][j]);
    I[0]=I[1]=1;
    rep(i,2,k*2) I[i]=1ll*(P-P/i)*I[P%i]%P;
    if(d==1){
        int ans=qpow(k,1ll*d*n);
        printf("%d\n",ans);
    } else if(d==2) {
        int ans=0;
        rep(i,0,k) ans=(ans+qpow((1ll*w[0]*i+1ll*w[1]*(k-i))%P,1ll*d*n)*C[k][i])%P;
        ans=ans*qpow(qpow(d,k))%P;
        printf("%d\n",ans);
    } else if(d==3) {
        int ans=0;
        rep(i,0,k) rep(j,0,k-i) 
            ans=(ans+qpow((1ll*w[0]*i+1ll*w[1]*j+1ll*w[2]*(k-i-j))%P,1ll*d*n)*C[i+j][i]%P*C[k][i+j])%P;
        ans=ans*qpow(qpow(d,k))%P;
        printf("%d\n",ans);
    } else if(d==4) {
        int ans=0;
        rep(i,-k,k) rep(j,-k,k) if(abs(i)+abs(j)<=k && (k-i-j)%2==0) {
            ll x=qpow((1ll*w[0]*i+1ll*w[1]*j)%P,1ll*d*n);
            ll y=1ll*C[k][(abs(i-j)+k)/2]*C[k][(k+abs(i+j))/2]%P;
            ans=(ans+x*y)%P;
        }
        ans=(ans+P)*qpow(qpow(d,k))%P;
        printf("%d\n",ans);
    } else {
        static int F[N*2][N*2];
        int ans=0;
        rep(i,0,k*2) rep(j,max(k-i,0),min(2*k,3*k-i)) {
            if(i==0) F[i][j]=C[k][j-k];
            else if(j==0) F[i][j]=C[k][i-k];
            else {
                int s=(2ll*k*(i>1?F[i-2][j]+F[i-2][j+1]:0)+1ll*k*(F[i-1][j-1]+F[i-1][j+1]))%P;
                int t=((i>1?1ll*(i-2)*(F[i-2][j]+F[i-2][j+1]):0)+
                    1ll*(i-1)*F[i-1][j-1]+1ll*i*F[i][j-1]+1ll*(i-1)*F[i-1][j+1])%P;
                F[i][j]=1ll*(s-t+P)*I[i]%P;
            }
            ans=(ans+qpow((1ll*w[0]*(i-k)+1ll*w[1]*(j-k))%P,1ll*d*n)*F[i][j])%P;
        }
        ans=(ans+P)*qpow(qpow(d,k))%P;
        printf("%d\n",ans);
    }
}

Posted on Edited on In

[BZOJ1852] [MexicoOI06]最长不下降序列(贪心)

考虑如下贪心

（我将问题反过来考虑，也就是要满足 \(A_i > \max_{j=1}^{j < i}{B_j}\) ）

首先对于读入的 \((A,B)\) ，按照 \(B\) 的值递增排序

(选出的答案序列不一定是其中一个有序的子序列)

答案序列存在若干个 \(B\) 递增的位置，设它们是 \(\{a_i\},a_{i-1}<a_i\)

合法的递增序列需要满足的限制是 \(A_{a_i}>B_{a_{i-1}}\)

考虑剩下的部分即 \(j\in[a_{i-1}+1,a_i-1]\) ，那么这些点放在 \(a_i\) 后面一定是最优的(因为此时不会改变最大的 \(B\) )，此时限制它们的 \(B\) 就是 \(B_{a_i}\)

即这一部分中满足 \(A_j>B_{a_i}\) 的 \(j\) 均可以选出来

为了便于表示，设 \(C(l,r)=|\{i|i\in[l,r],A_i>B_{r+1}\}|\) ，可以通过一个主席树维护

定义 \(dp_i\) 表示当前以 \(i\) 为最大的 \(B\) 的答案，特别的， \(dp_0\) 表示序列为空 \(A_0=B_0=-\infty\)

枚举每个 \(i\) 进行转移

朴素的转移就是可以枚举上一个位置 \(j\) ， \(O(n^2)\) 转移

需要找到前面第一个能把 \(i\) 接上去的 \(j\) 即可，即第一个 \(B_r<A_i(j\ge 0)\) 的位置，那么合法的决策位置就是 \([0,r]\)

则 \(dp_i=\max_0^r\{dp_j+1+C(j+1,i-1)\}\)

设最优决策点为 \(k\in[0,r]\) ，影响最优决策点位置的有两个方面

从 \([j+1,i-1]\) 这一段点考虑， \(j\) 越小时，就会有越多的点对被 \(B_i\) 限制，也就是说 \(j\) 越大，对于中间这一段来说越优

但是从 \(dp_j\) 的角度考虑，并不是 \(j\) 越大越好，因为可能存在一个 \(A_j\) 特别小限制了前面递增点列的选择

推论：如果前面存在一个 \(A_j>B_i\) ，那么 \(k\ge j\)

事实上应该说成 \(dp_j+C(j+1,i-1)\ge\forall d\in[0,j-1],dp_d+C(d+1,i-1)\)

综合上面两条来看，如果 \(A_j>B_i\) 意味着把它放进递增序列里绝对是优的，因为不会对前面的点产生不良限制

消除了不良限制之后，就满足最优性了

所以发现最优决策点的范围缩小到了 \([l=max\{k|A_k>B_i\},r]\)

发现决策范围内的 \(C(l+1,i-1)=C(l+2,i-1)=\cdots=C(r+1,i-1)\)

所以 \(dp_i=\max_l^r\{dp_j\}+C(r+1,i-1)=max_0^r\{dp_j\}+C(r+1,i-1)\)

所以可以直接维护一个前缀最大值，每次二分找到那个 \(r\) ，求出 \(C(r+1,i-1)\) 即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
const int N=2e5+10,K=15;
 
int n;
 
struct Node{
    int a,b;
    bool operator < (const Node __) const {
        return b<__.b;
    }
}A[N];
int h[N],hc;
int dp[N];
 
int rt[N*K],s[N*K],ls[N*K],rs[N*K],cnt;
void Add(int p,int pre,int l,int r,int x) {
    s[p]=s[pre]+1;
    if(l==r) return;
    int mid=(l+r)>>1;
    ls[p]=ls[pre],rs[p]=rs[pre];
    x<=mid?Add(ls[p]=++cnt,ls[pre],l,mid,x):Add(rs[p]=++cnt,rs[pre],mid+1,r,x);
}
 
int Que(int pl,int pr,int l,int r,int ql,int qr) {
    if(ql>qr) return 0;
    if(ql==l&&qr==r) return s[pr]-s[pl];
    int mid=(l+r)>>1;
    if(qr<=mid) return Que(ls[pl],ls[pr],l,mid,ql,qr);
    else if(ql>mid) return Que(rs[pl],rs[pr],mid+1,r,ql,qr);
    else return Que(ls[pl],ls[pr],l,mid,ql,mid)+Que(rs[pl],rs[pr],mid+1,r,mid+1,qr);
}
 
int main(){
    rep(i,1,n=rd()) {
        A[i].a=rd(),A[i].b=rd();
        h[++hc]=A[i].a,h[++hc]=A[i].b;
    }
    sort(h+1,h+hc+1),hc=unique(h+1,h+hc+1)-h-1;
    sort(A+1,A+n+1);
    rep(i,1,n) {
        A[i].a=lower_bound(h+1,h+hc+1,A[i].a)-h;
        A[i].b=lower_bound(h+1,h+hc+1,A[i].b)-h;
        Add(rt[i]=++cnt,rt[i-1],1,hc,A[i].a);
    }
    rep(i,1,n) {
        int l=1,r=i-1,res=0;
        while(l<=r) {
            int mid=(l+r)>>1;
            if(A[mid].b<A[i].a) l=mid+1,res=mid;
            else r=mid-1;
        }
        dp[i]=max(dp[i-1],dp[res]+1+Que(rt[res],rt[i-1],1,hc,A[i].b+1,hc));
    }
    int ans=dp[n];
    printf("%d\n",ans);
}

附：离线做法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include<cstdio>
#include<cctype>
#include<algorithm>
using namespace std;

#define reg register
typedef long long ll;
#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;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) f|=(IO=='-');
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=2e5+10,K=15;

int n;
struct Node{
	int a,b;
	bool operator < (const Node __)const{ return b<__.b; }
}A[N];
struct Query{
	int x,p,id,k;
	bool operator < (const Query __)const{ return p<__.p; }
}Q[N<<1];

int qc,h[N],hc,dp[N];
int L[N],Ans[N];

int s[N];
void Add(int p,int x){ while(p) s[p]+=x,p-=p&-p; }
int Que(int p){
	int res=0;
	while(p<=hc) res+=s[p],p+=p&-p;
	return res;
}

int main(){
	rep(i,1,n=rd()) {
		A[i].a=rd(),A[i].b=rd();
		h[++hc]=A[i].a,h[++hc]=A[i].b;
	}
	sort(h+1,h+hc+1),hc=unique(h+1,h+hc+1)-h-1;
	sort(A+1,A+n+1);
	rep(i,1,n) {
		A[i].a=lower_bound(h+1,h+hc+1,A[i].a)-h;
		A[i].b=lower_bound(h+1,h+hc+1,A[i].b)-h;
	}
	rep(i,1,n) {
		int l=1,r=i-1,res=0;
		while(l<=r) {
			int mid=(l+r)>>1;
			if(A[mid].b<A[i].a) l=mid+1,res=mid;
			else r=mid-1;
		}
		L[i]=res;
		if(i-1>L[i]) Q[++qc]=(Query){A[i].b+1,i-1,i,1},Q[++qc]=(Query){A[i].b+1,L[i],i,-1};
	}
	sort(Q+1,Q+qc+1);
	int p=1;
	rep(i,0,n) {
		if(i) Add(A[i].a,1);
		while(p<=qc && Q[p].p<=i) Ans[Q[p].id]+=Q[p].k*Que(Q[p].x),p++;
	}
	rep(i,1,n) dp[i]=max(dp[i-1],dp[L[i]]+1+Ans[i]);
	int ans=dp[n];
	printf("%d\n",ans);
}

Posted on Edited on In

[学军20201104CSP模拟赛] 二维码

简要题意：

对于 \(n\times m\) 的网格图，初始时全部为白色，现在 通过下面的方法染色

每次选择一个行或者列，把它全部染成黑色或者全部染成白色

求任意操作的情况下，可以得到的不同网格图的数量 \(\mod 998244353\)

\[ \ \]

判定一个染色方案是否有解的条件是：

染色完成的矩阵不包含一个子矩阵满足四个角分别为

01

10

或者

10

01

但是这样看这个条件似乎比较抽象，如果具体对于一个行上考虑，就是满足

每一行所包含的1的位置的集合之间 互为子集

显然一个方案可以任意交换行/列，不妨把按照每一行1的个数将每一行排序，设每一行有 \(a_i\) 个1，边界条件为 \(a_0=0\)

那么对于行上的1考虑排列，方案数为 \(\begin{aligned} \prod \binom{m-a_{i-1}}{a_i-a_{i-1}}\end{aligned}\) ，即从空的 \(m-a_{i-1}\) 个位置里选出多出的 \(a_i-a_{i-1}\) 个位置

而对于列之间的排列需要考虑 \(a_i\) 与 \(a_{i+1}\) 的关系，因为如果 \(a_i=a_{i+1}\) 时，必然满足这两行相同

设所有的 \(a_i\) 构成若干个相同的组，每一组包含 \(b_i(i\in[1,k])\) 个元素，则方案数显然为 \(\begin{aligned} \frac{n!}{\prod b_i!}\end{aligned}\)

而组内的 \(a_i\) 之间显然是没有 \(\begin{aligned} \sum \binom{m-a_{i-1}}{a_i-a_{i-1}}\end{aligned}\) 的贡献的，可以跳过

由此，不妨令 \(dp_{i,j}\) 表示 \(dp\) 了前 \(i\) 行，最后一行 \(a_i=j\) ，每次枚举每个组 \(b_i\) 转移

复杂度为 \(O(n^4)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#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=2010,P=998244353;
int n,m,C[N][N],dp[N][N],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 main(){
	rep(i,J[0]=1,N-1) J[i]=1ll*J[i-1]*i%P;
	I[N-1]=qpow(J[N-1]);
	drep(i,N-2,0) I[i]=1ll*I[i+1]*(i+1)%P;
	rep(i,0,N-1) rep(j,C[i][0]=1,i) C[i][j]=(C[i-1][j]+C[i-1][j-1])%P;

	scanf("%d%d",&n,&m);
	rep(i,0,n) dp[i][0]=I[i]; // 第一个块
	rep(i,1,n) rep(j,1,m) {
		rep(a,0,i-1) rep(b,0,j-1) {
			dp[i][j]=(dp[i][j]+1ll*dp[a][b]*I[i-a]%P*C[m-b][j-b])%P;
		}
	}
	int ans=0;
	rep(i,0,m) ans=(ans+dp[n][i])%P;
	ans=1ll*ans*J[n]%P;
	printf("%d\n",ans);
}

优化：

发现两维 \(dp\) 之间是相互独立的，分别是把 \(n,m\) 分组

令 \(F_{i,j}\) 表示把 \(n\) 分成了 \(i\) 个组，当前总和为 \(j\) 的方案数， \(G_{i,j}\) 表示把 \(m\) 分成 \(i\) 组，当前总和为 \(j\)

按照上面的系数转移，最后 \(O(n)\) 合并，复杂度为 \(O(n^3)\)

\[ \ \]

进一步优化：

为了方便下面的叙述，不妨先整理一下 \(a_i\) 之间转移的系数，不妨设边界 \(a_{n+1}=m\)

\(\begin{aligned} \prod \binom{m-a_{i-1}}{a_i-a_{i-1}}=\prod_{i=1}^n \frac{(m-a_{i-1})!}{(a_i-a_{i-1})!(m-a_{i})!}= \frac{m!}{\prod_{i=1}^{n+1} (a_{i}-a_{i-1})!}\end{aligned}\)

发现实际上和列之间的系数是类似的，每次枚举 \(a_i-a_{i-1}\) 即可

而实际上只有 \(k\) 个 \(b_i\) 直接相交的位置 \(a_{i}-a_{i-1}\) 有效，因此行和列实际上分别是将 \(n,m\) 分成了 \(k\) 组

观察上面的转移系数，行构成的块，首个块大小可以为 \(0\) ，而列构成的块最后一个块大小可以为 \(0\) ，所以这个并不是严格分成 \(k\) 组，下面会讨论这个问题

我们计算答案的复杂度消耗在计算分成若干块的方案，而实际上，把 \(n\) 分成 \(k\) 块的方案数就是 \(\begin{Bmatrix} n\\ k\end{Bmatrix}\cdot k!\)

用 \(n^2\) 递推第二类斯特林数的方法即可计算

对于并不是严格分成 \(k\) 组的问题，可以考虑把开头/结尾那一个大小为 \(0\) 的块删掉，即同时还要考虑 \(\begin{Bmatrix}n \\ k-1\end{Bmatrix}(k-1)!\)

最后再枚举 \(k\) ，复杂度为 \(O(n^2)\)

更优化的就是 \(\text{NTT}\) 计算斯特林数，带入通项公式

\(\begin{aligned} \begin{Bmatrix}n\\ m\end{Bmatrix}m!=\sum_{i=0}^m i^n(-1)^{m-i}\binom{m}{i} \end{aligned}\)

显然把组合数拆开 \(\text{NTT}\) 即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#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)

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

int n,m,ans;
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 I[N],J[N];

int rev[N];
int Init(int n){ 
	int R=2,c=0;
	while(R<=n) R<<=1,c++;
	rep(i,0,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
	return R;
}

void NTT(int n,int *a,int f){
	static int e[N>>1];
	rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=e[0]=1;i<n;i<<=1) {
		ll t=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
		for(int j=i-2;j>=0;j-=2) e[j+1]=(e[j]=e[j>>1])*t%P;
		for(int l=0;l<n;l+=i*2) {
			for(int j=l;j<l+i;++j) {
				int 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) {
		ll base=qpow(n);
		rep(i,0,n-1) a[i]=a[i]*base%P;
	}
}

int A[N],B[N],C[N];

int main(){
	scanf("%d%d",&n,&m);
	if(n<m) swap(n,m);
	rep(i,J[0]=1,N-1) J[i]=1ll*J[i-1]*i%P;
	I[N-1]=qpow(J[N-1]);
	drep(i,N-2,0) I[i]=1ll*I[i+1]*(i+1)%P;

	int R=Init(n+n+2);
	rep(i,0,n) A[i]=qpow(i,n)*I[i]%P;
	rep(i,0,n) B[i]=(i&1)?P-I[i]:I[i];
	NTT(R,A,1),NTT(R,B,1);
	rep(i,0,R-1) A[i]=1ll*A[i]*B[i]%P;
	NTT(R,A,-1);
	rep(i,n+1,R) A[i]=0;

	rep(i,0,m) C[i]=qpow(i,m)*I[i]%P;
	NTT(R,C,1);
	rep(i,0,R-1) C[i]=1ll*C[i]*B[i]%P;
	NTT(R,C,-1);
	rep(i,m+1,R) C[i]=0;

	int ans=0;
	rep(i,0,min(n,m)) ans=(ans+1ll*(1ll*A[i]*J[i]%P+1ll*A[i+1]*J[i+1]%P)*(1ll*C[i]%P*J[i]%P+1ll*C[i+1]*J[i+1]%P))%P;
	printf("%d\n",ans);
}
0%