zoukankan      html  css  js  c++  java
  • 【刷题】BZOJ 2154 Crash的数字表格

    Description

    今天的数学课上,Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b,LCM(a, b)表示能同时被a和b整除的最小正整数。例如,LCM(6, 8) = 24。回到家后,Crash还在想着课上学的东西,为了研究最小公倍数,他画了一张N*M的表格。每个格子里写了一个数字,其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下: 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格,Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题:这个表格中所有数的和是多少。当N和M很大时,Crash就束手无策了,因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大,Crash只想知道表格里所有数的和mod 20101009的值。

    Input

    输入的第一行包含两个正整数,分别表示N和M。

    Output

    输出一个正整数,表示表格中所有数的和mod 20101009的值。

    Sample Input

    4 5

    Sample Output

    122

    HINT

    100%的数据满足N, M ≤ 10^7。

    Solution

    莫比乌斯反演,膜拜PoPoQQQ

    这题题意还是很清楚的,话不多说,推式子吧

    哦,我之前(too young too simple)了,分块的时候判上界每次都会取min,结果慢了很多,这次保证了N比M小,所以在推式子的时候不要疑惑

    [ans=sum_{i=1}^Nsum_{j=1}^Mlcm(i,j) ]

    [=sum_{i=1}^Nsum_{j=1}^Mfrac{ij}{gcd(i,j)} ]

    [=sum_{d=1}^Ndsum_{i=1}^{lfloor frac{N}{d} floor}sum_{j=1}^{lfloor frac{M}{d} floor}[gcd(i,j)=1]ij ]

    这里解释一下,为什么枚举公约数,后面还要是一个([gcd(i,j)=1]):因为当(i)(j)除了(d)以外还有公因子(k)的话,那么(gcd(i,j)=dk),而不是(d),这个式子就错了

    然后我们设(g(x,y)=sum_{i=1}^{x}sum_{j=1}^{y}[gcd(i,j)=1]ij)(s(x,y)=sum_{i=1}^{x}sum_{j=1}^{y}ij=frac{(x+1)x}{2}frac{(y+1)y}{2})

    那么

    [ans=sum_{d=1}^Ndcdot g(lfloor frac{N}{d} floor,lfloor frac{M}{d} floor) ]

    假设我们知道(g(x,y))的值,那么算这个式子,整除分块就可以了吧

    那么我们接下来的问题就是算(g(x,y))

    对于(g(x,y)),这里的求法完全是单独的,用莫比乌斯反演。

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

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

    然后我们来推(F(n))

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

    ( =sum_{ni=1}^{x}sum_{nj=1}^{y}[n|gcd(ni,nj)]nicdot nj)

    ( =sum_{i=1}^{lfloor frac{x}{n} floor}sum_{j=1}^{lfloor frac{y}{n} floor}ni cdot nj)

    ( =n^2sum_{i=1}^{lfloor frac{x}{n} floor}sum_{j=1}^{lfloor frac{y}{n} floor}ij)

    ( =n^2cdot s(lfloor frac{x}{n} floor,lfloor frac{y}{n} floor))

    反演得(f(n))

    (f(n)=sum_{n|d}mu(frac{d}{n})F(d))

    ( =sum_{d=1}^{N}mu(d)F(nd))

    ( =sum_{d=1}^{N}mu(d)cdot n^2d^2cdot s(lfloor frac{x}{nd} floor,lfloor frac{y}{nd} floor))

    那么(g(x,y)=f(1)=sum_{d=1}^{N}mu(d)cdot d^2cdot s(lfloor frac{x}{d} floor,lfloor frac{y}{d} floor))

    于是,又可以整除分块了。。。

    预处理(mu(d)cdot d^2)前缀和,分块求(g(x,y))

    所以这个题目的算法就是一个大整除分块加一个小整除分块,最后复杂度两个根号乘起来正好一个(O(N))

    #include<bits/stdc++.h>
    #define ll long long
    const int MAXN=10000000+10,Mod=20101009;
    int N,M,prime[MAXN],cnt,mu[MAXN],vis[MAXN];
    ll res,two,mus[MAXN];
    template<typename T> inline void read(T &x)
    {
    	T data=0,w=1;
    	char ch=0;
    	while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
    	if(ch=='-')w=-1,ch=getchar();
    	while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
    	x=data*w;
    }
    template<typename T> inline void write(T x,char c='')
    {
    	if(x<0)putchar('-'),x=-x;
    	if(x>9)write(x/10);
    	putchar(x%10+'0');
    	if(c!='')putchar(c);
    }
    template<typename T> inline void chkmin(T &x,T y){x=(y<x?y:x);}
    template<typename T> inline void chkmax(T &x,T y){x=(y>x?y:x);}
    template<typename T> inline T min(T x,T y){return x<y?x:y;}
    template<typename T> inline T max(T x,T y){return x>y?x:y;}
    inline ll qexp(ll a,ll b)
    {
    	ll res=1;
    	while(b)
    	{
    		if(b&1)res=res*a%Mod;
    		a=a*a%Mod;
    		b>>=1;
    	}
    	return res;
    }
    inline void init()
    {
    	two=qexp(2,Mod-2);
    	memset(vis,1,sizeof(vis));
    	vis[0]=vis[1]=0;
    	mu[1]=1;
    	for(register int i=2;i<MAXN;++i)
    	{
    		if(vis[i])
    		{
    			prime[++cnt]=i;
    			mu[i]=-1;
    		}
    		for(register int j=1;j<=cnt&&i*prime[j]<MAXN;++j)
    		{
    			vis[i*prime[j]]=0;
    			if(i%prime[j])mu[i*prime[j]]=-mu[i];
    			else break;
    		}
    	}
    	for(register int i=1;i<MAXN;++i)mus[i]=(mus[i-1]+(ll)mu[i]*(ll)i*(ll)i)%Mod;
    }
    inline ll s(int x,int y)
    {
    	return (ll)(x+1)*(ll)x%Mod*two%Mod*(ll)(y+1)%Mod*(ll)y%Mod*two%Mod;
    }
    inline ll g(int x,int y)
    {
    	ll ans=0;
    	if(x>y)std::swap(x,y);
    	for(register int i=1;;)
    	{
    		if(i>x)break;
    		int j=min(x/(x/i),y/(y/i));
    		(ans+=(mus[j]-mus[i-1]+Mod)%Mod*s(x/i,y/i)%Mod)%=Mod;
    		i=j+1;
    	}
    	return ans;
    }
    int main()
    {
    	read(N);read(M);
    	init();
    	if(N>M)std::swap(N,M);
    	for(register int i=1;;)
    	{
    		if(i>N)break;
    		int j=min(N/(N/i),M/(M/i));
    		(res+=(ll)(i+j)*(ll)(j-i+1)%Mod*two%Mod*g(N/i,M/i)%Mod)%=Mod;
    		i=j+1;
    	}
    	write(res,'
    ');
    	return 0;
    }
    
    
    
  • 相关阅读:
    新东西
    Xcode6新特性
    下载模拟器
    iOS定位和地图
    iOS,作死集锦
    ThreadLocal源码解析
    JSON Web令牌(JWT)介绍与使用
    docker已运行容器里的时区修改
    Docker图形界面管理
    ZooKeeper开机启动的俩种方式
  • 原文地址:https://www.cnblogs.com/hongyj/p/8552047.html
Copyright © 2011-2022 走看看