求(sum_{i=1}^Nsum_{j=1}^Mlcm(i,j),n,mleq 10^7)
解
设(N<M),显然有
[sum_{i=1}^Nsum_{j=1}^Mfrac{ij}{gcd(i,j)}=sum_{d=1}^Nfrac{1}{d}sum_{i=1}^Nsum_{j=1}^Mij(gcd(i,j)==d)
]
设
[f(k)=sum_{i=1}^Nsum_{j=1}^Mij(gcd(i,j)==k)
]
[F(k)=sum_{i=1}^{N}sum_{j=1}^Mij(k|gcd(i,j))=k^2sum_{i=1}^{[N/k]}sum_{j=1}^{[M/k]}ij
]
设(dc[k]=sum_{i=1}^ki=frac{(1+k) imes k}{2},F(k)=k^2dc(N/k)dc(M/k))
由Mobius反演定理我们有
[f(k)=sum_{k|d}F(d)mu(d/k)
]
代入有
[ans=sum_{d=1}^Nfrac{1}{d}sum_{d|x}^{}x^2dc(N/x)dc(M/x)mu(x/d)=
]
[sum_{d=1}^Ndsum_{x=1}^{[N/d]}dc(N/xd)dc(M/xd)x^2mu(x)
]
维护出后式(x^2mu(x)),两次整除分块即可,不难得知时间复杂度(O(n))。
顺便提一下，如果(N,M)很小，我们可以变成一下形式，变为(O(nlogn+Tsqrt{n}))(T为询问组数)。
[sum_{x=1}^Ndc(N/x)dc(M/x)sum_{d|x}mu(x/d)frac{x^2}{d}
]
参考代码：
#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define yyb 20101009
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[10000001];
int prime[750000],pt,mb[10000001];
void prepare(int);
il int min(int,int),dx(int,int);
int main(){
int n,m,nd,md,ndx,mdx,i,ij,j,
jj,ans1(0),ans2;
scanf("%d%d",&n,&m);
if(n>m)swap(n,m);prepare(m);
for(i=1;i<=n;i=ij+1){
ij=min(n/(n/i),m/(m/i));
ans2&=0,nd=n/i,md=m/i;
for(j=1;j<=nd;j=jj+1)
jj=min(nd/(nd/j),md/(md/j)),
(ans2+=(ll)(mb[jj]-mb[j-1])*dx(1,nd/j)%yyb*dx(1,md/j)%yyb)%=yyb;
(ans1+=(ll)ans2*dx(i,ij)%yyb)%=yyb;
}printf("%d",(ans1+yyb)%yyb);
return 0;
}
il int dx(int a,int b){
return (ll)(a+b)*(b-a+1)/2%yyb;
}
void prepare(int n){
int i,j;check[1]|=mb[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])prime[++pt]=i,mb[i]=-1;
for(j=1;j<=pt&&prime[j]<=n/i;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j]))break;
mb[i*prime[j]]=-mb[i];
}
}for(i=1;i<=n;++i)mb[i]=((ll)mb[i]*i%yyb*i+mb[i-1])%yyb;
}
il int min(int a,int b){
return a<b?a:b;
}