BZOJ 2693: jzptab [莫比乌斯反演 线性筛]

题意：求\(\sum_{i=1}^n \sum_{j=1}^m lcm(i,j)\)

---

就是$$\sum_{i=1}^n \sum_{j=1}^m \frac{i j}{gcd(i,j)}$$

套路推♂倒

\[\begin{align*} \sum_{i=1}^n \sum_{j=1}^m \frac{i j}{gcd(i,j)} &=\sum_{d=1}^n \sum_{i=1}^{\frac{n}{d}} \sum_{j=1}^{\frac{m}{d}}ijd[gcd(i,j)=1]\\ &=\sum_{d=1}^n \sum_{e=1}^n \mu(e) \sum_{i=1}^{\frac{n}{de}} \sum_{j=1}^{\frac{m}{de}} ijde^2 \\ &= \sum_{D=1}^n D\sum_{d|D} d\ \mu(d) \ f(\frac{n}{D}, \frac{m}{D}) \\ \end{align*} \]

\(f(n,m) = \frac{n(n+1)}{2} \frac{m(m+1)}{2}\) 就是等比数列求和


剩下就是要处理出\(D\sum_{d|D} d\ \mu(d)\)的前缀和了，就是\(D \cdot (id \cdot \mu) * 1\)，是积性函数

令\(g(i) = i \cdot ((id \cdot \mu) * 1)(i)\)
\(g(1) = 1,\ g(p) = p(1-p)\)
考虑 \(p \nmid i\) , 因为带着\(\mu\)，所以因数\(d\)必须是sf才有贡献，多一个\(p\)之后肯定是0啊...所以只有前面的系数\(i\)改变了，\(g(ip) = g(i)\cdot p\)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=1e7+5, P=1e8+9;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}
    return x*f;
}
 
int n, m, k;
int notp[N], p[N]; ll mu[N], g[N];
void sieve(int n) {
    mu[1] = 1; g[1] = 1;
    for(int i=2; i<=n; i++) {
        if(!notp[i]) p[++p[0]] = i, mu[i] = -1, g[i] = (ll)i*(1-i) %P;
        for(int j=1; j<=p[0] && i*p[j]<=n; j++) {
            notp[i*p[j]] = 1;
            if(i%p[j] == 0) {
                mu[i*p[j]] = 0;
                g[i*p[j]] = g[i]*p[j] %P;
                break;
            }
            g[i*p[j]] = g[i] * g[p[j]] %P;
        }
    }
    for(int i=1; i<=n; i++) g[i] = (g[i] + g[i-1]) %P;
}
ll f(ll a, ll b) {return (a*(a+1)/2 %P) * (b*(b+1)/2 %P) %P;}
ll cal(int n, int m) {
    ll ans=0; int r;
    for(int i=1; i<=n; i=r+1) {
        r = min(n/(n/i), m/(m/i));
        ans += (g[r] - g[i-1]) * f(n/i, m/i) %P;
    }
    return (ans + P) %P;
}
int main() {
    //freopen("in","r",stdin);
    sieve(N-1);
    int T=read();
    while(T--) {
        n=read(); m=read();
        if(n>m) swap(n, m);
        printf("%lld\n", ( cal(n, m) + P) %P);
    }
}