zoukankan      html  css  js  c++  java
  • 【XSY3042】石像 拓扑排序 状压DP 洲阁筛

    题目大意

      有 (n) 个整数 (a_1,a_2,ldots,a_n),每个数的范围是 ([1,m])。还有 (k) 个限制,每个限制 (x_i,y_i) 表示 (a_{x_i}leq a_{y_i})

      问有多少种不同的情况,以及所有情况中 ({sigma_0(gcd(a_1,a_2,ldots,a_n))}^3) 的和。

      (nleq 20,mleq {10}^{10})

    题解

      记 (f(x)) 为当 (m=x) 时第一问的答案。

      记 (g(x)) 为当 (a_1,a_2,ldots,a_n)(x) 种不同的值时第一问的答案。

      (g(x)) 要怎么求呢?

      tarjan 求强连通分量,然后缩点,状压DP。

      直接枚举子集是 (O(n3^n)) 的。

      因为缩点后是一个 DAG,所以按拓扑序来转移就好了。

      设 (h_{i,j,S}) 为 当前考虑的值为 (i),正在考虑 (a_j) 这个数,已经确定数的集合为 (S) 的方案数。

      每次决定 (a_j) 这个数的值是不是 (i) 就可以了。

      时间复杂度:(O(n^22^n))

      (f(x)) 要怎么求呢?

      注意到 (a_1,a_2,ldots,a_n) 最多只有 (n) 种不同的值,所以可以用 (g(1),g(2),ldots,g(n)) 计算出 (f(x))

    [f(x)=sum_{i=1}^ninom{x}{i}g(i) ]

      求组合数时要记下 (2) 的次数和非 (2) 的倍数的乘积两部分。

      或者你求出前 (O(n)) 项然后拉格朗日插值也是可以的。

      这样就求出第一问的答案了。

      时间复杂度:(O(n))

      那第二问要怎么求呢?

    [egin{align} &sum_{a_1=1}^msum_{a_2=1}^mcdotssum_{a_n=1}^msigma_0({gcd(a_1,a_2,ldots,a_n)}^3)\ =&sum_{i=1}^msigma_0(i^3)^3sum_{a_1=1}^msum_{a_2=1}^mcdotssum_{a_n=1}^m[gcd(a_1,a_2,ldots,a_n)=i]\ =&sum_{i=1}^m(sum_{jmid i}sigma_0(j^3)^3mu(frac{i}{j}))sum_{a_1=1}^msum_{a_2=1}^mcdotssum_{a_n=1}^m[imidgcd(a_1,a_2,ldots,a_n)]\ =&sum_{i=1}^m(sum_{jmid i}sigma_0(j^3)^3mu(frac{i}{j}))f(lfloorfrac{m}{i} floor) end{align} ]

      那么前面那部分要怎么求呢?

      容易发现 (h(n)=sum_{dmid n}sigma_0(d^3)^3mu(frac{n}{d})) 是一个积性函数,满足

    [egin{align} h(1)&=1\ h(p)&=1 imes (-1)+4^3 imes 1=63\ h(p^c)&={(3(c-1)+1)}^3(-1)+{(3c+1)}^3(1)\ &=(27c^3+27c^2+9c+1)-(27c^2-54c^2+36c-8)\ &=81c^2-27c+9 end{align} ]

      所以可以大力洲阁筛/min_25筛求。

      可能需要卡一下常。

      时间复杂度:(O(n^22^n+nsqrt m+frac{m^frac{3}{4}}{log m}))

    代码

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cstdlib>
    #include<ctime>
    #include<utility>
    #include<cmath>
    #include<functional>
    using namespace std;
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int,int> pii;
    typedef pair<ll,ll> pll;
    void sort(int &a,int &b)
    {
    	if(a>b)
    		swap(a,b);
    }
    void open(const char *s)
    {
    #ifndef ONLINE_JUDGE
    	char str[100];
    	sprintf(str,"%s.in",s);
    	freopen(str,"r",stdin);
    	sprintf(str,"%s.out",s);
    	freopen(str,"w",stdout);
    #endif
    }
    int rd()
    {
    	int s=0,c,b=0;
    	while(((c=getchar())<'0'||c>'9')&&c!='-');
    	if(c=='-')
    	{
    		c=getchar();
    		b=1;
    	}
    	do
    	{
    		s=s*10+c-'0';
    	}
    	while((c=getchar())>='0'&&c<='9');
    	return b?-s:s;
    }
    void put(int x)
    {
    	if(!x)
    	{
    		putchar('0');
    		return;
    	}
    	static int c[20];
    	int t=0;
    	while(x)
    	{
    		c[++t]=x%10;
    		x/=10;
    	}
    	while(t)
    		putchar(c[t--]+'0');
    }
    int upmin(int &a,int b)
    {
    	if(b<a)
    	{
    		a=b;
    		return 1;
    	}
    	return 0;
    }
    int upmax(int &a,int b)
    {
    	if(b>a)
    	{
    		a=b;
    		return 1;
    	}
    	return 0;
    }
    
    int fp(int a,ll b)
    {
    	int s=1;
    	for(;b;b>>=1,a*=a)
    		if(b&1)
    			s*=a;
    	return s;
    }
    
    namespace gao1
    {
    	const int N=30;
    	int n,m;
    	vector<int> g[N],g2[N];
    	int lx[N*N],ly[N*N];
    	int dfn[N],low[N],b[N],st[N],top,c[N],ti,ccnt;
    	int f[1<<20],a[N];
    	int s[N];
    	int d[N][N];
    	int inv[N],pw[N];
    	void dfs(int x)
    	{
    		dfn[x]=low[x]=++ti;
    		st[++top]=x;
    		b[x]=1;
    		for(auto v:g[x])
    			if(b[v]!=2)
    			{
    				if(!b[v])
    					dfs(v);
    				low[x]=min(low[x],low[v]);
    			}
    		if(low[x]>=dfn[x])
    		{
    			ccnt++;
    			int v;
    			do
    			{
    				v=st[top--];
    				b[v]=2;
    				c[v]=ccnt;
    			}
    			while(v!=x);
    		}
    	}
    	void add(int &a,int b)
    	{
    		a+=b;
    	}
    	void gao()
    	{
    		scanf("%d",&m);
    		for(int i=1;i<=m;i++)
    		{
    			scanf("%d%d",&lx[i],&ly[i]);
    			g[lx[i]].push_back(ly[i]);
    		}
    		for(int i=1;i<=n;i++)
    			if(!b[i])
    				dfs(i);
    		for(int i=1;i<=m;i++)
    			if(c[lx[i]]!=c[ly[i]])
    			{
    				g2[c[lx[i]]].push_back(c[ly[i]]);
    				a[c[lx[i]]]|=1<<(c[ly[i]]-1);
    			}
    		n=ccnt;
    		f[0]=1;
    		for(int i=1;i<=n;i++)
    		{
    			for(int j=1;j<=n;j++)
    				for(int k=0;k<1<<n;k++)
    					if(!((k>>(j-1))&1)&&(k&a[j])==a[j])
    						add(f[k|(1<<(j-1))],f[k]);
    			s[i]=f[(1<<n)-1];
    		}
    		for(int i=0;i<=n;i++)
    		{
    			d[i][0]=1;
    			for(int j=1;j<=i;j++)
    				d[i][j]=d[i-1][j-1]+d[i-1][j];
    		}
    		for(int i=n;i>=1;i--)
    		{
    			for(int j=i-1;j>=1;j--)
    				s[i]+=((i-j)&1?-1:1)*d[i][j]*s[j];
    		}
    		for(int i=1;i<=n;i++)
    		{
    			int x=i;
    			while(!(x&1))
    			{
    				pw[i]++;
    				x>>=1;
    			}
    			inv[i]=fp(x,(1ll<<31)-1);
    		}
    	}
    	int query(ll x)
    	{
    		int res=0;
    		int prod=1;
    		int num=0;
    		for(int i=1;i<=n&&i<=x;i++)
    		{
    			ll y=x-i+1;
    			while(num<=20&&!(y&1))
    			{
    				y>>=1;
    				num++;
    			}
    			num-=pw[i];
    			prod*=inv[i];
    			prod*=y;
    			res+=(prod<<num)*s[i];
    		}
    		return res;
    	}
    }
    
    namespace gao2
    {
    	const int M=100010;
    	const int m=100000;
    	ll n;
    	int pri[M],b[M],cnt;
    	int f1[M],f2[M],s1[M],s2[M];
    	void init()
    	{
    		for(int i=2;i<=m;i++)
    		{
    			if(!b[i])
    				pri[++cnt]=i;
    			for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
    			{
    				b[i*pri[j]]=1;
    				if(i%pri[j]==0)
    					break;
    			}
    		}
    	}
    	void gao()
    	{
    		init();
    		int mx=n/(m+1);
    		for(int i=2;i<=m;i++)
    			f1[i]=i-1;
    		for(int i=1;i<=mx;i++)
    			f2[i]=n/i-1;
    		for(int i=1;i<=cnt;i++)
    		{
    			int x=pri[i];
    			int y=f1[x-1];
    			int mx1=min(n/x/(m+1),n/x/x);
    			int mx2=min((ll)mx,n/x/x);
    			ll p2=(ll)pri[i]*pri[i];
    			ll o=n/x;
    			for(int j=1;j<=mx1;j++)
    				f2[j]-=f2[j*x]-y;
    			for(int j=mx1+1;j<=mx2;j++)
    				f2[j]-=f1[o/j]-y;
    			for(int j=m;j>=p2;j--)
    				f1[j]-=f1[j/x]-y;
    		}
    		for(int i=2;i<=m;i++)
    		{
    			f1[i]*=63;
    			s1[i]=f1[i];
    		}
    		for(int i=1;i<=mx;i++)
    		{
    			f2[i]*=63;
    			s2[i]=f2[i];
    		}
    		for(int i=cnt;i>=1;i--)
    		{
    			int x=pri[i];
    			int mx1=min((ll)mx,n/x/x);
    			for(int j=1;j<=mx1;j++)
    			{
    				s2[j]-=63;
    				ll x2=n/j/x;
    				for(int k=1;x2;x2/=x,k++)
    					s2[j]+=((x2<=m?s1[x2]-f1[min((int)x2,x)]:s2[n/x2]-f1[x])+1)*(81*k*k-27*k+9);
    			}
    			ll mn=(ll)pri[i]*pri[i];
    			for(int j=m;j>=mn;j--)
    			{
    				s1[j]-=63;
    				int x2=j/x;
    				for(int k=1;x2;x2/=x,k++)
    					s1[j]+=(s1[x2]-f1[min(x2,x)]+1)*(81*k*k-27*k+9);
    			}
    		}
    	}
    	int query(ll x)
    	{
    		if(x<=0)
    			return 0;
    		return (x<=m?s1[x]:s2[n/x])+1;
    	}
    }
    
    int main()
    {
    	open("statue");
    	scanf("%d%lld",&gao1::n,&gao2::n);
    	gao1::gao();
    	gao2::gao();
    	int ans=0;
    	ll n=gao2::n;
    	for(ll i=1,j;i<=n;i=j+1)
    	{
    		j=n/(n/i);
    		ans+=(gao2::query(j)-gao2::query(i-1))*gao1::query(n/i);
    	}
    	printf("%u %u
    ",gao1::query(n),ans);
    	return 0;
    }
    
  • 相关阅读:
    c# 三种取整方法 向上取整 向下取整 四舍五入
    Lambda表达式对DataRow处理
    Dapper数据库字段和model属性映射
    union limit
    北邮五十题
    搜索____深搜 学易错点
    动态规划____有重叠子问题的搜索,都可以转为记忆化搜索
    64位 __int 与 long long写法
    做做 卡特兰数 与 卡米歇尔数
    vector 有点麻烦啊 能简单点么?
  • 原文地址:https://www.cnblogs.com/ywwyww/p/9229200.html
Copyright © 2011-2022 走看看