杜教筛

杜教筛可以在 \(O(n^{3/4})\) 或预处理 \(O(n^{2/3})\) 的时间复杂度内求出某些积性函数前 \(n\) 项的和。

考虑积性函数 \(f(x)\)，我们要求出 \(s(x)=\sum_{i=1}^{x}f(i)\) 在 \(n\) 的取值 \(s(n)\)。考虑构造另一积性函数 \(g(x)\) 满足 \(g(x)\) 和 \((f*g)(x)\) 都可以快速求出前缀和。那么：

\[ \begin{align*} \sum_{i=1}^{n}(f*g)(i)&=\sum_{i=1}^{n}\sum_{j|i}f(j)g\Bigl(\frac ij\Big)\\ &=\sum_{i=1}^{n}\sum_{j=1}^{\lfloor n/i\rfloor}g(i)f(j)\\ &=\sum_{i=1}^{n}g(i)s\Bigl(\left\lfloor\frac ni\right\rfloor\Bigr) \end{align*} \]

移项：

\[ s(n)=\sum_{i=1}^{n}(f*g)(i)-\sum_{i=2}^{n}g(i)s\Bigl(\left\lfloor\frac ni\right\rfloor\Bigr) \]

后面一项用整除分块解决，并且本质不同的 \(s\bigl(\left\lfloor\frac ni\right\rfloor\bigr)\) 只有 \(O(\sqrt n)\) 种，用记忆化即可。算上转移，时间复杂度为 \(O(n^{3/4})\)。

证明

考虑递归过程，\(\big\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\big\rfloor=\big\lfloor\frac{n}{ab}\big\rfloor\)，因此不管递归几层，总有 \(n\in\big\{\lfloor\frac{n}{x}\rfloor\mid x\in\N^*\wedge x\le n\big\}\)。因此状态数是 \(O(\sqrt n)\) 的。

时间复杂度可以写成

\[ O\Bigl(\sum_{i\in\{\lfloor\frac{n}{x}\rfloor\}}\sqrt i~\Bigr)=O\Big(\sum_{i=1}^{\sqrt n}\sqrt\frac{n}{i}\Big)=O(\sqrt n\int_1^{\sqrt n}x^{-\frac 12}dx)=O(n^{3/4}) \]

提前筛出 \(s\) 的前 \(n^{2/3}\) 项，可以将时间复杂度降低到 \(O(n^{2/3})\)。具体分析可以看四月樱花的分析，和这个差不多。

\(f(x)\)	\(g(x)\)	\((f*g)(x)\)
\(\mu(x)\)	\(1(x)\)	\(\epsilon(x)\)
\(\varphi(x)\)	\(1(x)\)	\(\operatorname{id}(x)\)
\(\sigma_k(x)\)	\(\mu(x)\)	\(\operatorname{id}^k(x)\)
\(x\times f(x)\)	\(x\)	\(x\times (f*1)(x)\)
\(x\times f(x)\)	\(x\times g(x)\)	\(x\times (f*g)(x)\)

P13527 [OOI 2023] Another n-dimensional chocolate bar / n 维巧克力问题

时间复杂度分析和杜教筛很像。设 \(f_{i,j}\) 表示考虑了前 \(i\) 个维度，剩下的维度还需要切出 \(j\) 块，\(\cdots\) 的最大值。\(j\) 只有 \(\sqrt n\) 种取值，转移也只会从 \(\sqrt j\) 个位置转移，因此时间复杂度和上面一样，是 \(O(nk^{3/4})\)。

P3768 简单的数学题

给定 \(n\)，求

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

\(n\le 10^{10}\)

枚举 \(\gcd\)：

\[ \sum_{d=1}^{n}d\sum_{d|i}\sum_{d|j}ij[\gcd(i,j)=1]=\sum_{d=1}^{n}d^3\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ij[\gcd(i,j)=1] \]

转化成欧拉函数（互质数求和结果为 \(\frac{n}{2}\varphi(n)\)，当 \(n>1\)；证明考虑将互质的数字两两配对）：

\[ \sum_{d=1}^{n}d^3\big((-1)+\sum_{i=1}^{\lfloor\frac nd\rfloor}i^2\varphi(i)\big) \]

后面的东西是积性函数，并且 \(\lfloor\frac nd\rfloor\) 的取值集合和杜教筛的重合，因此直接杜教筛即可。具体的，构造函数 \(g(x)=x^2\)，记 \(f(x)=x^2\varphi(x)\)，那么 \((f*g)(x)=x^3\)。

是否用哈希表对常数影响不大。为了更加通用（二维整除分块），可以用哈希表。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 3e6 + 10;
int MOD, inv_4, inv_6;

inline int qpow(int a, int b) { int res = 1; while(b) { (b & 1) && (res = res * a % MOD); a = a * a % MOD; b >>= 1; } return res; }

int n;

bool npri[N]; int pri[N], pcnt;
int phi[N];

struct hash_table {
    static const int MOD = 200029;
    struct Node {
        int key, v, next;
    } pool[2 * N];
    int nn, head[MOD];
    void insert(int key, int v) {
        pool[++nn] = {key, v, head[key % MOD]};
        head[key % MOD] = nn;
    }
    int query(int key) {
        for(int i = head[key % MOD]; i; i = pool[i].next) {
            if(pool[i].key == key) return pool[i].v;
        }
        return -1;
    }
} mp;

inline int sum3(int n) { return n * n % MOD * (n + 1) % MOD * (n + 1) % MOD * inv_4 % MOD; }
inline int sum3(int l, int r) { return (sum3(r % MOD) - sum3((l - 1) % MOD) + MOD) % MOD; }
inline int sum2(int n) { return n * (n + 1) % MOD * (2 * n + 1) % MOD * inv_6 % MOD; }
inline int sum2(int l, int r) { return (sum2(r % MOD) - sum2((l - 1) % MOD) + MOD) % MOD; }

int calc(int n) {
    if(n < N) return phi[n];
    if(~mp.query(n)) return mp.query(n);
    int res = sum3(n % MOD);
    for(int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res = (res + MOD - sum2(l, r) * calc(n / l) % MOD) % MOD;
    }
    mp.insert(n, res);
    return res;
}

signed main() {

    cin >> MOD >> n;
    inv_4 = qpow(4, MOD - 2); inv_6 = qpow(6, MOD - 2);

    phi[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!npri[i]) pri[++pcnt] = i, phi[i] = i - 1;
        for(int j = 1; j <= pcnt; j++) {
            if(i * pri[j] >= N) break;
            npri[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            } else phi[i * pri[j]] = phi[i] * (pri[j] - 1);
        }
    }

    for(int i = 1; i < N; i++) phi[i] = (phi[i - 1] + phi[i] * (i * i % MOD) % MOD) % MOD;

    int ans = 0;
    for(int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans = (ans + sum3(l, r) * calc(n / l) % MOD) % MOD;
    }
    cout << ans << '\n';

    return 0;
}

P1587 [NOI2016] 循环之美

给定 \(n,m,k\)，请你统计 \(k\) 进制下

\[ \Bigl\{\frac{a}{b}\mid 1\le a\le n\wedge 1\le b\le m\Big\} \]

中有多少个数字是纯循环或者是整数。

\(1\le n,m\le 10^9,~2\le k\le 2000\)

考察有理数 \(\frac{a}{b}\) 在 \(k\) 进制下纯循环的条件：假设（最小）循环节长度为 \(p\)，那么

\[ \Bigl\{\frac{a}{b}k^p\Big\}=\Bigl\{\frac{a}{b}\Big\} \]

写成下取整函数也就是

\[ \frac{a}{b}k^p-\Bigl\lfloor\frac{a}{b}k^p\Big\rfloor=\frac{a}{b}-\Bigl\lfloor\frac{a}{b}\Big\rfloor \]

发现形式很像取模，两遍都乘以 \(b\)：

\[ \begin{align*} ak^p-b\Bigl\lfloor\frac{ak^p}{b}\Big\rfloor&=a-b\Bigl\lfloor\frac{a}{b}\Big\rfloor\\ ak^p\bmod b&=a\bmod b\\ ak^p&\equiv a\pmod b\\ a(k^p-1)&\equiv 0\pmod b \end{align*} \]

由于值相同的分数只统计一遍，因此钦定 \(a,b\) 互质即可，此时 \(a\) 在模 \(b\) 意义下存在逆元，因此式子等价于

\[ k^p\equiv 1\pmod b \]

因此 \(k\perp b\)。答案可以写成

\[ \sum_{i=1}^{n}\sum_{j=1}^{m}\bigl[(i,j)=1\wedge (j,k)=1\big] \]

考虑到 \([(x,k)=1]\) 是完全积性函数，因此对 \((i,j)=1\) 进行莫反：

\[ \sum_{d=1}\mu(d)\bigl[(d,k)=1\big]\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\bigl[(j,k)=1\big]=\sum_{d=1}\mu(d)\bigl[(d,k)=1\big]\Big\lfloor\frac{n}{d}\Big\rfloor f\Bigl(\Bigl\lfloor\frac{m}{d}\Big\rfloor\Big) \]

注意到 \(f(n)=\sum_{i=1}^{n}[(i,k)=1]=\varphi(k)\lfloor\frac{n}{k}\rfloor+f(n\bmod k)\)，可以 \(O(k)-O(1)\) 求出。问题转化为求 \(\mu(d) [(d,k)=1]\) 的前缀和。将它卷 \([(x,k)=1]\)，并求前缀和：

\[ \sum_{ij=1}^{n}\bigl[(ij,k)=1\big]\sum_{i|ij}\mu(i)=1 \]

杜教筛即可。

代码

#include<iostream>
#include<map>
#define int long long
#define ll long long
using namespace std;
const int N = 1000010;

int n, m, k;
ll ans;

int npri[N], pri[N], pcnt;
int mu[N], smu[N], perp[N];

int f(int n) {
    return perp[k] * (n / k) + perp[n % k];
}

map<int, int> mp;
int get_smu(int k) {
    if(k < N) return smu[k];
    if(mp.count(k)) return mp[k];
    int &res = mp[k]; res = 1;
    for(int l = 2, r; l <= k; l = r + 1) {
        r = k / (k / l);
        res -= (f(r) - f(l - 1)) * get_smu(k / l);
    }
    return res;
}

signed main() {

    cin >> n >> m >> k;

    mu[1] = perp[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!npri[i]) {
            pri[++pcnt] = i;
            mu[i] = -1; perp[i] = (bool)(k % i);
        }
        for(int j = 1; j <= pcnt; j++) {
            if(i * pri[j] >= N) break;
            npri[i * pri[j]] = 1;
            perp[i * pri[j]] = perp[i] * perp[pri[j]];
            if(i % pri[j] == 0) {
                mu[i * pri[j]] = 0;
                break;
            } else {
                mu[i * pri[j]] = -mu[i];
            }
        }
    }

    for(int i = 1; i < N; i++) smu[i] = smu[i - 1] + mu[i] * perp[i];
    for(int i = 1; i < N; i++) perp[i] += perp[i - 1];

    for(int l = 1, r; l <= n && l <= m; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += (get_smu(r) - get_smu(l - 1)) * ((ll)n / l * f(m / l));
    }

    cout << ans << '\n';

    return 0;
}