zoukankan      html  css  js  c++  java
  • [BZOJ3994][SDOI2015]约数个数和

    BZOJ
    Luogu
    题意:
    给定n,m,求(sum_{i=1}^{n}sum_{j=1}^{m}d(ij)),其中(d(x))表示x的约数个数。多组数据,n,m<=50000,T<=50000

    sol

    首先我们大胆猜想,

    [d(ij)=sum_{u|i}sum_{v|j}[gcd(u,v)==1] ]

    (证明晚上再补。。。)
    好然后我们看,我们已知

    [ans=sum_{i=1}^{n}sum_{j=1}^{m}sum_{u|i}sum_{v|j}[gcd(u,v)==1] ]

    显然这四个(sum)很不好搞,所以我们考虑枚举(u,v),计算每一对((u,v))的贡献。

    [ans=sum_{u=1}^{n}sum_{v=1}^{m}[gcd(u,v)==1]lfloor frac nu floorlfloor frac mv floor ]

    写出这个形式那就好办了,

    [f(d)=sum_{u=1}^{n}sum_{v=1}^{m}[gcd(u,v)==d]lfloor frac nu floorlfloor frac mv floor ]

    [F(d)=sum_{u=1}^{n}sum_{v=1}^{m}[d|gcd(u,v)]lfloor frac nu floorlfloor frac mv floor ]

    (F(d))的表达式中显然(u)(v)都是(d)的倍数,所以我们可以令(u=id,v=jd)然后

    [F(d)=sum_{i=1}^{n/d}sum_{j=1}^{m/d}lfloor frac {n}{id} floorlfloor frac {m}{jd} floor=sum_{i=1}^{n/d}lfloor frac {n/d}{i} floor * sum_{j=1}^{m/d}lfloor frac {m/d}{j} floor ]

    注意上面的那个结构,形如(sum_{i=1}^{n}frac ni),我们把它记作(sum(n))。如果你做过这道题[AHOI2005]约数研究就应该不难知道这是啥。
    (sum(n)=sum_{i=1}^{n} frac ni)表示1~n中每个数的约数个数和
    所以就是对每个数求一下约数个数再取前缀和即可。
    对一个数求约数个数使用唯一分解定理,复杂度(O(nsqrt n))(预处理)
    别忘了答案式

    [ans=f(1)=sum_{d=1}^{n}mu(d)F(d)=sum_{d=1}^{n}mu(d)sum(lfloor frac {n}{d} floor)sum(lfloor frac {m}{d} floor) ]

    (O(n))处理出(mu(d))的前缀和然后直接数论分块一波
    复杂度(O(Tsqrt n))

    code

    #include<cstdio>
    #include<algorithm>
    using namespace std;
    #define ll long long
    const int N = 50000;
    int gi()
    {
    	int x=0,w=1;char ch=getchar();
    	while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
    	if (ch=='-') w=0,ch=getchar();
    	while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
    	return w?x:-x;
    }
    int mu[N+5],pri[N+5],tot,zhi[N+5];
    ll s[N+5],f[N+5];
    void Mobius()
    {
    	zhi[1]=mu[1]=1;
    	for (int i=2;i<=N;i++)
    	{
    		if (!zhi[i]) pri[++tot]=i,mu[i]=-1;
    		for (int j=1;j<=tot&&i*pri[j]<=N;j++)
    		{
    			zhi[i*pri[j]]=1;
    			if (i%pri[j]) mu[i*pri[j]]=-mu[i];
    			else {mu[i*pri[j]]=0;break;}
    		}
    	}
    	for (int i=1;i<=N;i++)
    		s[i]=s[i-1]+mu[i];
    }
    int Divide(int x)
    {
    	int p[10]={0},k[10]={0},t=0;
    	for (int i=2;i*i<=x;i++)
    		if (x%i==0)
    		{
    			p[++t]=i;
    			while (x%i==0) k[t]++,x/=i;
    		}
    	if (x>1) p[++t]=x,k[t]=1;
    	int res=1;
    	for (int i=1;i<=t;i++)
    		res*=k[i]+1;
    	return res;
    }
    int main()
    {
    	Mobius();
    	for (int i=1;i<=N;i++)
    		f[i]=f[i-1]+Divide(i);
    	int T=gi();
    	while (T--)
    	{
    		int n=gi(),m=gi();
    		if (n>m) swap(n,m);
    		int i=1;ll ans=0;
    		while (i<=n)
    		{
    			int j=min(n/(n/i),m/(m/i));
    			ans+=(s[j]-s[i-1])*f[n/i]*f[m/i];
    			i=j+1;
    		}
    		printf("%lld
    ",ans);
    	}
    	return 0;
    }
    
  • 相关阅读:
    0045算法笔记——【随机化算法】舍伍德随机化思想搜索有序表
    精进~如何成为很厉害的人
    哪些小习惯一旦养成终生受用?
    2016第24周四
    2016第24周三
    2016第24周二
    2016第24周一
    2016第23周日
    前端资源汇总
    2016第23周五
  • 原文地址:https://www.cnblogs.com/zhoushuyu/p/8251403.html
Copyright © 2011-2022 走看看