Bzoj2154 crash的数字表格

Posted by yjjr's blog on December 15, 2017

标签:莫比乌斯反演

题目

题目传送门

Description

今天的数学课上,Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b,LCM(a, b)表示能同时被a和b整除的最小正整数。例如,LCM(6, 8) = 24。回到家后,Crash还在想着课上学的东西,为了研究最小公倍数,他画了一张NM的表格。每个格子里写了一个数字,其中第i行第j列的那个格子里写着数为LCM(i, j)。一个45的表格如下: 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

【数据规模和约定】

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

分析

因为定理:$a∗b=gcd(a,b)∗lcm(a,b)$(不会的可以参考我的博文http://blog.csdn.net/qwerty1125/article/details/78714085

所以:$ ans=\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)=\sum_{i=1}^n\sum_{j=1}^m (i*j)/gcd(i,j)$

我们可以枚举$d=gcd(i,j)$

令$F(x,y)=\sum_{1<=i<=x;1<=j<=y;d=1;} ij$ 令$Sum(x,y)=\sum_{i=1}^x\sum_{j=1}^y ij=(x(x+1)/2)(y*(y+1)/2)$

然后套用莫比乌斯反演$ans=\sum_{d=1}^{min(n,m)} d*F(\lfloor n/d \rfloor \lfloor m/d \rfloor)$

$F(x,y)=\sum_{i=1}^{min(x,y)} i^2\mu(i)Sum(\lfloor x/i \rfloor ,\lfloor y/i \rfloor)$

然后两次分块计算前缀和就好了

code

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define dep(i,a,b) for(int i=a;i>=b;i--)
#define ll long long
#define mem(x,num) memset(x,num,sizeof x)
#define mod 20101009
using namespace std;
inline ll read()
{
    ll f=1,x=0;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
const int maxn=1e7+6;
int cnt,mu[maxn],prime[maxn],s[maxn];
ll ans=0,n,m;
bool is_prime[maxn];
ll sum(ll x,ll y)
{
	return ((x*(x+1)/2)%mod)*((y*(y+1)/2)%mod)%mod;
}
void gets()
{
	mu[1]=1;
	rep(i,2,min(n,m))
	{
		if(!is_prime[i])prime[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&prime[j]*i<=min(n,m);j++)
		{
			is_prime[prime[j]*i]=1;
			mu[prime[j]*i]=-mu[i];
			if(i%prime[j]==0){mu[prime[j]*i]=0;break;}
		}
	}
	for(ll i=1;i<=min(n,m);i++)
 		s[i]=(s[i-1]+(i*i*mu[i])%mod)%mod;
}

ll F(ll x,ll y)
{
	ll ans=0,last;
	for(ll i=1;i<=min(x,y);i=last+1)
	{
		last=min(x/(x/i),y/(y/i));
		ans=(ans+(s[last]-s[i-1])*sum(x/i,y/i)%mod)%mod;
	}
	return ans;
}
int main()
{
	n=read();m=read();
	gets();
	ll last;
	for(ll d=1;d<=min(n,m);d=last+1)
	{
		last=min(n/(n/d),m/(m/d));
		ans=(ans+(d+last)*(last-d+1)/2%mod*F(n/d,m/d)%mod)%mod;
	}
	printf("%lld",(ans+mod)%mod);
	return 0;
}
本文可以转载,但必须附上原文链接,否则你会终生找不到妹子!!!欢迎关注我的CSDN: ahyjjr