BZOJ2154: Crash的数字表格
题目描述
题目分析
这题就是要求:
[Ans=sum_{i=1}^nsum_{j=1}^m{lcm(i,j)}
]
首先我们知道(lcm(a,b)=frac{ab}{gcd(a,b)}),然后就可把它代进去。
[Ans=sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)}
]
然后根据套路，我们枚举(gcd)
[Ans=sum_{d=1}^{min(n,m)}sum_{i=1}^nsum_{j=1}^mfrac{ij}{d}{[gcd(i,j)=d]}
]
把(d)提出去，可以得到：
[Ans=sum_{d=1}^{min(n,m)}dsum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}ij[gcd(i,j)=1]
]
根据莫比乌斯函数有趣的性质，可以化成：
[Ans=sum_{d=1}^{min(n,m)}dsum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{i=1}^{lfloorfrac{m}{d}
floor}ijsum_{xmid gcd(i,j)}mu(x)
]
再把(x)放到外层枚举
[Ans=sum_{d=1}^{min(n,m)}dsum_{x=1}^{min(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)}mu(x)sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}ij[xmid gcd(i,j)]
]
然后再把原式由枚举(i,j)变成枚举其(x)的倍数
[Ans=sum_{d=1}^{min(n,m)}dsum_{x=1}^{min(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)}mu(x)x^2sum_{i=1}^{lfloorfrac{n}{dx}
floor}sum_{j=1}^{lfloorfrac{m}{dx}
floor}ij
]
然后发现后面的(i,j)此时可以分开计算了。
[Ans=sum_{d=1}^{min(n,m)}dsum_{x=1}^{min(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)}mu(x)x^2(sum_{i=1}^{lfloorfrac{n}{dx}
floor}i)(sum_{j=1}^{lfloorfrac{m}{dx}
floor}j)
]
这个式子看上去是能过的，实际上某谷也是可以过的，但是BZOJ会T。
这不禁让我们思考是否更加优秀的做法。
可以发现其实刚才分开的部分是一个范围内的数对之和，我们设
[f(n,m)=sum_{i=1}^{n}sum_{j=1}^{m}ij=frac{n imes (n+1)}{2}frac{m imes (m+1)}{2}
]
这样后面的一坨就变成了
[sum_{x=1}^{min(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)}mu(x)x^2f(lfloorfrac{n}{x}
floor,lfloorfrac{m}{x}
floor)
]
明显这个式子是可以整数分块解的。
我们设这类式子为(sum(n,m))
则我们的原式就变成了
[Ans=sum_{d=1}^{min(n,m)}d·sum(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)
]
这个整数分块解的。
那么使用整数分块套整数分块即可以在优秀的时间复杂度内完成本题。
是代码呢
#include <bits/stdc++.h>
using namespace std;
const int MAXN=1e7+7;
#define ll long long
const int mo=20101009;
const ll inv=10050505;
int n,m,prime[MAXN],mu[MAXN];
ll sum[MAXN],ans;
bool vis[MAXN];
inline void get_mu(int N)
{
mu[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0];j++){
if(1ll*i*prime[j]>N) break;
vis[i*prime[j]]=1;
if(i%prime[j]==0) break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=N;i++) sum[i]=(sum[i-1]+1ll*i*i%mo*mu[i])%mo;
}
inline ll G(int x,int y)
{
return (1ll*x*(x+1)/2)%mo*(1ll*y*(y+1)/2%mo)%mo;
}
inline ll Sum(int x,int y)
{
ll summ=0;
for(int l=1,r;l<=min(x,y);l=r+1){
r=min(x/(x/l),y/(y/l));
summ=(summ+1ll*(sum[r]-sum[l-1]+mo)*G(x/l,y/l)%mo)%mo;
}
return summ;
}
inline int read()
{
int x=0,c=1;
char ch=' ';
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
while(ch=='-')c*=-1,ch=getchar();
while(ch>='0'&&ch<='9')x=x*10+ch-'0',ch=getchar();
return x*c;
}
int main()
{
n=read();m=read();
get_mu(min(n,m)+1);
int mx=min(n,m);
for(int l=1,r;l<=mx;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(r-l+1)*(l+r)/2%mo*Sum(n/l,m/l)%mo)%mo;
}
printf("%lld
", (ans%mo+mo)%mo);
}