bzoj 2154 莫比乌斯反演

题目链接：哆啦A梦传送门


题意：给n，m。
求：${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$。


参考博客链接：https://www.cnblogs.com/cjyyb/p/8253033.html)
设 n<m:$$\begin{array}{rl}ans& =\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)\\ & =\sum _{i=1}^n\sum _{j=1}^m(\frac{i\ast j}{gcd(i,j)})\\ & =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d](\frac{i\ast j}{d})\\ & =\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}[gcd(i,j)=1](i\ast j)\end{array}$$


$我们再设f(d)={\sum }_{i=1}^x{\sum }_{j=1}^y[gcd(i,j)=d](i\ast j)$


我们根据莫比乌斯第二描述:
$$\begin{gathered} F(n)=\sum _{n|d}f(d) \\ f(n)=\sum _{n|d}u(\frac{d}{n})F(d) \end{gathered}$$
可得：
f(d)为满足gcd(i,j)=d的i,j的乘积和。
那F(1) = f(1) + f(2) + f(3) + …

F(2) = f(2) + f(4) + f(6) +…

我们可以看出F(d)就是满足gcd(x,y)为d的倍数的x,y的乘积和

那F(d)的公式就容易求了

因为

F(1) = f(1) + f(2) + f(3) + …

所以

f(1) = μ(1)*F(1) + μ(2)*F(2) + μ(3)*F(3) + …


即：
$\begin{array}{rl}F(n)& =\sum _{n|d}f(d)\\ & =\sum _{i=1}^x\sum _{j=1}^y[n|gcd(i,j)=d](i\ast j)\\ & =n^2\sum _{i=1}^{x/n}\sum _{j=1}^{y/n}[1|gcd(i,j)=d](i\ast j)\end{array}$


我们再设:$sum={\sum }_{i=1}^a{\sum }_{j=1}^b(i\ast j)=\frac{a\ast (a+1)}{2}\ast \frac{b\ast (b+1)}{2}$
显然永远满足$[1|gcd(i,j)=d]。$


最后我们反演一下：
$$\begin{array}{rl}f(1)& =\sum _{n|d}u(d)F(d)\\ & =\sum _{d=1}^{x}u(d)F(d)\\ & =\sum _{d=1}^{x}u(d)d^2\sum _{i=1}^{x/d}\sum _{j=1}^{y/d}(i\ast j)\end{array}$$

即:
$$\begin{gathered} ans=\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}[gcd(i,j)=1](i\ast j) \\ f(1)=\sum _{d=1}^{x}u(d)d^2\sum _{i=1}^{x/d}\sum _{j=1}^{y/d}(i\ast j) \end{gathered}$$
最终得：$ans={\sum }_{d=1}^{n}d{\sum }_{i=1}^{n/d}u(i)\ast i^{2}sum(\frac{n}{d\ast i},\frac{m}{d\ast i})$
根据这两个式子，我们先预处理$u(i)\ast i^2$，再分块下就解决了。

#include<bits/stdc++.h>
using namespace std;

typedef long long LL;
const int N = 1e7 + 5;
int mu[N], vis[N], prime[N];
int tot;
int n,m;
const LL mod=20101009;
LL sum[N];
void init(){
    mu[1] = 1;
    for(int i = 2; i<=n; i ++){
        if(!vis[i]){
            prime[tot ++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot && i * prime[j] <=n; j ++){
            vis[i * prime[j]] = 1;
            if(i % prime[j]) mu[i * prime[j]] = -mu[i];
            else{
                mu[i * prime[j]] = 0;
                break;
            }
        }
    }
    
    ///预处理u(i)*i*i的前缀和
    sum[0]=0LL;
    for(int i=1;i<=n;i++)
        sum[i]=(sum[i-1]+(1LL*mu[i]*i%mod*i%mod))%mod;
}


LL ans(LL x,LL y)
{
    x=(x*(x+1)/2)%mod;
    y=(y*(y+1)/2)%mod;
    return x*y%mod;
}

///整数分块求函数f(1)
LL f(int x,int y)
{
    int r;
    LL re=0;
    for(int l=1;l<=x;l=r+1){
        r=min(x/(x/l),y/(y/l));
        re=(re+(sum[r]-sum[l-1]+mod)%mod*ans((LL)x/l,(LL)y/l)%mod)%mod;
    }
    return re;
}
int main(){

    scanf("%d%d",&n,&m);

    if(n>m) swap(n,m);
     init();

     LL re=0;
     int r;
     ///整数分块求sum
     for(int l=1;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        re=(re+1LL*(l+r)*(r-l+1)/2%mod*f(n/l,m/l)%mod)%mod;
     }
     printf("%lld\n",re);
     return 0;
}

我的标签：做个有情怀的程序员。