luogu P1587 [NOI2016] 循环之美

https://www.luogu.com.cn/problem/P1587

首先思考我们要求的是什么？

\(\large \begin{Bmatrix}\frac{xk^l}{y}\end{Bmatrix} =\begin{Bmatrix}\frac{x}{y}\end{Bmatrix}\)

\(\large xk^l-y \begin{bmatrix}\frac{xk^l}{y} \end{bmatrix}=x-y \begin{bmatrix}\frac{x}{y} \end{bmatrix}\)
\(k^l \equiv 1 (\mod y)\)

可以得到我们实际上要求的是

\[\sum_{i=1}^n\sum_{j=1}^m [i\perp j] [j \perp k] \]

按照套路

\(\large \sum\limits_{d=1}^{min(n, m)} \mu(d) \sum\limits_{i=1}^{n/d} \sum\limits_{j=1}^{m/d} [jd \perp k]\)

\(\large = \sum\limits_{d=1}^{min(n, m)} \mu(d) (n/d) \sum\limits_{j=1}^{m/d} [j \perp k][d \perp k]\)

设\(g(n)=\sum\limits_{i=1}^n [i\perp k]\)

容易发现\(g(n)=(n/k) \phi(k)+g(n \mod k)\)
上式子

\(\large = \sum\limits_{d=1}^{min(n, m)} (n/d) g(m/d) \mu(d) [d \perp k]\)
前两项显然可以整除分块，现在要计算的是后两项的区间和

记\(f(n)=\sum\limits_{i=1}^n \mu(i)[i \perp k]\)

区间和就可以表示为\(f(r)-f(l-1)\)

考虑怎么快速计算\(f\)

发现是两个积性函数卷在一起,还是一个积性函数

我们考虑让它卷上\([i \perp k]\)

得到杜教筛的那条式子

\(\large f(n)=\sum\limits_{i=1}^n[i \perp k] f(n/i)-\sum\limits_{i=2}^n[i\perp k] f(n/i)\)

考虑前面那一部分怎么算？

\(\sum\limits_{i=1}^n[i \perp k] \sum\limits_{j=1}^{n/i}\mu(j)[j\perp k]\)
\(=\sum\limits_{i=1}^n\sum\limits_{j=1}^{n/i}\mu(j)[ij\perp k]\)

考虑枚举\(ij\)

\(\sum\limits_{i=1}^n \sum\limits_{d|i}\mu(d)[i\perp k]\)
\(\sum\limits_{i=1}^n[i\perp k] \sum\limits_{d|i} \mu(d)=1\)

于是乎就得到了

\(\large f(n)=1-\sum\limits_{i=2}^n[i\perp k] f(n/i)\)

而\(f(n/i)\)前面那个在整除分块的时候就是\(g\)
大力杜教筛即可

code:

#include<bits/stdc++.h>
#define ll long long
#define N 5000050
using namespace std;
int gcd(int x, int y) {
    return y? gcd(y, x % y) : x;
}
int g[N], f[N], mu[N], n, m, k, prime[N], sz, vis[N];
ll G(int n) {
    return 1ll * (n / k) * g[k] + g[n % k];
} 
const int lim = 5000000;
void init(int n) {
    for(int i = 1; i <= k; i ++) g[i] = g[i - 1] + (gcd(i, k) == 1);
    f[1] = mu[1] = 1;
    for(int i = 2; i <= n; i ++) {
        if(!vis[i]) {
            prime[++ sz] = i;
            mu[i] = -1;
        }
        for(int j = 1; j <= sz && i * prime[j] <= n; j ++) {
            vis[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
            mu[i * prime[j]] = - mu[i];
        }
        f[i] = f[i - 1] + mu[i] * (G(i) - G(i - 1));
    }
}
map<int, ll> sf;
ll F(int n) {
    if(n <= lim) return f[n];
    if(sf[n]) return sf[n];
    ll ret = 1;
    for(int l = 2, r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        ret -= F(n / l) * (G(r) - G(l - 1)); 
    } return sf[n] = ret;
}
int main() {
    scanf("%d%d%d", &n, &m, &k);
    init(lim);
    ll ans = 0;
    for(int l = 1, r = 0; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += 1ll * (n / l) * G(m / l) * (F(r) - F(l - 1));
    }
    printf("%lld", ans);
    return 0;
}