BZOJ 2154 + 2693 莫比乌斯反演 （好题

BZOJ 2154题目链接
BZOJ 2693题目链接（权限题）

再送上一道类似推导思路的题：
NJU 1017题目链接

参考文档
（用了一整天 终于刷完文档例题了 qwq）

---

题意：
给定$n,m(<=1e7)$ 求$$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)$$
第一题是单组数据，第二题是多组数据。

---

思路：

开始写了一个$O(n\sqrt{n})$的算法 果断T。
下面试着推导一个$O(n)$甚至$O(\sqrt{n})$的计算表达式。

设答案为$Ans$，由$lcm$定义，可得：
$$Ans=\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{gcd(i,j)}$$

但这么直接计算复杂度是$O(nm)$严重爆炸
故考虑枚举$gcd(i,j)$(设为$g$)的值，想办法应用莫比乌斯反演进行化简。
$$Ans=\sum _{g=1}^{min(n,m)}(\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{g}(gcd(i,j)==g))$$

进一步化简：(转换枚举变量，使$gcd=1$)
$$Ans=\sum _{g=1}^{min(n,m)}(\sum _{i=1}^{n/g}\sum _{j=1}^{m/g}\frac{g^2\ast ij}{g}(gcd(i,j)==1))$$
得：
$$Ans=\sum _{g=1}^{min(n,m)}g(\sum _{i=1}^{n/g}\sum _{j=1}^{m/g}ij(gcd(i,j)==1))$$

令$$f(x,y)=\sum _{i=1}^x\sum _{j=1}^{y}ij(gcd(i,j)==1)$$

则：
$$Ans=\sum _{g=1}^{min(n,m)}gf(\frac{n}{g},\frac{m}{g})$$

故问题关键便在于如何来求解$f(x,y)$
此处介绍两种推导方式，一种是利用莫比乌斯函数的性质，另一种则是直接使用变换式。

---

利用${\sum }_{d|n}\mu (d)=(n==1)$

因$$f(x,y)=\sum _{i=1}^x\sum _{j=1}^{y}ij(gcd(i,j)==1)$$
将$gcd(i,j)$代入性质表达式中的$n$，得：
$$f(x,y)=\sum _{i=1}^x\sum _{j=1}^{y}ij\sum _{d|gcd(i,j)}\mu (d)$$

因$d|gcd(i,j)$，故$i,j$一定是$d$的倍数
设$i=k_{1}dj=k_{2}d$
转换枚举变量为$d$，则：
$$f(x,y)=\sum _{d=1}^{min(x,y)}\sum _{k_1=1}^{\frac{x}{d}}\sum _{k_2=1}^{\frac{y}{d}}{k}_1{k}_2{d}^2\mu (d)$$
即：
$$f(x,y)=\sum _{d=1}^{min(x,y)}{d}^2\mu (d)\sum _{k_1=1}^{\frac{x}{d}}\sum _{k_2=1}^{\frac{y}{d}}{k}_1{k}_2$$

设$$sum(x,y)=\sum _{i=1}^x\sum _{j=1}^{y}ij=\sum _{i=1}^{x}i\sum _{j=1}^{y}j=\frac{x(x+1)}{2}\frac{y(y+1)}{2}$$

故
$$f(x,y)=\sum _{d=1}^{min(x,y)}{d}^2\mu (d)sum(\frac{x}{d},\frac{y}{d})$$

故我们可以$O(n)$预处理$d^2\mu (d)$的前缀和，然后两次$O(\sqrt{n})$的分块除法，总复杂度$O(n)$下解决问题。

下面再说说如何直接利用反演表达式推出此公式。

---

利用$F(n)={\sum }_{n|d}f(d)=>f(n)={\sum }_{n|d}\mu (\frac{d}{n})F(d)$

因$$f(x,y)=\sum _{i=1}^x\sum _{j=1}^{y}ij(gcd(i,j)==1)$$
不妨设
$$g(n)=\sum _{i=1}^x\sum _{j=1}^{y}ij(gcd(i,j)==n)$$
故$g(1)=f(x,y)$

又设
$$G(n)=\sum _{i=1}^x\sum _{j=1}^{y}ij(n|gcd(i,j))$$

故可得：
$$G(n)=\sum _{n|d}g(d)=>g(n)=\sum _{n|d}\mu (\frac{d}{n})G(d)$$
所以：
$$g(1)=\sum _{d=1}^{min(x,y)}\mu (d)G(d)$$
故此时的关键在于$G$ 函数的求解

同上第一个推导，设：
$$sum(x,y)=\sum _{i=1}^x\sum _{j=1}^{y}ij=\sum _{i=1}^{x}i\sum _{j=1}^{y}j=\frac{x(x+1)}{2}\frac{y(y+1)}{2}$$

由$G$定义，易得：
$$G(d)=\sum _{i=1}^{x/d}\sum _{j=1}^{y/d}ij\ast d^2=d^{2}sum(\frac{x}{d},\frac{y}{d})$$

综上知：
$$f(x,y)=g(1)=\sum _{d=1}^{min(x,y)}{d}^2\mu (d)sum(\frac{x}{d},\frac{y}{d})$$

---

两种推导方式推出的结果已经足够解决第一题：
代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;

const ll mod = 20101009;
const int A = 1e7 + 10;
bool vis[A];
int pri[A],mu[A],tot;
int n,m;
ll sum[A];

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

ll S(ll x,ll y){
    x %= mod;y %= mod;
    ll res1 = x*(x+1)/2%mod;
    ll res2 = y*(y+1)/2%mod;
    return res1*res2%mod;
}

ll solve(ll x,ll y){
    ll last,res = 0;
    for(ll i=1 ;i<=min(x,y) ;i=last+1){
        last = min(x/(x/i),y/(y/i));
        res = (res + (sum[last] - sum[i-1])*S(x/i,y/i)%mod)%mod;
    }
    if(res<0) res+=mod;
    return res;
}

int main(){
    scanf("%d%d",&n,&m);if(n>m) swap(n,m);
    init();
    ll last,ans = 0;
    for(ll i=1 ;i<=n ;i=last+1){
        last = min(n/(n/i),m/(m/i));
        ans = (ans + (i+last)*(last-i+1)/2%mod*solve(n/i,m/i)%mod)%mod;
    }
    if(ans < 0) ans += mod;
    printf("%lld\n",ans);
}

---

但对于有多组输入的第二题，$O(n)$的复杂度还是不能在时限内得出答案。
让我们再仔细观察一下最终的结果表达式：

$$Ans=\sum _{g=1}^{min(n,m)}g\sum _{d=1}^{min(\frac{n}{g},\frac{m}{g})}{d}^2\mu (d)sum(\frac{n}{dg},\frac{m}{dg})$$

令$D=dg$我们再一次转化枚举变量。
$$Ans=\sum _{D=1}^{min(n,m)}sum(\frac{n}{D},\frac{m}{D})\sum _{d|D}{d}^2\mu (d)\frac{D}{d}$$

易发现，后面一个求和的表达式是一个可以线性筛出的积性函数。
故我们可以$O(n)$预处理，然后每组数据$O(\sqrt{n})$得出答案。
至此 该题大功告成。

同样该转化式也可以用在优化第一题：
代码：（第一题 BZOJ 2154）

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
 
const ll mod = 20101009;
const int A = 1e7 + 10;
bool vis[A];
int pri[A],tot;
ll sum[A],h[A];
int n,m;
 
void init(){
    tot = 0;h[1] = 1;
    for(ll i=2 ;i<=n ;i++){
        if(vis[i] == 0){pri[++tot] = i;h[i] =(i-1LL*i*i)%mod;}
        for(int j=1 ;j<=tot && i*pri[j] <= n ;j++){
            vis[i*pri[j]] = 1;
            if(i%pri[j] == 0){h[i*pri[j]] = 1LL*pri[j]*h[i]%mod;break;}
            h[i*pri[j]] = h[pri[j]]*h[i]%mod;
        }
    }
    for(int i=1 ;i<=n ;i++) sum[i] = ((h[i] + sum[i-1])%mod + mod)%mod;
}
 
ll S(ll n,ll m){
    n %= mod;m %= mod;
    return ((n*(n+1)/2)%mod) * ((m*(m+1)/2)%mod)%mod;
}
 
int main(){
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    init();
    ll ans = 0;
    int last;
    if(n>m) swap(n,m);
    for(int i=1 ;i<=n ;i=last+1){
        last = min(n/(n/i),m/(m/i));
        ans = (ans + S(n/i,m/i)*(sum[last] - sum[i-1])%mod)%mod;
    }
    printf("%lld\n",(ans%mod+mod)%mod);
    return 0;
}

---

代码：（第二题 BZOJ 2693）

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
 
const ll mod = 100000009;
const int A = 1e7 + 10;
bool vis[A];
int pri[A],tot;
ll sum[A],h[A];
 
void init(){
    tot = 0;h[1] = 1;
    for(ll i=2 ;i<A ;i++){
        if(vis[i] == 0){pri[++tot] = i;h[i] =(i-1LL*i*i)%mod;}
        for(int j=1 ;j<=tot && i*pri[j] < A ;j++){
            vis[i*pri[j]] = 1;
            if(i%pri[j] == 0){h[i*pri[j]] = 1LL*pri[j]*h[i]%mod;break;}
            h[i*pri[j]] = h[pri[j]]*h[i]%mod;
        }
    }
    for(int i=1 ;i<A ;i++) sum[i] = ((h[i] + sum[i-1])%mod + mod)%mod;
}
 
ll S(ll n,ll m){
    n %= mod;m %= mod;
    return ((n*(n+1)/2)%mod) * ((m*(m+1)/2)%mod)%mod;
}
 
int main(){
    init();
    int T;
    scanf("%d",&T);
    while(T--){
        int n,m;
        scanf("%d%d",&n,&m);
        ll ans = 0;
        int last;
        if(n>m) swap(n,m);
        for(int i=1 ;i<=n ;i=last+1){
            last = min(n/(n/i),m/(m/i));
            ans = (ans + S(n/i,m/i)*(sum[last] - sum[i-1])%mod)%mod;
        }
        printf("%lld\n",(ans%mod+mod)%mod);
    }
    return 0;
}