【luogu P1829】Crash的数字表格 (数论分块，莫比乌斯反演）

· · 题解

题意

求解式子

ANS(n,m) = \sum_{i=1}^n\sum_{j=1}^mlcm(i,j)

其中1 \leq n,m \leq 1e7

分析

因为是单次询问，所需复杂度不高于O(n)

推导的思路：

将lcm转化成gcd

ANS(n,m) = \sum_{i=1}^n\sum_{j=1}^m \frac{i\cdot j}{gcd(i,j)}

求和中有gcd, 转变求和顺序，先以相同的gcd分组，枚举gcd。

由gcd的性质可以知道gcd(i,j)| i 并且gcd(i,j)|j

ANS(n,m) = \sum_{d=1}^n (\frac 1 d\sum_{d|i} \sum_{d|j} ([gcd(i,j)==d]\cdot i \cdot j))

观察发现i的取值是d,2d,...,\lfloor\frac n d \rfloor d，

j$的取值是$d,2d,...,\lfloor\frac m d \rfloor d

换元 i = \frac i d,j = \frac j d

ANS(n,m) = \sum_{d=1}^n d\sum_{i=1}^{\lfloor \frac n d \rfloor} \sum_{j=1}^{\lfloor \frac m d \rfloor}([gcd(i,j)==1]\cdot i \cdot j )

令

S(n,m) = \sum_{i=1}^{n} \sum_{j=1}^{m}([gcd(i,j)==1]\cdot i \cdot j )

这里是在求解互质的数对的乘积和

则

ANS(n,m) = \sum_{d=1}^n( d \cdot S(\lfloor \frac n d\rfloor,\lfloor \frac m d\rfloor))

考虑如何求解S(n,m)

因为\epsilon(n)在n为1的时候为1，其他时候为0

S(n,m) = \sum_{i=1}^{n} \sum_{j=1}^{m}(\epsilon(gcd(i,j))\cdot i \cdot j )

因为1 * \mu = \epsilon

所以用\mu 替代\epsilon

S(n,m) = \sum_{i=1}^{n} \sum_{j=1}^{m}(\sum _{d|gcd(i,j)}\mu(d))\cdot i \cdot j )

于是这里满足的关系是d | gcd(i,j) 且gcd(i,j)|i,gcd(i,j)|j，即d|i,d|j

和刚刚一样的思路，转化为先枚举因子

S(n,m) = \sum_{d=1}^{n} \mu (d) \sum_{d|i} \sum _{d|j} i \cdot j

换元 i = \frac i d,j = \frac j d

S(n,m) = \sum_{d=1}^n \mu(d) d^2\sum_{i=1}^{\lfloor \frac n d \rfloor} \sum_{j=1}^{\lfloor \frac m d \rfloor}( i \cdot j )

令

g(n,m) = \sum_{i=1}^{n} \sum_{j=1}^{m}( i \cdot j )

这里可以用公式计算出来

g(n,m) = \frac {n(n+1)} 2 \cdot \frac {m(m+1)} 2

而前半部分的\mu(d) d^2可以线性筛筛出莫比乌斯函数后，预处理前缀和求解

则

S(n,m) = \sum_{d=1}^n \mu(d)\cdot d^2\cdot g(\lfloor \frac n d \rfloor,\lfloor \frac m d \rfloor)

可以分块求解

ANS(n,m) = \sum_{d=1}^n( d \cdot S(\lfloor \frac n d\rfloor,\lfloor \frac m d\rfloor))

再次分块求解

求得最后答案

代码

#include<bits/stdc++.h>
using namespace std;
const int mod = 20101009;
typedef long long ll;

bool not_prime[10000001];
int mu[10000001];
int prime[1000001];
void mus(int n){
    mu[1] = 1;int tot=0;
    for(int i=2;i<=n;i++) {
        if(not_prime[i]==0) prime[++tot] = i,mu[i]=-1;
        for(int j=1;j<=tot&& prime[j]*i<=n;j++) {
            not_prime[prime[j]*i] = 1;
            if(i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]] = -mu[i];
        }
    }

    for(int i=1;i<=n;i++) {
        mu[i] = mu[i] *1ll* i % mod * i % mod;
        mu[i] += mu[i-1];
        mu[i] %= mod;
    }
}

int sum2(int n,int m){
    int a = (n*1ll*(n+1)/2)%mod;
    int b = (m*1ll*(m+1)/2)%mod;
    return a*1ll*b%mod;
}

int sum(int n,int m) {

    if(n>m) swap(n,m);
    int res = 0;
    for(int d=1;d<=n;) {
        int en = min(n/(n/d),m/(m/d));
        int a = sum2(n/d,m/d);
        int b = (mu[en] - mu[d-1] + mod)%mod;
        res += (a*1ll*b%mod);
        res %=mod;
        d = en+1;
    }
    return res;
}
int solve(int n,int m) {
    if(n>m) swap(n,m);
    int res = 0;
    for(int d=1;d<=n;) {
        int en = min(n/(n/d),m/(m/d));
        int a = sum(n/d,m/d);

        int b = ((d+en)*1ll*(en-d+1)/2)%mod;
        res += (a*1ll*b%mod);
        res %= mod;
        d = en+1;
    }
    return res;
}

int main() {
    int n,m;
    cin>>n>>m;
    mus(n);
    cout<<solve(n,m)<<endl;
    return 0;
    int res =0;
    for(int i=1;i<=n;i++) {
        for(int j=1;j<=m;j++) {
            int g = __gcd(i,j);
            res += ((i*1ll*j / g) %mod);
            res %=mod;
        }
    }
    cout<<res<<endl;
}