魔术球问题

题解 魔术球问题

前言

这道题提交之后看到满屏的ac，顿时流下纵横老泪，喜极而泣，不能自已，欣喜若狂，语无伦次sdljfaldsjflgjaccnznfdjajqioe

这一道题搞了快一个下午没搞出来，疯狂wa,tle,re，最后却能用36ms收尾，真是感慨万千

正在前往 洛谷

题目誊录如下:

假设有n根柱子，现要按下述规则在这n根柱子中依次放入编号为1，2，3，...的球。

（1）每次只能在某根柱子的最上面放球。

（2）在同一根柱子中，任何2个相邻球的编号之和为完全平方数。

试设计一个算法，计算出在n根柱子上最多能放多少个球。例如，在4 根柱子上最多可放11 个球。

«编程任务：

对于给定的n，计算在n根柱子上最多能放多少个球。

step 1 一脸懵逼

由于上一道题刚刚做了“最小路径覆盖问题”，一看到题就有了一点点思路。

因为每一个球能放在另外一个合法的球的可能性并不是很高。（暴力跑了一下，好像球的编号到了1000左右，也只有约4，5个球可以合法地放在其上）那么就分点，每一种可以放的组合就对应建一条边，然后跑一下最小路径覆盖的模型即可。

然而我再稍微想想就有点懵逼。题目只给了柱子数，并没有给出确切的球数？？那怎么建图？怎么二分图最大匹配？

emm，本着刷题（懒惰）的精神，为了追求速度更快的结束去摸鱼我点开了题解。

嗯，果然还是太菜了。考虑到当球数增加，最少需要的柱子数必定是不下降序列，也就是说满足单调性。那么只需要不断加球，直到某一次加球后，最少需要的柱子数大于了n，就停止，回到上一个时间点就ok嘛

step 2 万念俱灰 12005ms 0pts ~ 23ms+7re 36pts

然后，按照“最小路径覆盖问题”的做法，老老实实还原了二分图的匈牙利。

但是？wtf？第一代评测记录 12005ms 0pts

emm，咋回事妮？

呼唤了某l姓oj，偷了他的数据之后，跑了一遍，也不知所云。

再一次翻看，发现了蹊跷

这个判完全平方数的for循环好像也不太相同，难道是说的我的两遍sqrt()导致精度丢失的太严重？事实证明这也确实是一方面

于是再原来的基础上，修改了一下判断完全平方数的for循环。

第二代评测记录23ms+7re 36pts

至此，二分图匈牙利原版算法完整代码：

// 二分图匈牙利原版
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <ctime>
using namespace std;

const int MAX=2e4+5;
int n,m;
int ecnt,edge[MAX],head[MAX],nxt[MAX];
int mx[MAX],my[MAX]; bool vis[MAX];

clock_t s,e;

void insert(int,int,int);
int magnolia();
int point(int);
void print(int);

int main(){
	freopen("test.in","r",stdin);
	// freopen("test.out","w",stdout);
	s=clock();
	cin>>n;

	m=0; ecnt=0;
	while(true){
		m++; int cnt=0;
		for(int i=sqrt(m)+1;i*i<(m<<1);++i){
			int cur=i*i-m;
			if(cur<=0||cur>m) continue;
			insert(i*i-m,m,++ecnt);
			cnt++;
		}

		if(m-magnolia()>n) break;
		// if(m==1567){
			// break;
		// }
	}

	m--;
	// for(int i=1;i<=m;++i) if(edge[head[i]]==m+1) head[i]=nxt[head[i]];
	// magnolia();
	// cout<<m<<endl;
	// for(int i=1;i<=m;++i){
	// 	if(!my[i]) print(i);
	// }
	e=clock();
	printf("%lf",(double)(e-s)/CLOCKS_PER_SEC);
	return 0;
}


void insert(int from,int to,int id){
	nxt[id]=head[from]; head[from]=id; edge[id]=to;
}

int magnolia(){
	int ans=0;
	for(int i=1;i<=m;++i) mx[i]=my[i]=0;
	for(int i=1;i<=m;++i){
		if(mx[i]) continue;
		for(int j=1;j<=m;++j) vis[j]=false;
		ans+=point(i);
	}
	return ans;
}

int point(int u){
	for(int i=head[u];i;i=nxt[i]){
		int v=edge[i]; if(vis[v]) continue;
		vis[v]=true;
		if(!my[v]||point(my[v])){
			mx[u]=v;
			my[v]=u;
			return 1;
		}
	}
	return 0;
}

void print(int u){
	printf("%d ",u);
	if(mx[u]) print(mx[u]);
	else printf("\n");
}

step 3 一线希望 1474ms 55pts ~ 1076ms 100pts

然后妮，接着摸题解，又发现了一条诡异的特征

题解里只要讲的只要不是贪心算法，用的都是dinic，竟然没看到用匈牙利的。后来资料和实验双重证明dinic确实比匈牙利更加优秀

也难怪“最小路径覆盖”这道题的范围是\(<=150\)，但是这道题的范围最大为\(<=1567\),是不凉，熟可凉？

于是，重构代码！！

由于dinic和匈牙利完全不是一回事儿，需要完全推翻重写。不过好在匈牙利留了下来，这在将来起到了重大作用。

二分图里写dinic可比写magnolia难写多了。过程中不断碰到大大小小的bug，一路打下来，硬是把dinic打了下来

第三代评测记录 1474ms 55pts

emm？为什么还是没过？接着是漫长的找bug...

期间将第二个数据下了下来，运行结果，人脑spj。一度怀疑有着1.6k ac 的spj是不是有问题...

最后还是发现了问题所在。原来是为了最一次的dinic，做出的初始化中，出现了小bug。

emm，好了，到此，终于残喘ac

第四代评测记录 1076ms 100pts

到此为止，网络流dinic算法完整代码：

//网络流 dinic算法

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cmath>
#include <ctime>
using namespace std;
#define x(x) ((x+1)<<1)
#define y(x) (((x+1)<<1)|1)

const int MAX=5e4+5,INF=0x3f3f3f3f;

int n,m,src,tar;
int ecnt,edge[MAX],head[MAX],nxt[MAX],wei[MAX];
int dep[MAX];
int cx[MAX],cy[MAX];
clock_t s,e;

void insert(int,int,int,int);
int dinic();
bool bfs();
int dfs(int,int);
void print(int);

int main(){
	freopen("test.in","r",stdin);
	freopen("test.out","w",stdout);
	cin>>n;

	src=1; tar=2; ecnt=1;
	while(true){
		m++;
		insert(src,x(m),1,++ecnt); insert(x(m),src,0,++ecnt);
		insert(y(m),tar,1,++ecnt); insert(tar,y(m),0,++ecnt);
		for(int i=sqrt(m)+1;i*i<(m<<1);++i){
			int tmp=i*i-m;
			insert(x(tmp),y(m),1,++ecnt); insert(y(m),x(tmp),0,++ecnt);
		}

		for(int i=2;i<=ecnt;++i){
			if(i%2==0) wei[i]=1;
			else wei[i]=0;
		}

		if(m-dinic()>n) break;
	}
	m--; cout<<m<<endl;

	for(int i=1;i<=m;++i){
		int u=x(i);
		if(edge[head[u]]==y(m+1)) head[u]=nxt[head[u]];
	}
	for(int i=2;i<=ecnt;++i){
		if(i%2==0) wei[i]=1;
		else wei[i]=0;
	}

	dinic();
	for(int i=1;i<=m;++i){
		int u=x(i);
		for(int j=head[u];j;j=nxt[j]){
			int v=edge[j]; if(wei[j]||v==src) continue;
			int to=(((v-1)>>1)-1);
			cx[i]=to; cy[to]=i;
			break;
		}
	}

	for(int i=1;i<=m;++i){
		if(!cy[i]) print(i);
	}
	return 0;
}

void insert(int from,int to,int w,int id){
	nxt[id]=head[from]; head[from]=id; edge[id]=to; wei[id]=w;
}

int dinic(){
	int ans=0;
	while(bfs()) ans+=dfs(src,INF);
	return ans;
}

bool bfs(){
	int line[MAX],l=0,r=1; line[r]=src;
	for(int i=1;i<=y(m);++i) dep[i]=0;
	dep[src]=1;
	while(l<r){
		int u=line[++l];
		for(int i=head[u];i;i=nxt[i]){
			int v=edge[i];
			if(dep[v]||!wei[i]) continue;
			dep[v]=dep[u]+1;
			line[++r]=v;
		}
	}
	if(dep[tar]) return true;
	else return false;
}

int dfs(int u,int flow){
	int ans=0;
	if(u==tar||!flow) return flow;
	for(int i=head[u];i;i=nxt[i]){
		int v=edge[i]; if(dep[v]!=dep[u]+1) continue;
		int add=dfs(v,min(flow,wei[i]));
		if(add){
			wei[i]-=add;
			wei[i|1]+=add;
			flow-=add;
			ans+=add;
		}
	}
	return ans;
}

void print(int u){
	printf("%d ",u);
	if(cx[u]) print(cx[u]);
	else printf("\n");
}	

step 4 豁然开朗 38ms 100pts ~ 36ms 100pts

打完之后，感觉身体被掏空。做了一下午的题总算是踩线做了出来。但是看了最优解之后真是感觉心塞——我的程序所需时间相较最优解的时间之差差不多就是题目的时间限制。

感觉被生活欺骗。

然后，突然就有一个新的思路！！

想一想匈牙利算法。它的思路是挨个为每一个点找增广路。假如我把\(u\)->\(v\)的边反着建成\(v\)->\(u\)，那么就可以保证\(v<u\),也就是说在为前面的点寻找增广路的时候，和比当前点大的点没有任何联系！

那么，只需要每一次建点，添边的时候，只是跑一边增广路即可，而不是从头到尾跑一边magnolia！效率大大的提升！

第五代评测记录 38ms 100pts

然后妮，吸氧跑了一遍，并没有明显优化，可能我的代码不是好氧性代码吧？ （；´д｀）ゞ

第六代评测记录 36ms 100pts

此致，二分图magnolia优化版本：

//二分图 匈牙利的改进版
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <ctime>
using namespace std;

const int MAX=2e4+5;
int n,m;
int ecnt,edge[MAX],head[MAX],nxt[MAX];
int mx[MAX],my[MAX]; bool vis[MAX];

clock_t s,e;

void insert(int,int,int);
int magnolia();
int point(int);
void print(int);

int main(){
	freopen("test.in","r",stdin);
	// freopen("test.out","w",stdout);
	s=clock();
	cin>>n;

	int mi=0;
	while(true){
		m++;
		for(int i=sqrt(m)+1;i*i<(m<<1);++i){
			int cur=i*i-m;
			if(cur<=0||cur>m) continue;
			insert(m,i*i-m,++ecnt);
			mi+=point(m);
		}

		if(m-mi>n) break;
	}

	m--; cout<<m<<endl;

	for(int i=1;i<=m;++i){
		if(edge[head[i]]==m+1) head[i]=nxt[head[i]];
	}
	magnolia();
	for(int i=1;i<=m;++i){
		if(!mx[i]) print(i);
	}

	// e=clock();
	// printf("%lf",(double)(e-s)/CLOCKS_PER_SEC);
	return 0;
}


void insert(int from,int to,int id){
	nxt[id]=head[from]; head[from]=id; edge[id]=to;
}

int magnolia(){
	int ans=0;
	for(int i=1;i<=m;++i) mx[i]=my[i]=0;
	for(int i=1;i<=m;++i){
		if(mx[i]) continue;
		for(int j=1;j<=m;++j) vis[j]=false;
		ans+=point(i);
	}
	return ans;
}

int point(int u){
	for(int i=head[u];i;i=nxt[i]){
		int v=edge[i]; if(vis[v]) continue;
		vis[v]=true;
		if(!my[v]||point(my[v])){
			mx[u]=v;
			my[v]=u;
			return 1;
		}
	}
	return 0;
}

void print(int u){
	printf("%d ",u);
	if(my[u]) print(my[u]);
	else printf("\n");
}

后记

这道题足足做了我半天，整整半天。但是做完之后真是成就感满满。“最小路径覆盖”将会像“a+b problem”一样刻在我的脑子里。

愿匈牙利如 \(magnolia\) 般有高尚的灵魂,dinic什么的可以去死吧。(。・∀・)ノ