【BZOJ2693】jzptab（莫比乌斯反演）（数学）

题目大意：

求

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

n,m≤10000000
多组数据 T<=

题解：

与这道题↓一样，但因为是多组数据，所以需要优化。
【BZOJ2154】Crash的数字表格
题解前面请看上面的博客，这里只讲优化。
因为

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

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

则
∵当kd>n时，⌊nkd⌋=0，所以下式∑min(n,m)k=1不需要写为∑min(⌊nd⌋,⌊md⌋)k=1

ans=∑d=1min(n,m)d×∑k=1min(n,m)μ(k)k2⌊nkd⌋(⌊nkd⌋+1)2⌊mkd⌋(⌊mkd⌋+1)2

把k和d合成为D=kd

ans=∑D=1min(n,m)⌊nD⌋(⌊nD⌋+1)2⌊mD⌋(⌊mD⌋+1)2∑k|DDkμ(k)k2

ans=∑D=1min(n,m)⌊nD⌋(⌊nD⌋+1)2⌊mD⌋(⌊mD⌋+1)2∑k|DDk×μ(k)

∴我们可以求出∑k|DDk×μ(k)的前缀和

∑k|DDk×μ(k)=D∑k|Dk×μ(k)

设

g(D)=∑k|Dk×μ(k)

可知g一定是积性函数，可用线性筛来求。

代码：

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