bzoj2693 - jzptab（莫比乌斯函数，数论分块）

题目

求\(\sum\limits^{n}_{i=1}\sum\limits^{m}_{j=1}{lcm(i, j)}\)
答案模100000009输出

题解

注意100000009是1e8级别，而且不是质数，不要想逆元之类的东西。

老老实实推式子，一步一步小心不要出错。注意在提出项时可能要平方之类的。

\[\sum\limits_{d_1}\sum\limits_{i=1}^n\sum\limits_{j=1}^m{\frac{ij}{d_1}[\gcd(i,j)=d_1]} \]

提出\(d_1\)

\[\sum\limits_{d_1}\sum\limits_{i=1}^{\frac{n}{d_1}}\sum\limits_{j=1}^{\frac{m}{d_1}}{ijd_1[\gcd(i,j)=1]} \]

\[\sum\limits_{d_1}\sum\limits_{i=1}^{\frac{n}{d_1}}\sum\limits_{j=1}^{\frac{m}{d_1}}{ijd_1\sum\limits_{d_2|\gcd(i,j)}{\mu(d_2)}} \]

\[\sum\limits_{d_1}\sum\limits_{d_2}{\mu(d_2)d_1}\sum\limits_{i=1}^{\frac{n}{d_1}}{i[d_2|i]}\sum\limits_{j=1}^{\frac{m}{d_1}}{j[d_2|j]} \]

提出\(d_2\)

\[\sum\limits_{d_1}\sum\limits_{d_2}{\mu(d_2)d_1d_2^2}\sum\limits_{i=1}^{\frac{n}{d_1d_2}}{i}\sum\limits_{j=1}^{\frac{m}{d_1d_2}}{j} \]

令\(d_1d_2=T\)

\[\sum\limits_{T=1}^{\min(n,m)}(\sum\limits_{d|T}{\mu(d)dT})\sum\limits_{i=1}^{\frac{n}{T}}{i}\sum\limits_{j=1}^{\frac{m}{T}}{j} \]

令\(F(T)=\sum\limits_{d|T}{\mu(d)dT}\)，求其前缀和，就可以数论分块了。

#include <bits/stdc++.h>

#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N) 
typedef long long ll;

using namespace std;
/*-----------------------------------------------------------------*/

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f

const int N = 1e7 + 10;
const double eps = 1e-5;
const int M = 1e8 + 9;

int mu[N];
ll f[N];
int prime[N];
int cnt;
int vis[N];

void init() {
    mu[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!vis[i]) {
            mu[i] = -1;
            prime[cnt++] = i;
        }
        for(int j = 0; j < cnt; j++) {
            int p = prime[j];
            if(i * p > N) break;
            vis[i * p] = 1;
            if(i % p == 0) {
                mu[i * p] = 0;
                break;
            }
            mu[i * p] = -mu[i];
        }
    }
    for(int i = 1; i < N; i++) {
        for(int j = i; j < N; j += i) {
            f[j] += 1ll * mu[i] * i % M;
            f[j] %= M;
        }
    }
    for(int i = 2; i < N; i++) {
        f[i] = f[i] * i % M;
        f[i] += f[i - 1];
        f[i] %= M;
    }
}

inline ll sum(int n) {
    return 1ll * (1 + n) * n / 2 % M;
}




int main() {
    IOS;
    init();
    int t;
    cin >> t;
    while(t--) {
        ll ans = 0;
        int n, m;
        cin >> n >> m;
        int cur = 1;
        while(cur <= min(n, m)) {
            int p1 = n / (n / cur);
            int p2 = m / (m / cur);
            int nt = min(min(n, m), min(p1, p2));
            ans += 1ll * sum(n / cur) * sum(m / cur) % M * (f[nt] - f[cur - 1]) % M;
            ans %= M;
            cur = nt + 1;
        }
        cout << (ans + M) % M << endl;
    }
}