NanoApe's Blog

既是咸鱼又是辣鸡

烷烃同分异构体个数的计算——BZOJ 4271:Chemistry

NanoApe posted @ 2016年1月08日 20:19 in 蒟蒻整天被神题虐哭 , 7467 阅读

源于一节化学课的思想跳跃


化学课终于教到有机了,第一次知道什么是甲烷

看到烷烃的结构,尼玛这不就是棵树嘛

化学老师也是闲,叫我们来推推戊烷己烷的同分异构体个数

这对一个有OI基础的人来说基本不难

然后就顺势想到了如何编程求烷烃的同分异构体个数

心想着:这毕竟是一棵树啊,树的话OI就能考。。。

激动万分

然而想了想发现自己连口胡算法都不会。。。

为自己的智商感动得泪流满面

 

下课上网找题解,抱着不大的希望在BZOJ搜索了下Chemistry,就找到原题了

泪流满面,然而才8人A,其中大部分都是打表的。。。

找到ACM原题也没能找到题解。。。我擦勒

问了下Claris,知道了貌似要枚举树重心去DP一发。。。听都没听过的新思路

维基搜了下,没找到公式,OEIS的列表倒是有一点用(别跟我说公式其实就藏在OEIS下面一堆英文中QwQ

还找到了范浩强神犇的Blog,解法也是与重心有关,然而题解的网站挂了。。。

我擦勒什么题解都没有,而且听到算法梗概就觉得已经超出了我的思考能力了

泪流满面

 

然后隔没多久,接触到了有根树无根树的计数方法

我擦勒无根树计数不就是和求烷烃的数量啊

想了想,度数限制不大于4。。。感动

接下来就是一系列不服输的行为了

照着公式乱啃找计数和去重原理,自己琢磨出度数限制的方案

两节理综荒废,最后推出了一个什么鬼的递推公式。。。

上机实验,发现各种错误,在草稿纸上改公式啊改公式

我还把35种壬烷都画了出来。。。

然后弄出了个$O(n^3)$的算法,顺利过掉ACM的$N\le50$规模

BZOJ的$N\le500$由于高精度和算法复杂度太大,最慢9s

睡觉的时候才发现自己硬是把$O(n^2)$的算法写成$O(n^3)$

智商感人

 

 

思路就是用背包DP计算满足度数限制(每个点的度数不大于4)的有根树的数量

然后去重计算满足度数限制的无根树的数量

原理类似普通有根树和无根树的计数原理,建议结合看

 

#include <cstdio>
#include <algorithm>

#define rep(i, l, r) for(int i=l; i<=r; i++)
#define down(i, l, r) for(int i=l; i>=r; i--)
#define clr(x, c) memset(x, c, sizeof(x))
#define travel(x) for(edge *p=fir[x]; p; p=p->n)
#define pb push_back
#define ll long long
#define ull unsigned long long

#define maxn 501
#define maxl 29
#define mx 100000000
#define inf 0x7fffffff

using namespace std;
inline int read()
{
	int x=0, f=1; char ch=getchar();
	while (ch<'0' || ch>'9') {if (ch=='-') f=-1; ch=getchar();}
	while ('0'<=ch && ch<='9') x=x*10+ch-'0', ch=getchar();
	return x*f;
}

struct node{int v[maxl], l;} f[maxn], k[maxn], h[maxn], g[5][maxn], sum, I, Nul, m[5];
int n, T; ll tmpc[maxl];

inline void Add(node &a, node b)
{
	rep(i, 0, b.l-1) a.v[i]+=b.v[i];
	a.l=max(a.l,b.l);
	rep(i, 0, a.l-1) if (a.v[i]>=mx) a.v[i+1]+=a.v[i]/mx, a.v[i]%=mx;
	while (a.v[a.l]!=0) 
	{
		if (a.v[a.l]>=mx) a.v[a.l+1]+=a.v[a.l]/mx, a.v[a.l]%=mx;
		a.l++;
	}
}

inline void Sub(node &a, node b)
{
	rep(i, 0, a.l-1) a.v[i]-=b.v[i];
	rep(i, 0, a.l-1) while (a.v[i]<0) a.v[i+1]--, a.v[i]+=mx;
	while (a.v[a.l-1]==0 && a.l>1) a.l--;
}

inline void Add(node &a)
{
	a.v[0]++;
	rep(i, 0, a.l-1) if (a.v[i]>=mx) a.v[i+1]+=a.v[i]/mx, a.v[i]%=mx; else return;
	while (a.v[a.l]!=0) 
	{
		if (a.v[a.l]>=mx) a.v[a.l+1]+=a.v[a.l]/mx, a.v[a.l]%=mx;
		a.l++;
	}
}

inline node Sub(node a)
{
	a.v[0]--;
	rep(i, 0, a.l-1) if (a.v[i]<0) a.v[i+1]--, a.v[i]+=mx; else return a;
	while (a.v[a.l-1]==0 && a.l>1) a.l--;
	return a;
}

inline node Div(node a, int x)
{
	int tmp=0;
	down(i, a.l-1, 0) tmp=tmp*mx+a.v[i], a.v[i]=tmp/x, tmp%=x;
	while (a.v[a.l-1]==0 && a.l>1) a.l--;
	return a;
}

node operator * (node a, node b)
{
	node c=Nul; c.l=a.l+b.l-1;
	//clr(tmpc,0); 
	rep(i, 0, c.l) tmpc[i]=0;
	rep(i, 0, a.l-1) rep(j, 0, b.l-1) tmpc[i+j]+=1LL*a.v[i]*b.v[j];
	
	rep(i, 0, c.l-1) if (tmpc[i]>=mx) tmpc[i+1]+=tmpc[i]/mx, c.v[i]=tmpc[i]%mx; else c.v[i]=tmpc[i];
	while (tmpc[c.l]!=0) 
	{
		if (tmpc[c.l]>=mx) tmpc[c.l+1]=tmpc[c.l]/mx, c.v[c.l]=tmpc[c.l]%mx; else tmpc[c.l+1]=0, c.v[c.l]=tmpc[c.l];
		c.l++;
	}
	return c;
}

node C(node x, int y)
{
	node now=I;
	rep(i, 1, y) now=now*x, Add(x);
	rep(i, 1, y) now=Div(now,i);
	return now;
}

inline void Put(node &a)
{
	printf("%d", a.v[a.l-1]);
	down(i, a.l-2, 0) printf("%08d", a.v[i]);
}

int main()
{
	n=read();
	I.l=1; I.v[0]=1;
	f[1]=k[1]=h[1]=I; rep(o, 0, 4) g[o][o]=I;
	rep(a, 2, n)
	{
		//printf("%d %d\n", a, T);
		sum=Nul;
		rep(o, 1, 3) Add(sum,g[o][a-1]);
		f[a]=sum;
		
		Add(sum,g[4][a-1]);
		rep(c, 1, 4) if (a*c<=n) m[c]=C(f[a],c);
		
		down(o, 4, 1) rep(c, 1, o) rep(i, a*c+o-c, n-1) if (g[o-c][i-a*c].l!=0)
		{
			node &A=g[o][i], B=g[o-c][i-a*c]*m[c];
			rep(i, 0, B.l-1) A.v[i]+=B.v[i];
			A.l=max(A.l,B.l);
			rep(i, 0, A.l-1) if (A.v[i]>=mx) A.v[i+1]+=A.v[i]/mx, A.v[i]%=mx;
			while (A.v[A.l]!=0) 
			{
				if (A.v[A.l]>=mx) A.v[A.l+1]+=A.v[A.l]/mx, A.v[A.l]%=mx;
				A.l++;
			}
		}
		k[a]=sum;
		
		rep(i, 1, (a-1)/2) if (i!=a-i) Sub(sum,f[i]*f[a-i]);
		if (a%2==0) Sub(sum,Div(f[a/2]*Sub(f[a/2]),2));
		h[a]=sum;
	}
	Put(h[n]);
    return 0;
}

 

有老司机知道枚举重心DP的做法么。。。现在依旧不知道QwQ

 

update 16.01.10:

后来在书中翻到了一句话:

假如两棵无根树同构,那么以各自的重心作为根的有根树也同构

也许这就是正解了。。。

有根树背包DP的时候只要使得子树大小小于$\frac n2$就行了

748ms AC

 

一些链接:

OEIS序列

范浩强Blog——烷烃同分异构体个数的计算

ACM原题

Avatar_small
zhouzixuan 说:
Feb 10, 2016 03:43:13 PM

大神,求问在哪里学习的有根数计数方法
我再网上看到一篇,貌似运用了高数知识TAT

Avatar_small
NanoApe 说:
Feb 10, 2016 05:55:16 PM

就背包一发就可以了啊。。。
貌似记得是自己YY出来的。。。
嗯记得有根树计数是有公式的,你要我可以找找

Avatar_small
zhouzixuan 说:
Feb 10, 2016 08:55:12 PM

%%%
能不能大致说说做法捏?

Avatar_small
NanoApe 说:
Feb 12, 2016 10:11:36 AM

有根树的话就背包f(i,j)表示结点数为i且森林数为j的方案数是多少(不计顺序)
然后我们要求结点数为a的有根树数目,就从f(a-1,k)得到,k任取

WerKeyTom_FTD 说:
Mar 04, 2016 08:21:32 PM

%%%,时隔9个月大神写的东西终于让我再去切了这道题。
“假如两棵无根树同构,那么以各自的重心作为根的有根树也同构”,受益匪浅。

Avatar_small
NanoApe 说:
Mar 05, 2016 02:22:53 PM

啊我貌似是从哪本书看来的这句话(撕烤)

BPM136 说:
Mar 21, 2016 08:41:57 AM

感觉之前在某GD最坑爹学校集训写的某道名字是化学的题一模一样= =做法差不多= =
子树大小不超过1/2不应该是最小表示法的问题吗,否则就可以旋转变得相同
然后树的重心也是一个DP
因为树的重心的性质,可以logn枚举
然后最小表示法表示子树
这样就是n^3logn了= =
并不会n^3做法


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter