[bzoj] 2693 jzptab || 莫比乌斯反演

原题

求\(\sum^{n}_{x=1}\sum^{m}_{y=1}lcm(x,y)\)

---

\(\sum^{n}_{x=1}\sum^{m}_{y=1}lcm(x,y)\)

\(=\sum_{p}\sum^{\lfloor \frac{n}{p} \rfloor}_{x=1}\sum^{\lfloor \frac{m}{p} \rfloor}_{y=1}i*j*p[gcd(x,y)==1]\)

\(=\sum_{p}\sum^{\lfloor \frac{n}{p} \rfloor}_{x=1}\sum^{\lfloor \frac{m}{p} \rfloor}_{y=1}\sum_{d|gcd(x,y)}\mu(d)*x*y*p\)

\(=\sum_{p}p*\sum_d\sum^{\lfloor \frac{n}{p} \rfloor}_{d|x}\sum^{\lfloor \frac{m}{p} \rfloor}_{d|y}\sum_{d|gcd(x,y)}\mu(d)*x*y\)

\(=\sum_{p}p*\sum_d\mu(d)*d^2*\sum^{\lfloor \frac{n}{p} \rfloor}_{d|x}\sum^{\lfloor \frac{m}{p} \rfloor}_{d|y}\frac{x}{d}*\frac{y}{d}\)

设i=x/d,j=y/d

原式\(=\sum_{p}p*\sum_d\mu(d)*d^2*\sum^{\lfloor \frac{n}{pd} \rfloor}_{i=1}\sum^{\lfloor \frac{m}{pd} \rfloor}_{j=1}i*j\)

因为i，j是连续的，所以为等差数列求和：

\(sum(i,j)=\sum^{\lfloor \frac{n}{pd} \rfloor}_{i=1}\sum^{\lfloor \frac{m}{pd} \rfloor}_{j=1}i*j=((\lfloor \frac{n}{pd} \rfloor)*(\lfloor \frac{n}{pd} \rfloor+1)/2)*((\lfloor \frac{m}{pd} \rfloor)*(\lfloor \frac{m}{pd} \rfloor+1)/2)\)

\(=\sum_{p}p*\sum_d\mu(d)*d^2*sum(\lfloor \frac{n}{pd} \rfloor,\lfloor \frac{m}{pd} \rfloor)\)

对于\(\sum_{p}p*\sum_d\mu(d)*d^2\)这一部分可以预处理出前缀和，\(sum(\lfloor \frac{n}{pd} \rfloor,\lfloor \frac{m}{pd} \rfloor)\)可以\(O(\sqrt(n)\)枚举。

#include<cstdio>
#include<algorithm>
#define N 10000010
#define mod 100000009ll
typedef long long ll;
using namespace std;
ll pri[N>>3],tot,g[N],sum[N],t,n,m,ans;
bool vis[N];

void init()
{
    g[1]=sum[1]=1;
    for (ll i=2;i<N;i++)
    {
	if (!vis[i])
	{
	    pri[++tot]=i;
	    g[i]=((-(ll)(i-1)*i)%mod+mod)%mod;
	}
	for (ll j=1;j<=tot && pri[j]*i<N;j++)
	{
	    vis[i*pri[j]]=1;
	    if (i%pri[j]==0)
	    {
		g[i*pri[j]]=g[i]*pri[j]%mod;
		break;
	    }
	    g[i*pri[j]]=g[i]*g[pri[j]]%mod;
	}
	sum[i]=sum[i-1]+g[i];
    }
}

ll F(ll x,ll y)
{
    return ((ll)(x+1)*x/2)%mod*(((ll)(y+1)*y/2)%mod)%mod;
}

int main()
{
    init();
    scanf("%lld",&t);
    while (t--)
    {
	scanf("%lld%lld",&n,&m);
	if (n>m) swap(n,m);
	ans=0;
	for (int i=1,last=0;i<=n;i=last+1)
	{
	    last=min(n/(n/i),m/(m/i));
	    ans+=F(n/i,m/i)*(sum[last]-sum[i-1]+mod)%mod;
	    ans%=mod;
	}
	printf("%lld\n",ans);
    }
    return 0;
}