bzoj2154: Crash的数字表格

链接

　　http://www.lydsy.com/JudgeOnline/problem.php?id=2154

题解

　　这道题当时觉得太难了不敢做，今天在__stdcall大爷的帮助下学会了。

　　其实我就是不懂套路，莫比乌斯题的套路有两种，一种是用∑d|xμ(d)，即莫比乌斯函数的性质。还有一种就是这道题的莫比乌斯反演。

　　下面推式子。

　　ans=∑i=1n∑j=1mlcm(i,j)　　=∑d=1n1d∑d|in∑d|jm[gcd(i,j)==d]

　　注意这里千万别化成[gcd(i,j)==1]的形式，否则后面就化不动了。

　　令f(n,m,d)=∑d|in∑d|jm[gcd(i,j)==d]　　s(n,m,d)=∑d|in∑d|jmij　　=∑i=1⌊nd⌋∑j=1⌊md⌋idjd　　=d2∑i=1⌊nd⌋∑j=1⌊md⌋ij　　=d2⌊nd⌋(⌊nd⌋+1)2⌊md⌋(⌊md⌋+1)2

　　显然

　　s(n,m,d)=∑d|xnf(n,m,x)

　　莫比乌斯反演

　　f(n,m,d)=∑d|xnμ(xd)s(n,m,x)　　=∑i=1⌊nd⌋μ(i)s(n,m,id)　　=∑i=1⌊nd⌋μ(i)i2d2⌊nid⌋(⌊nid⌋+1)2⌊mid⌋(⌊mid⌋+1)2

　　那么

　　ans=∑d=1n1df(n,m,d)　　=∑d=1n1d∑i=1⌊nd⌋μ(i)i2d2⌊nid⌋(⌊nid⌋+1)2⌊mid⌋(⌊mid⌋+1)2　　=∑d=1nd∑i=1⌊nd⌋μ(i)i2⌊nid⌋(⌊nid⌋+1)2⌊mid⌋(⌊mid⌋+1)2

　　显然μ(i)i2可以前缀和预处理。最后的时间复杂度是O(n)(这个是我找will7101用微积分算出来的)

代码

//莫比乌斯反演
#include <cstdio>
#include <algorithm>
#define maxn 10000001
#define mod 20101009
#define ll long long
using namespace std;
int mu[maxn], prime[maxn/10], N, M, s[maxn];
bool mark[maxn];
ll f(ll x){return x*(x+1)/2%mod;}
void shai()
{
int i, j;
mu[1]=1;
for(i=2;i<=N;i++)
{
if(!mark[i])prime[++prime[0]]=i,mu[i]=-1;
for(j=1;j<=prime[0] and i*prime[j]<=N;j++)
{
mark[i*prime[j]]=1;
if(i%prime[j]==0){mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
for(i=1;i<=N;i++)s[i]=((ll)s[i-1]+(ll)mu[i]*i%mod*i)%mod;
}
void work()
{
ll ans=0, t;
int d, i, last, n, m, dd;
for(d=1;d<=N;d=dd+1)
{
dd=min(N/(N/d),M/(M/d));
n=N/d, m=M/d;
for(i=1,t=0;i<=n;i=last+1)
{
last=min(n/(n/i),m/(m/i));
t=(t+(s[last]-s[i-1])*f(n/i)%mod*f(m/i))%mod;
}
ans=(ans+(f(dd)-f(d-1))*t)%mod;
}
printf("%lld\n",(ans+mod)%mod);
}
int main()
{
scanf("%d%d",&N,&M);
if(N>M)swap(N,M);
shai();
work();
return 0;
}