题链:
http://www.lydsy.com/JudgeOnline/problem.php?id=2154
题解:
莫比乌斯反演。
题意还是很清楚的,就不赘述了。
显然有
$ANS=sum_{i=1;j=1}^{n;m} lcm(i,j)$
化为较为熟悉的gcd形式:
$quadquad=sum_{i=1;j=1}^{n;m} frac{i imes j}{gcd(i,j)}$
令$g$为gcd的值,$F(n,m)=sum i imes j,满足1leq i leq lfloor frac{n}{g} floor,1leq j leq lfloor frac{m}{g} floor,且gcd(i,j)=1$
那么,$ANS=sum_{g=1}^{min(n,m)} g imes F(n,m)$
如果已知$F(n,m)$的值,那么ANS就可以在$Osqrt N$的复杂度内求出。
接下来看看对于一个确定的g,$F(n,m)$怎么求
即我们要求满足$1leq i leq lfloor frac{n}{g} floor,1leq j leq lfloor frac{m}{g} floor$,且$gcd(i,j)=1$的所有$i imes j$的和
还记得这个题么?●HDU 1695 GCD,让求的是满足上述条件的$(i,j)$的对数。
我们仍可以类似地去求:(令$x=lfloor frac{n}{g} floor,y=lfloor frac{m}{g} floor$)
令$t(k)为gcd(i,j)=k的i imes j的和$
$T(k)为gcd(x,y)=lambda k的i imes j的和$
显然$T(k)=sum_{k|d}{f(d)}$,即T为t的倍数关系和函数
考虑一下如何计算T(k)。
首先i的取值有:$k,2k,3k,cdots,lfloor frac{x}{k} floor k$
首先j的取值有:$k,2k,3k,cdots,lfloor frac{y}{k} floor k$
任意的两两组合都是T(k)的一部分。
令$sum(a,b)=frac{(1+a)a}{2} imesfrac{(1+b)b}{2}$
所以$T(k)=k^2 imes sum(lfloor frac{x}{k} floor,lfloor frac{y}{k} floor)$
那么由莫比乌斯反演公式的形式二(倍数关系那个式子):
$t(k)=sum_{k|d}mu(frac{d}{k})T(d)$
$quadquad=sum_{k|d}mu(frac{d}{k})d^2 imes sum(lfloor frac{x}{d} floor,lfloor frac{y}{d} floor)$
而我们要求的是t(1),所以
$t(1)=sum_{d=1}^{min(lfloor frac{n}{g} floor,lfloor frac{m}{g} floor)}mu(d)d^2 imes sum(lfloor frac{x}{d} floor,lfloor frac{y}{d} floor)$
显然,这个式子也可以在$Osqrt N$的复杂度内求出。
所以综上,时间复杂度为$O(N)$
代码:
#include<cstdio> #include<cstring> #include<iostream> #define MAXN 10000050 using namespace std; const int mod=20101009; int mu[MAXN],pmui2[MAXN]; void Sieve(int n){ static bool np[MAXN]; static int prime[MAXN],pnt; mu[1]=pmui2[1]=1; for(int i=2;i<=n;i++){ if(!np[i]) prime[++pnt]=i,mu[i]=-1; for(int j=1;j<=pnt&&i<=n/prime[j];j++){ np[i*prime[j]]=1; if(i%prime[j]) mu[i*prime[j]]=-mu[i]; else{mu[i*prime[j]]=0; break;} } pmui2[i]=(pmui2[i-1]+1ll*mu[i]*i%mod*i%mod)%mod; } } int sum(int n,int m){ return (1ll*(1+n)*n/2%mod)*(1ll*(1+m)*m/2%mod)%mod; } int F(int n,int m){ int mini=min(n,m),ret=0; for(int d=1,last;d<=mini;d=last+1){ last=min(n/(n/d),m/(m/d)); ret=(1ll*ret+1ll*(pmui2[last]-pmui2[d-1]+mod)%mod*sum(n/d,m/d)%mod)%mod; } return ret; } int main(){ int n,m,mini,ans=0; scanf("%d%d",&n,&m); mini=min(n,m); Sieve(mini); for(int g=1,last;g<=mini;g=last+1){ last=min(n/(n/g),m/(m/g)); ans=(1ll*ans+1ll*(g+last)*(last-g+1)/2%mod*F(n/g,m/g)%mod)%mod; } printf("%d ",(ans+mod)%mod); return 0; }