【BZOJ2154】Crash的数字表格

【BZOJ2154】Crash的数字表格

题面

bzoj

洛谷

题解

不妨设(nleq m)

题目是求：

[sum_{i=1}^nsum_{j=1}^mlcm(i,j) ]

还是照常推式子qaq

[sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)}\ sum_{d=1}^nsum_{i=1}^nsum_{j=1}^m[gcd(i,j)==d]frac{ij}{d}\ sum_{d=1}^nsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ijd\ sum_{d=1}^ndsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ij ]

暂时不管前面的，先搞后面的一坨

设(x=n/d,y=m/d)

并且设

[f(k)=sum_{i=1}^xsum_{j=1}^y[gcd(i,j)==k]ij\ g(k)=sum_{k|d}f(d) ]

常见套路

[g(k)=sum_{i=1}^xsum_{j=1}^y[k|gcd(i,j)]ij ]

提出(k)

[g(k)=k^2sum_{i=1}^{x/k}sum_{j=1}^{y/k}[1|gcd(i,j)]ij\ g(k)=k^2sum_{i=1}^{x/k}sum_{j=1}^{y/k}ij ]

这个就是两个等差数列相乘。

要求的东西变为

[f(1)=sum_{i=1}^xmu(i)g(i) ]

问题解决很多了，但是还必须要往下

[ans=sum_{d=1}^ndsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ij ]

看看(f(1))是什么

[f(1)=sum_{i=1}^xmu(i)i^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}ij ]

里外都可以数论分块，总复杂度(O(n))

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring> 
#include <cmath> 
#include <algorithm> 
using namespace std;
const int Mod = 20101009;
const int MAX_N = 1e7 + 5; 
int N = 1e7, M, mu[MAX_N], prime[MAX_N], num; 
bool is_prime[MAX_N];
int s[MAX_N]; 
void sieve() {
    for (int i = 1; i <= N; i++) is_prime[i] = 1; 
    is_prime[1] = 0, mu[1] = 1; 
    for (int i = 2; i <= N; i++) {
        if (is_prime[i]) { prime[++num] = i, mu[i] = -1; }
        for (int j = 1; j <= num && i * prime[j] <= N; j++) { 
            is_prime[i * prime[j]] = 0; 
            if (i % prime[j]) mu[i * prime[j]] = -mu[i]; 
            else { mu[i * prime[j]] = 0; break; } 
        } 
    } 
}
int solve(int x, int y) { 
    int ans = 0;
    for (int l = 1, r; l <= x; l = r + 1) { 
        r = min(x / (x / l), y / (y / l)); 
        int res = 1ll * (1ll * (1 + x / l) * (x / l) / 2 % Mod) * (1ll * (1 + y / l) * (y / l) / 2 % Mod) % Mod;
        ans = (1ll * (s[r] - s[l - 1]) % Mod * res % Mod + ans) % Mod; 
    }
    return (ans + Mod) % Mod; 
} 
int main () {
    sieve(); 
    cin >> N >> M;
    if (N > M) swap(N, M); 
    for (int i = 1; i <= N; i++) s[i] = (s[i - 1] + 1ll * i * i % Mod * mu[i] % Mod + Mod) % Mod;
    int ans = 0; 
    for (int l = 1, r; l <= N; l = r + 1) { 
        r = min(N / (N / l), M / (M / l)); 
        int res = 1ll * (l + r) * (r - l + 1) / 2 % Mod; 
        ans = (ans + 1ll * solve(N / l, M / l) * res % Mod) % Mod; 
    }
    printf("%d
", ans); 
    return 0; 
}