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

题目大意：

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

n,m≤10000000n,m≤10000000

多组数据 T<=T<=

题解：

与这道题↓一样，但因为是多组数据，所以需要优化。

【BZOJ2154】Crash的数字表格

题解前面请看上面的博客，这里只讲优化。

因为

ans=∑d=1min(n,m)d×f(⌊nd⌋,⌊md⌋,1)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)2f(x,y,1)=∑k=1min(x,y)μ(k)k2⌊xk⌋(⌊xk⌋+1)2⌊yk⌋(⌊yk⌋+1)2

则

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

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

把kk和dd合成为D=kdD=kd

ans=∑D=1min(n,m)⌊nD⌋(⌊nD⌋+1)2⌊mD⌋(⌊mD⌋+1)2∑k|DDkμ(k)k2ans=∑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)ans=∑D=1min(n,m)⌊nD⌋(⌊nD⌋+1)2⌊mD⌋(⌊mD⌋+1)2∑k|DDk×μ(k)

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

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

设

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

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

代码：

#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;
}