bzoj2154: Crash的数字表格

传送门

首先lcm(i,j)=i∗jgcd(i,j)

所以答案=∑ni=1∑mj=1i∗jgcd(i,j)

设d=gcd(i,j)

则答案=∑d=1∑i=1nd∑j=1mdi∗j(gcd(i,j)==1)

=∑d=1∑i=1ndi∑j=1mdj(gcd(i,j)==1)

显然(gcd(i,j)==1)是可以用莫比乌斯反演在sqrt(i,j)的时间内求出来的。

再对d进行分块就行了。

#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
const ll mo=20101009;
ll n,m,l,r,ans,mou[10000005];
ll sum(ll x,ll y){return (x*(x+1)/2%mo)*(y*(y+1)/2%mo)%mo;}
ll F(ll n,ll m){
ll ans=0,l,r;
for (l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(mou[r]-mou[l-1])%mo*sum(n/l,m/l)%mo;
ans%=mo;
}
return ans;
}
int main(){
scanf("%lld%lld",&n,&m);
if (n>m) swap(n,m);
for (ll i=2;i<=m;i++)
if (!mou[i]){
mou[i]=i;
for (ll j=i*i;j<=m;j+=i) mou[j]=i;
}
mou[1]=1;
for (ll i=2;i<=m;i++)
if (i/mou[i]%mou[i]==0) mou[i]=0;
else mou[i]=-mou[i/mou[i]];
for (ll i=1;i<=m;i++)
mou[i]=(mou[i-1]+mou[i]*i%mo*i%mo)%mo;
for (l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(l+r)*(r-l+1)/2%mo*F(n/l,m/l)%mo;
ans%=mo;
}
if (ans<0) ans+=mo;
printf("%lld",ans);
}