[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} floor}_{x=1}sum^{lfloor frac{m}{p} floor}_{y=1}i*j*p[gcd(x,y)==1])

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

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

(=sum_{p}p*sum_dmu(d)*d^2*sum^{lfloor frac{n}{p} floor}_{d|x}sum^{lfloor frac{m}{p} floor}_{d|y}frac{x}{d}*frac{y}{d})

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

原式(=sum_{p}p*sum_dmu(d)*d^2*sum^{lfloor frac{n}{pd} floor}_{i=1}sum^{lfloor frac{m}{pd} floor}_{j=1}i*j)

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

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

(=sum_{p}p*sum_dmu(d)*d^2*sum(lfloor frac{n}{pd} floor,lfloor frac{m}{pd} floor))

对于(sum_{p}p*sum_dmu(d)*d^2)这一部分可以预处理出前缀和，(sum(lfloor frac{n}{pd} floor,lfloor frac{m}{pd} floor))可以(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
",ans);
    }
    return 0;
}