zoukankan      html  css  js  c++  java
  • 【BZOJ2154】Crash的数字表格(莫比乌斯反演)

    【BZOJ2154】Crash的数字表格(莫比乌斯反演)

    题面

    BZOJ
    简化题意:
    给定(n,m)
    求$$sum_{i=1}^nsum_{j=1}^mlcm(i,j)$$

    题解

    以下的一切都默认(n<m)
    我们都知道(lcm(i,j)=frac{ij}{gcd(i,j)})
    所以所求化简

    [sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)} ]

    看到(gcd(i,j))很不爽,于是就再提出来

    [sum_{d=1}^{n}sum_{i=1}^nsum_{j=1}^m[gcd(i,j)==d]frac{ij}{d} ]

    也就是

    [sum_{d=1}^{n}sum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ijd} ]

    (d)提出来

    [ans=sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ij} ]

    前面这一堆看起来管不了了
    看后面的一段

    [sum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ij} ]

    看到(n/d)这种东西很不爽呀
    就写成这样吧。。

    [sum_{i=1}^{x}sum_{j=1}^{y}[gcd(i,j)==1]{ij} ]

    这种东西怎么求?

    令$$f(d)=sum_{i=1}^{x}sum_{j=1}^{y}[gcd(i,j)==d]{ij}$$
    根据莫比乌斯反演的常见套路

    [G(d)=sum_{i=1}^{x}sum_{j=1}^{y}[d|gcd(i,j)]{ij} ]

    直接把(d)提出来

    [G(d)=d^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}[1|gcd(i,j)]{ij} ]

    (1|gcd(i,j))是显然成立的
    所以$$G(d)=d^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}{ij}$$
    这玩意明显可以(O(1))求(相当于两个等差数列相乘)

    所以,要求的东西就是$$f(1)=sum_{i=1}^xmu(i)G(i)$$

    这道题就解决了一大半了
    现在我们的复杂度是(O(nsqrt n))(O(n^2))之间
    需要继续优化

    很显然的

    [ans=sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]{ij} ]

    这个式子可以数论分块一波,复杂度少了(O(sqrt n))

    还不够

    继续看,

    [f(1)=sum_{i=1}^xmu(i)G(i) ]

    这个式子把(G(x))展开

    [f(1)=sum_{i=1}^xmu(i)i^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}{ij} ]

    还是可以数论分块
    但是要预处理(mu(i)*i^2)的前缀和

    然后复杂度就成了(O(n))
    注释掉的是没用数论分块的式子

    #include<iostream>
    #include<cstdio>
    #include<cstdlib>
    #include<cstring>
    #include<cmath>
    #include<algorithm>
    #include<set>
    #include<map>
    #include<vector>
    #include<queue>
    using namespace std;
    #define MOD 20101009
    #define MAX 12000000
    #define ll long long
    inline int read()
    {
        int x=0,t=1;char ch=getchar();
        while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
        if(ch=='-')t=-1,ch=getchar();
        while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
        return x*t;
    }
    int mu[MAX],pri[MAX],tot;
    bool zs[MAX];
    int n,m;
    int G[MAX],ans;
    int smu[MAX],sqr[MAX];
    void Getmu()
    {
    	zs[1]=true;mu[1]=1;
    	for(int i=2;i<=n;++i)
    	{
    		if(!zs[i]){pri[++tot]=i;mu[i]=-1;}
    		for(int j=1;j<=tot&&i*pri[j]<=n;++j)
    		{
    			zs[i*pri[j]]=true;
    			if(i%pri[j])mu[i*pri[j]]=-mu[i];
    			else{mu[i*pri[j]]=0;break;}
    		}
    	}
    	for(int i=1;i<=n;++i)smu[i]=(smu[i-1]+mu[i]+MOD)%MOD;
    }
    int Solve(int xx,int yy)
    {
    	long long ans=0;
    	//for(int i=1;i<=xx;++i)G[i]=1ll*i*i%MOD*1ll*(1ll*(1+xx/i)*(xx/i)/2%MOD)*(1ll*(1+yy/i)*(yy/i)/2%MOD)%MOD;
    	//for(int i=1;i<=xx;++i)ans=(ans+1ll*mu[i]*G[i]%MOD+MOD)%MOD;
    	int i=1,j;
    	while(i<=xx)
    	{
    	 	j=min(xx/(xx/i),yy/(yy/i));
    		int GG=1ll*(1ll*(1+xx/i)*(xx/i)/2%MOD)*(1ll*(1+yy/i)*(yy/i)/2%MOD)%MOD;
    		ans+=1ll*(sqr[j]-sqr[i-1])%MOD*GG%MOD;
    	 	ans%=MOD;
    	 	i=j+1;
    	}
    	return (ans+MOD)%MOD;
    }
    int main()
    {
    	n=read();m=read();
    	if(n>m)swap(n,m);
    	Getmu();
    	for(int i=1;i<=n;++i)sqr[i]=(sqr[i-1]+1ll*i*i%MOD*mu[i]%MOD+MOD)%MOD;
    	//for(int i=1;i<=n;++i)ans=1ll*((ans+1ll*i*Solve(n/i,m/i))%MOD)%MOD;
    	int i=1,j;
    	while(i<=n)
    	{
    		j=min(n/(n/i),m/(m/i));
    		int t=1ll*(i+j)*(j-i+1)/2%MOD;
    		ans=(ans+1ll*Solve(n/i,m/i)*t%MOD)%MOD;
    		i=j+1;
    	}
    	printf("%d
    ",ans);
    	return 0;
    }
    
    
  • 相关阅读:
    使用脚本改变树控件的行为 (转)点文本 收..
    (面包屑)SiteMapPath控件简化Web网站导航 (转)
    Web.config详解(转)
    SiteMap(站点地图)示例(转)
    url传递中文的解决方案总结(转)
    INI文件编程,WINAPI函数WritePrivateProfileString,GetPrivateProfileString(转帖)
    hdu 2708
    hdu 1709
    hdu 1045
    hdu 2714
  • 原文地址:https://www.cnblogs.com/cjyyb/p/8253033.html
Copyright © 2011-2022 走看看