【BZOJ2154】Crash的数字表格（莫比乌斯反演）（数学）

题目大意：

求

∑i=1n∑j=1mlcm(i,j)

题解：

ans=∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mijgcd(i,j)

设

f(x,y,d)=∑1≤i≤n1≤j≤mgcd(i,j)=dij

F(x,y,d)=∑d|kf(x,y,k)=∑1≤i≤n1≤j≤md|gcd(i,j)ij=d2∑1≤i≤⌊nd⌋1≤j≤⌊md⌋ij=d2⌊nd⌋(⌊nd⌋+1)2⌊md⌋(⌊md⌋+1)2

∴

ans=∑d=1min(n,m)d2×f(⌊nd⌋,⌊md⌋,1)d=∑d=1min(n,m)d×f(⌊nd⌋,⌊md⌋,1)

f(x,y,1)=∑d=1min(x,y)μ(k)F(x,y,k)=∑k=1min(x,y)μ(k)k2⌊xk⌋(⌊xk⌋+1)2⌊yk⌋(⌊yk⌋+1)2

可以做了

代码：

#include<cstdio>#include<algorithm>using namespace std;const long long MOD=20101009LL,MAXN=10000005LL;long long sum_mu[MAXN],prime[MAXN/3LL],pcnt;bool noprime[MAXN];void init_mu(int);long long solve(long long,long long);long long f(long long,long long);long long sum(long long,long long);int main(){    long long n,m;    scanf("%lld%lld",&n,&m);    init_mu(min(n,m));    printf("%lld\n",solve(n,m));    return 0;}void init_mu(int N){    noprime[1]=1;    sum_mu[1]=1LL;    for(long long i=2LL;i<=N;i++)    {        if(!noprime[i])        {            prime[++pcnt]=i;            sum_mu[i]=-1LL;        }        for(long long j=1;j<=pcnt&&i*prime[j]<=N;j++)        {            noprime[i*prime[j]]=1;            if(i%prime[j]==0)            {                sum_mu[i*prime[j]]=0LL;                break;            }            sum_mu[i*prime[j]]=-sum_mu[i];        }    }    for(long long i=1LL;i<=N;i++)        sum_mu[i]=(sum_mu[i-1]+((((i*i)%MOD)*sum_mu[i])%MOD))%MOD;}long long solve(long long n,long long m){    long long ret=0LL,last;    if(n>m)        swap(n,m);    for(long long d=1LL;d<=n;d=last+1)    {        last=min(n/(n/d),m/(m/d));        ret=(ret+((sum(d,last)*f(n/d,m/d))%MOD))%MOD;    }    return ret;}long long f(long long x,long long y){    if(x>y)        swap(x,y);    long long ret=0LL,last,tmp;    for(long long k=1LL;k<=x;k=last+1LL)    {        last=min(x/(x/k),y/(y/k));        tmp=(((((x/k)*(x/k+1LL))>>1LL)%MOD)        *((((y/k)*(y/k+1LL))>>1LL)%MOD))%MOD;        ret=(ret+(((sum_mu[last]-sum_mu[k-1]+MOD)%MOD)*tmp)%MOD)%MOD;    }    return ret;}long long sum(long long a,long long b){    long long ret;    ret=(((a+b)*(b-a+1LL))>>1LL)%MOD;    return ret;}