[数论记录]P1587 [NOI2016] 循环之美

学了两个方法，好强！

题意：

求 \(\displaystyle \sum\limits_{i=1}^n \sum_{j=1}^{m} [i \perp j] [j \perp k]\)

\(n,m \leq 10^9,k \leq 2000\)

\[\newcommand{\frc}[2] { \lfloor \frac{#1}{#2} \rfloor } \begin{aligned} & \sum\limits_{i=1}^n \sum_{j=1}^{m} [i \perp j] [j \perp k] \\ = & \sum_{j=1}^m [j \perp k] \sum_{d|i,d|j} \mu(d)\\ = & \sum_{d=1}^n \frc{n}{d} [d \perp k] \mu(d) \cdot \sum_{j=1}^{\frc{n}{d}} [j \perp k]\\ = & \sum_{d=1}^n f(d,k) \cdot \frc{n}{d} \cdot g(\frc{m}{d},k) \\ \end{aligned} \]

其中

\[\begin{aligned} f(n,k) =& \sum_{i=1}^n [i \perp k] \mu(i) \\ g(n,k) =& \sum_{i=1}^n [i \perp k] \end{aligned} \]

要弄出 \(f\) 的前缀和和 \(g\) 的特殊点值。

\(g\) 是非常轻易的，现在来考虑 \(f\) 的前缀和。

两个思路。

其一，考虑处理 \(k\)。

不妨假设 \(k\) 没有平方因子。

记 \(k' = \dfrac{k}{p}\)，\(p\) 是 \(k\) 最小因子。

\[\begin{aligned} f(n,k) &= \sum_{i=1}^n [i \perp k] \mu(i) \\ &= \sum_{i=1}^n [i \perp k'] \mu(i) - \sum_{i=1}^n [p|i] [i \perp k'] \mu(i) \\ &= f(n,k') - \sum_{i=1}^{\frc{n}{p}} \mu(ip) [ip \perp k'] \\ &= f(n, k') - \sum_{i=1}^{\frc{n}{p}} \mu(i) \mu(p) [i \perp p] [i \perp k'] [p \perp k'] \\ &= f(n, k') + \sum_{i=1}^{\frc{n}{p}} \mu(i) [i \perp k] \\ &= f(n, k') + f(\frc{n}{p},k) \end{aligned} \]

边界 \(k=1\) 和 \(n=0\)，容易通过杜教筛处理。预处理或记忆化均可，后者感觉好写些。

其二，考虑继续反演

\[\begin{aligned} f(n,k) &= \sum_{i=1}^n [i \perp k] \mu(i) \\ &= \sum_{i=1}^n \mu(i) \sum_{d|i,d|k} \mu(d) \\ &= \sum_{d|k} \mu(d) \sum_{i=1}^{\frc{n}{p}} \mu(id) \\ &= \sum_{d|k} \mu^2(d) \sum_{i=1}^{\frc{n}{p}} \mu(i) [i \perp d] \\ &= \sum_{d|k} \mu^2(d) f(\frc{n}{d},d) \end{aligned} \]

很漂亮！

直接做极不好做，尤其是我第一步想反演 \([y \perp k]\)。我想做法归功于递归。

#include <cstdio>
#include <map>
#include <iostream>
#include <algorithm>
#include <cmath>
#define LL long long
using namespace std;
const int M = 1e5 + 5, N = 1e6 + 5;
int s[M], fac[M], cnt;
map<pair<int, int>, LL> tf;
int summu[N], mu[N], p[N], tt; bool f[N];
void pre(int n){
    summu[1] = f[1] = mu[1] = 1;
    for(int i = 2; i <= n; i++){
        if(!f[i]) p[++tt] = i, mu[i] = -1;
        for(int j = 1; j <= tt; j++){
            int k = p[j] * i; if(k > n) break;
            if(i % p[j] == 0) {f[k] = 1; mu[k] = 0; break;}
            f[k] = 1; mu[k] = -mu[i];
        }
        summu[i] = summu[i-1] + 1ll * mu[i];
    } 
}
int T, mx = 1000000; map<int, LL> smu;
LL Summu(int n){
    if(n <= mx) return summu[n];
    if(smu[n]) return smu[n];
    LL ans = 1ll;
    for(LL l = 2, r; l <= n; l = r+1){
        r = 1ll * n/(n/l); ans -= (r-l+1) * Summu(n/l);
    } return smu[n] = ans;
}
LL F(int n, int k) {
    if(k == 1) return Summu(n);
    if(!n) return 0;
    if(tf.find(make_pair(n, k)) != tf.end()) return tf[make_pair(n, k)];
    LL ans = 0;
    for(int i = 1; i*i <= k; i++) {
        if(k % i == 0) {
            if(mu[i]) ans += F(n/i, i);
            if(i*i != k && mu[k/i]) ans += F(n / (k / i), k / i);
        }
    }
    return tf[make_pair(n, k)] = ans;
}
LL G(int n, int k) {
    return (n / k) * s[k] + s[n % k];
}
LL solve(int n, int m, int k) {
    LL ans = 0;
    for(int l = 1, r; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += (F(r, k) - F(l-1, k)) * (n / l) * G(m / l, k);
    }
    return ans;
}
int main(){
    pre(1e6); int n, m, k;
    scanf("%d %d %d", &n, &m, &k);
    for(int i = 1; i <= k; i++){
        s[i] = s[i-1] + (__gcd(i, k) == 1);
        if(k % i == 0 && mu[i] != 0) fac[++cnt] = i;
    }
    printf("%lld\n", solve(n, m, k));
}