数论函数相关

数论函数相关

数论函数

定义：

  • 数论函数：定义域为正整数的函数称为数论函数。
  • 加性函数：满足 \(\forall a, b \in \mathbb{N^*}, a \bot b, f(ab) = f(a) + f(b)\) 的函数。
  • 积性函数：满足 \(\forall a, b \in \mathbb{N^*}, a \bot b, f(ab) = f(a) \times f(b)\) 的函数。
  • 完全积性函数：满足 \(\forall a, b \in \mathbb{N^*}, f(ab) = f(a) \times f(b)\) 的函数。

数论函数的运算：

  • 加法：\((f + g)(x) = f(x) + g(x)\) 。
  • 数乘：\((c \cdot f)(x) = c \times f(x)\) 。
  • 点乘：\((f \cdot g)(x) = f(x) \times g(x)\) 。

常见数论函数（记 \(n\) 的唯一分解为 \(\prod_{i = 1}^m p_i^{c_i}\) ）：

  • 单位函数： \(\epsilon(n) = [n - 1]\) ，为完全积性函数。

  • 常数函数：\(1(n) = 1\) ，为完全积性函数。

  • 幂函数：\(\mathrm{id}_k(n) = n^k\) ，为完全积性函数。

  • 除数函数：\(\sigma_k(n) = \sum_{d \mid n} d^k\) ，为积性函数。

    • \(\sigma_0(n) = d(n) = \prod_{i = 1}^m (c_i + 1)\) 表示 \(n\) 的约数个数。

      • \(d(ij) = \sum_{x | i} \sum_{y | j} [\gcd(x, y) = 1]\) 。

        证明：考虑一个质因子，若都没被枚举则不统计，否则若被 \(x\) 统计，则贡献为 \(x\) 里对应的次幂，否则贡献为 \(i\) 里对应的次幂加上 \(j\) 里对应的次幂。

      • \(d(ijk) = \sum_{a | i} \sum_{b | j} \sum_{c | k} [\gcd(a, b) = 1] [\gcd(b, c) = 1] [\gcd(a, c) = 1]\) 。

      • \(\sigma_0(n^2) = \sum_{d \mid n} \mu(d) \sigma_0(\frac{n}{d})^2 = \sum_{d \mid n} \mu(d)^2 \sigma_0(\frac{n}{d})\) 。

        证明：只需证明 \(n = p^k\) 成立即可，原式等价于 \(2k + 1 = (k + 1)^2 - k^2 = k + (k + 1)\) 。

    • \(\sigma_1(n) = \sigma(n) = \prod_{i = 1}^m \frac{p_i^{c_i + 1} - 1}{p_i - 1}\) 表示 \(n\) 的约数和，具有如下性质

      • \(\sum_{i = 1}^n \sigma(i) = \sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor\) 。

        证明：考虑每个 \(i\) 的贡献即可。

    • \(k \ge 1\) 时 \(\sigma_k(n) = \prod_{i = 1}^n \frac{p_i^{(c_i + 1)k} - 1}{p_i - 1}\) ，这是因为 \(n\) 的所有约数的 \(k\) 次方和可以写作 \(\prod_{i = 1}^m \sum_{j = 0}^{c_i} p_i^{jk}\) ，等比数列求和即可。

  • 本质不同质因子个数函数：\(\omega(n) = \sum_{p \in \mathbb{P}} [p \mid n]\) ，为加性函数。

线性筛

基础的用法是筛素数：

inline void sieve(int n) {
    memset(isp, true, sizeof(isp));
    isp[1] = false;
    
    for (int i = 2; i <= n; ++i) {
        if (isp[i])
            pri[++pcnt] = i;
        
        for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
            isp[i * pri[j]] = false;
            
            if (!(i % pri[j]))
                break;
        }
    }
}

线性筛可以 \(O(n)\) 求出某些积性函数在 \(1 \sim n\) 处的所有取值，需要满足 \(f(p^c)\) 的值便于分析。

线性筛时若找到一个质数 \(p\) ，则将所有 \(p\) 的幂处的函数取值都求出。记 \(low_i\) 最小质因子的次幂，则对于合数 \(i\) 有 \(f(i) = f(\frac{i}{low_i}) \times f(low_i)\) 。

inline void sieve(int n) {
    memset(isp, true, sizeof(isp));
    isp[1] = false, f[1] = 1;
    
    for (int i = 2; i <= n; ++i) {
        if (isp[i])
            pri[++pcnt] = i, low[i] = i, f[i] = ...;

        for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j]) {
                low[i * pri[j]] = low[i] * pri[j];

                if (i == low[i])
                    f[i * pri[j]] = ...;
                else
                    f[i * pri[j]] = f[i / low[i]] * f[low[i * pri[j]]];
            } else {
                low[i * pri[j]] = low[j], f[i * pri[j]] = f[i] * f[pri[j]];
                break;
            }
        }
    }
}

Dirichlet 卷积

定义狄利克雷卷积：

\[h(n) = \sum_{d \mid n} f(d) g(\frac{n}{d}) \]

简记为 \(h = f * g\) ，按定义式计算可以做到 \(O(n \ln n)\) 。

由于 \(\epsilon * f = f\) ，因此 \(\epsilon\) 称为狄利克雷卷积的单位元。

定义 \(f\) 的逆元 \(f^{-1}\) 满足 \(f *f^{-1} = \epsilon\) ，其中 \(f\) 可逆当且仅当 \(f(1) \ne 0\) ，可以计算得出：

\[f^{-1}(n) = \begin{cases} \frac{1}{f(1)} & n = 1 \\ - \frac{\sum_{d \mid n, d \ne n} f^{-1}(d) f(\frac{n}{d})}{f(1)} & n \ge 1 \end{cases} \]

性质：

  • 交换律：\((f * g)(n) = (g * f)(n)\) 。
  • 结合律：\(((f * g) * h)(n) = (f * (g * h))(n)\) 。
  • 分配律：\(((f + g) * h)(n) = (f * h + g * h)(n)\) 。
  • \(f = g \iff f * h = g * h (h(1) \ne 0)\) 。
  • 积性函数的狄利克雷卷积是积性函数。
  • 积性函数的逆元是积性函数。

常见的狄利克雷卷积式：

  • \(\sigma_0 = 1 * 1\) 。
  • \(\mathrm{id} = \varphi * 1\) 。
  • \(\sigma = 1 * \mathrm{id} = \sigma_0 * \varphi\) 。
  • \(1 = \sigma_0 * \mu\) 。
  • \(\epsilon = \mu * 1\) 。
  • \(\varphi = \mu * \mathrm{id}\) 。

整除分块

注意到 \(\lfloor \frac{n}{x} \rfloor\) 在 \(x \in [1, n]\) 时的取值不超过 \(2 \sqrt{n}\) 种，对于某些问题，对于一段 \(\lfloor \frac{n}{x} \rfloor\) 均相等的区间，可以直接得出该区间里的所有 \(x\) 对答案的贡献。

如果要实现整块一起统计，则需要求出每一块的块头 \(l\) 和块尾 \(r\)，得到答案为：

\[\sum_{i=1}^{k} \lfloor \frac{n}{i} \rfloor = \sum_{[l,r]} (r-l+1)(\lfloor \frac{n}{l} \rfloor) \]

注意到每一块的 \(l\) 都可以由上一块的 \(r\) 推出，故我们继续讨论如何在已知 \(l\) 的情况下推出 \(r\) 。令 \(t = \lfloor \frac{n}{l} \rfloor\) ，容易得到：

\[\begin{cases} r = \min(\lfloor \frac{n}{t} \rfloor, n) & (t \ne 0) \\ r = n & (t = 0) \end{cases} \]

直接计算即可，时间复杂度 \(O(\sqrt{n})\) 。

for (int l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    // Calculate the contribution of [l, r]
}

扩展到二维甚至多维也是类似的，时间复杂度 \(O(k \sqrt{n})\) ，其中 \(k\) 为维度。

二维整除分块：

for (int l = 1, r; l <= min(n, m); l = r + 1) {
    r = min(n / (n / l), m / (m / l));
    // Calculate the contribution of [l, r]
}

一类题目的整除分块形式为向上取整，只需令 \(r \gets \lfloor \frac{n - 1}{\lceil \frac{n}{l} \rceil - 1} \rfloor\) 即可，这是因为：

\[\frac{n}{r} > \lceil \frac{n}{l} \rceil - 1 \iff r (\lceil \frac{n}{l} \rceil - 1) < n \iff r < \frac{n}{\lceil \frac{n}{l} \rceil - 1} \iff r \le \frac{n - 1}{\lceil \frac{n}{l} \rceil - 1} \]

注意特判 \(\lceil \frac{n}{l} \rceil = 1\) ，此时 \(r\) 的上界为 \(+ \infty\) ，需要取实际上界。

欧拉函数

欧拉函数 \(\varphi(n)\) 定义为 \(1 \sim n\) 内与 \(n\) 互质的数的个数，具有如下性质：

  • 欧拉函数是积性函数，特别的有 \(\varphi(2n) = \varphi(n)\) 。
  • \(\sum_{d | n} \varphi(d) = n\) 。
    • 证明：设 \(f(x)\) 表示 \(\gcd(k, n) = x (k \in [1, n])\) 的个数，则 \(n = \sum_{d | n} f(d) = \sum_{d | n} \varphi(\frac{n}{d}) = \sum_{d | n} \varphi(d)\) 。
  • 当 \(n\) 为质数时， \(\varphi(n) = n - 1\) 。
  • 当 \(n = p^k (p \in \mathbb{P})\) 时，\(\varphi(n) = p^k - p^{k - 1}\) 。
  • 由唯一分解定理，设 \(n = \prod p_i^{k_i}\) ，则 \(\varphi(n) = n \times \prod \frac{p_i - 1}{p_i}\) 。
  • $n > 2 \iff 2 | \varphi(n) $ ，证明考虑对称性即可。
  • \(\varphi(mn) \varphi(\gcd(m, n)) = \varphi(m) \varphi(n) \gcd(m, n)\) 。

单个数的求解：

inline int euler_phi(int n) {
    int phi = n;
    
    for (int i = 2; i * i <= n; ++i)
        if (!(n % i)) {
            phi = phi / i * (i - 1);
            
            while (!(n % i))
                n /= i;
        }
    
    if (n > 1)
        phi = phi / n * (n - 1);
    
    return phi;
}

线性筛求解：

inline void sieve(int n) {
    memset(isp, true, sizeof(isp));
    isp[1] = false, phi[1] = 1;
    
    for (int i = 2; i <= n; ++i) {
        if (isp[i])
            pri[++pcnt] = i, phi[i] = i - 1;
        
        for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
            isp[i * pri[j]] = false;
            
            if (i % pri[j])
                phi[i * pri[j]] = phi[pri[j]] * phi[i];
            else {
                phi[i * pri[j]] = pri[j] * phi[i];
                break;
            }
        }
    }
}

利用 \(n = \sum_{d | n} \varphi(d)\) 可以对数论式子做一些变换，如：

\[\gcd(i, j) = \sum_{d | \gcd(i, j)} \varphi(d) = \sum_{d | i} \sum_{d | j} \varphi(d) \]

由此有结论：

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

GCDMAT - GCD OF MATRIX

给定 \(n, m\) ，\(T\) 次询问，每次给出 \(i_1, j_1, i_2, j_2\) ，求：

\[\sum_{i = i_1}^{i_2} \sum_{j = j_1}^{j_2} \gcd(i, j) \pmod{10^9 + 7} \]

\(n, m \le 5 \times 10^4\)

前缀和差分转化为求：

\[\begin{align} S(n, m) &= \sum_{i = 1}^n \sum_{j = 1}^m \gcd(i, j) \\ &= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{d | i} \sum_{d | j} \varphi(d) \\ &= \sum_{d = 1}^n \varphi(d) \sum_{i = 1}^n \sum_{j = 1}^m [d | i] [d | j] \\ &= \sum_{d = 1}^n \varphi(d) \lfloor \dfrac{n}{d} \rfloor \lfloor \dfrac{m}{d} \rfloor \end{align} \]

数论分块可以做到 \(O(n) - O(\sqrt{n})\) 。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 1e9 + 7;
const int N = 5e4 + 7;

int pri[N], phi[N], sum[N];
bool isp[N];

int T, n, m, pcnt;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, phi[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, phi[i] = i - 1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[pri[j] * i] = false;

            if (i % pri[j])
                phi[i * pri[j]] = phi[i] * phi[pri[j]];
            else {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = add(sum[i - 1], phi[i]);
}

inline int S(int n, int m) {
    if (n > m)
        swap(n, m);
    
    int ans = 0;

    for (int l = 1, r = 0; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans = add(ans, 1ll * (sum[r] - sum[l - 1]) * (n / l) % Mod * (m / l) % Mod);
    }

    return ans;
}

signed main() {
    sieve();
    scanf("%d%d%d", &T, &n, &m);

    while (T--) {
        int i1, j1, i2, j2;
        scanf("%d%d%d%d", &i1, &j1, &i2, &j2);
        printf("%d\n", add(dec(dec(S(i2, j2), S(i1 - 1, j2)), S(i2, j1 - 1)), S(i1 - 1, j1 - 1)));
    }

    return 0;
}

P5221 Product

求：

\[\prod_{i = 1}^n \prod_{j = 1}^n \dfrac{\operatorname{lcm}(i, j)}{\gcd(i, j)} \pmod{104857601} \]

\(n \le 10^6\)

\[\begin{align} & \prod_{i = 1}^n \prod_{j = 1}^n \dfrac{\operatorname{lcm}(i, j)}{\gcd(i, j)} \\ =& \prod_{i = 1}^n \prod_{j = 1}^n \dfrac{ij}{\gcd(i, j)^2} \\ =& \dfrac{(n!)^{2n}}{(\prod_{i = 1}^n \prod_{j = 1}^n \gcd(i, j))^2} \\ \end{align} \]

单独考虑下面的部分得到 ：

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

单独考虑指数部分得到：

\[\begin{align} & \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} [\gcd(i, j) = 1] \\ =& \left( \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \left(2 \times \sum_{j = 1}^i [\gcd(i, j) = 1] \right) \right) - 1 \\ =& 2 \times \left( \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \varphi(i) \right) - 1 \\ \end{align} \]

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 104857601;
const int N = 1e6 + 7;

int pri[N], phi[N], sum[N];
bool isp[N];

int n, pcnt;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline void sieve(int n) {
    memset(isp, true, sizeof(isp));
    isp[1] = false, phi[1] = 1;

    for (int i = 2; i <= n; ++i) {
        if (isp[i])
            pri[++pcnt] = i, phi[i] = i - 1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                phi[i * pri[j]] = phi[i] * phi[pri[j]];
            else {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
        }
    }

    for (int i = 1; i <= n; ++i)
        sum[i] = (sum[i - 1] + phi[i]) % (Mod - 1);
}

signed main() {
    scanf("%d", &n);
    sieve(n);
    int fac = 1;

    for (int i = 1; i <= n; ++i)
        fac = 1ll * fac * i % Mod;

    int ans = 1ll * mi(fac, n * 2) * fac % Mod * fac % Mod;

    for (int i = 1; i <= n; ++i)
        ans = 1ll * ans * mi(mi(i, 4ll * sum[n / i] % (Mod - 1)), Mod - 2) % Mod;

    printf("%d", ans);
    return 0;
}

莫比乌斯函数

定义莫比乌斯函数：

\[\mu(n) = \begin{cases} 1 & n = 1 \\ 0 & n 含有平方因子 \\ (-1)^k & k 为 n 的本质不同质因子个数 \end{cases} \]

性质：

  • \(\varphi(n) = \sum_{d | n} \mu(d) \times \frac{n}{d}\) 。
  • \(\sum_{d | n} \mu(d) = [n = 1]\) 。
  • \(\mu^2(n) = \sum_{d^2 | n} \mu(d)\) 。
    • \(d = 1\) 贡献为 \(1\) ，剩下贡献用二项式定理推即可。

单个数的求解：

inline int Mu(int n) {
    int mu = 1;
    
    for (int i = 2; i * i <= n; ++i)
        if (!(n % i)) {
            if (!(n % (i * i)))
                return 0;

            mu *= -1, n /= i;
        }
    
    return n > 1 ? -mu : mu;
}

线性筛求解：

inline void sieve(int n) {
    memset(isp, true, sizeof(isp));
    isp[1] = false, mu[1] = 1;
    
    for (int i = 2; i <= n; ++i) {
        if (isp[i])
            pri[++pcnt] = i, mu[i] = -1;
        
        for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
            isp[i * pri[j]] = false;
            
            if (i % pri[j])
                mu[i * pri[j]] = -mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }
}

事实上可以通过魔力筛在 \(O(n \log \log n)\) 的时间复杂度求出任意数论函数与 \(\mu\) 的狄利克雷卷积。

设 \(g_{i, n} = \sum_{d \mid n} f(d) \mu(\frac{n}{d})\)，其中 \(d\) 只含前 \(i\) 种质因子，则：

\[g_{i,n} = \begin{cases} g_{i - 1,n} & p_i \nmid n \\ g_{i - 1,n} - g_{i - 1, \frac{n}{p_i}} & p_i \mid n \end{cases} \]

需要滚动数组优化。

inline void calc(int n) {
    memcpy(g + 1, f + 1, sizeof(int) * n);
    
    for (int i = 1, i <= pcnt; ++i)
        for (int j = n / pri[i]; j; --j)
            g[j * pri[i]] -= g[j];
}

莫比乌斯反演

设 \(f(n), g(n)\) 为两个数论函数，则：

\[f(n) = \sum_{d \mid n} g(d) \iff g(n) = \sum_{d \mid n} f(d) \mu(\frac{n}{d}) \\ f(n) = \sum_{n \mid d} g(d) \iff g(n) = \sum_{n \mid d} f(d) \mu(\frac{d}{n}) \]

本质就是对质因子做二项式反演。

常见的应用：

\[\begin{aligned} &\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = 1] \\ =& \sum_{i = 1}^n \sum_{j = 1}^m \sum_{d \mid i} \sum_{d \mid j} \mu(d) \\ =& \sum_{d = 1}^n \mu(d) \sum_{i = 1}^n \sum_{j = 1}^m [d \mid i] [d \mid j] \\ =& \sum_{d = 1}^n \mu(d) \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \end{aligned} \]

推广：

\[\begin{align} & \sum_{i = 1}^n \sum_{j = 1}^m f(\gcd(i, j)) \\ =& \sum_{d = 1}^n f(d) \sum_{i = 1}^{n} \sum_{j = 1}^{m} [\gcd(i, j) = d] \\ =& \sum_{d = 1}^n f(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{d} \rfloor} [\gcd(i, j) = 1] \\ =& \sum_{d = 1}^n f(d) \sum_{k = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(k) \lfloor \dfrac{n}{d k} \rfloor \lfloor \dfrac{m}{d k} \rfloor \\ \end{align} \]

令 \(d k = T\) ，得到：

\[\sum_{T = 1}^n \lfloor \dfrac{n}{T} \rfloor \lfloor \dfrac{m}{T} \rfloor \sum_{d | T} f(d) \mu(\dfrac{T}{d}) \]

P2257 YY的GCD / PGCD - Primes in GCD Table / P2568 GCD

求：

\[\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) \in \mathbb{P}] \]

\(n, m \le 10^7\)

推式子得到答案为：

\[\sum_{T = 1}^n \lfloor \dfrac{n}{T} \rfloor \lfloor \dfrac{m}{T} \rfloor \sum_{k | T, k \in \mathbb{P}} \mu(\dfrac{T}{k}) \]

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;

ll f[N], sum[N];
int pri[N], mu[N];
bool isp[N];

int T, n, m, pcnt;

inline void sieve() {
    memset(isp, true, sizeof(isp));
    mu[1] = 1, isp[1] = false;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, mu[i] = -1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[pri[j] * i] = false;

            if (i % pri[j])
                mu[i * pri[j]] = -mu[i];
            else
                break;
        }
    }
    
    for (int i = 1; i <= pcnt; ++i)
        for (int j = 1; pri[i] * j < N; ++j)
            f[j * pri[i]] += mu[j];

    for (int i = 1; i < N; ++i)
        sum[i] = sum[i - 1] + f[i];
}

signed main() {
    sieve();
    scanf("%d", &T);

    while (T--) {
        scanf("%d%d", &n, &m);
        
        if (n > m)
            swap(n, m);
        
        ll ans = 0;

        for (int l = 1, r = 0; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
        }

        printf("%lld\n", ans);
    }

    return 0;
}

P4449 于神之怒加强版

求：

\[\sum_{i = 1}^n \sum_{j = 1}^m \gcd(i, j)^k \pmod{10^9 + 7} \]

\(n, m, k \le 5 \times 10^6\)

答案即为：

\[\sum_{T = 1}^n \lfloor \dfrac{n}{T} \rfloor \lfloor \dfrac{m}{T} \rfloor \sum_{d | T} d^k \mu(\dfrac{T}{d}) \]

线性筛 \(g(T) = \sum_{d | T} d^k \mu(\dfrac{T}{d})\) 即可。

#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 5e6 + 7;

int pri[N], low[N], g[N], sum[N];
bool isp[N];

int T, n, m, k, pcnt;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, low[1] = 1, g[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i]) {
            pri[++pcnt] = i, low[i] = i, g[i] = dec(mi(i, k), 1);
            int x = i;

            while (1ll * x * i < N)
                x *= i, low[x] = x, g[x] = dec(mi(x, k), mi(x / i, k));
        }

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            int x = i * pri[j];
            isp[x] = false;
            low[x] = i % pri[j] ? pri[j] : low[i] * pri[j];
            g[x] = 1ll * g[x / low[x]] * g[low[x]] % Mod;

            if (!(i % pri[j]))
                break;
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = add(sum[i - 1], g[i]);
}

signed main() {
    scanf("%d%d", &T, &k);
    sieve();

    while (T--) {
        scanf("%d%d", &n, &m);
        int ans = 0;

        if (n > m)
            swap(n, m);

        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans = add(ans, 1ll * (n / l) * (m / l) % Mod * dec(sum[r], sum[l - 1]) % Mod);
        }

        printf("%d\n", ans);
    }

    return 0;
}

P3327 [SDOI2015] 约数个数和

求：

\[\sum_{i = 1}^n \sum_{j = 1}^m d(ij) \]

\(T, n, m \le 5 \times 10^4\)

\[\begin{aligned} & \sum_{i = 1}^n \sum_{j = 1}^m d(ij) \\ =& \sum_{i = 1}^n \sum_{j = 1}^m \sum_{x | i} \sum_{y | j} [\gcd(x, y) = 1] \\ =& \sum_{x = 1}^n \sum_{y = 1}^m \lfloor \dfrac{n}{x} \rfloor \lfloor \dfrac{m}{y} \rfloor \sum_{d|x, d | y} \mu(d) \\ =& \sum_{d = 1}^n \mu(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{d} \rfloor} \lfloor \dfrac{n}{id} \rfloor \lfloor \dfrac{m}{jd} \rfloor \\ =& \sum_{d = 1}^n \mu(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \lfloor \dfrac{n}{id} \rfloor \sum_{j = 1}^{\lfloor \frac{m}{d} \rfloor} \lfloor \dfrac{m}{jd} \rfloor \\ \end{aligned} \]

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 5e4 + 7;

ll f[N];
int pri[N], mu[N], sum[N];
bool isp[N];

int T, n, m, pcnt;

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, mu[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, mu[i] = -1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                mu[i * pri[j]] = -mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = sum[i - 1] + mu[i];

    for (int i = 1; i < N; ++i)
        for (int l = 1, r; l <= i; l = r + 1) {
            r = i / (i / l);
            f[i] += 1ll * (i / l) * (r - l + 1);
        }
}

signed main() {
    sieve();
    scanf("%d", &T);

    while (T--) {
        scanf("%d%d", &n, &m);

        if (n > m)
            swap(n, m);

        ll ans = 0;

        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans += f[n / l] * f[m / l] * (sum[r] - sum[l - 1]);
        }

        printf("%lld\n", ans);
    }

    return 0;
}

P6055 [RC-02] GCD

求：

\[\sum_{i = 1}^{n} \sum_{j = 1}^n \sum_{p = 1}^{\lfloor \frac{n}{j} \rfloor} \sum_{q = 1}^{\lfloor \frac{n}{j} \rfloor} [\gcd(i, j) = 1] [\gcd(p, q) = 1] \pmod{998244353} \]

\(n \le 2 \times 10^9\)

\[\begin{align} & \sum_{i = 1}^{n} \sum_{j = 1}^n \sum_{p = 1}^{\lfloor \frac{n}{j} \rfloor} \sum_{q = 1}^{\lfloor \frac{n}{j} \rfloor} [\gcd(i, j) = 1] [\gcd(p, q) = 1] \\ =& \sum_{i = 1}^{n} \sum_{j = 1}^n [\gcd(i, j) = 1] \sum_{p = 1}^{n} \sum_{q = 1}^{n} [\gcd(p, q) = j] \\ =& \sum_{i = 1}^{n} \sum_{p = 1}^{n} \sum_{q = 1}^{n} [\gcd(i, p, q) = 1] \\ =& \sum_{d = 1}^n \mu(d) \lfloor \dfrac{n}{d} \rfloor^3 \end{align} \]

杜教筛处理 \(\mu\) 的前缀和即可。

#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353;
const int N = 1e7 + 7;

map<int, int> mp;

int pri[N], mu[N], sum[N];
bool isp[N];

int n, pcnt;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, mu[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, mu[i] = Mod - 1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                mu[i * pri[j]] = Mod - mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = add(sum[i - 1], mu[i]);
}

int Sum(int n) {
    if (n < N)
        return sum[n];

    if (mp.find(n) != mp.end())
        return mp[n];

    int res = 1;

    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res = dec(res, 1ll * (r - l + 1) * Sum(n / l) % Mod);
    }

    return mp[n] = res;
}

signed main() {
    sieve();
    scanf("%d", &n);
    int ans = 0;

    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans = add(ans, 1ll * (n / l) * (n / l) % Mod * (n / l) % Mod * dec(Sum(r), Sum(l - 1)) % Mod);
    }

    printf("%d", ans);
    return 0;
}

P3768 简单的数学题

求：

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

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

\[\begin{align} & \sum_{i = 1}^n \sum_{j = 1}^n ij \gcd(i, j) \\ =& \sum_{d = 1}^n \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} d \times id \times jd \times [\gcd(i, j) = 1] \\ =& \sum_{d = 1}^n d^3 \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} ij \sum_{k | i, k | j} \mu(k) \\ =& \sum_{d = 1}^n d^3 \sum_{k = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(k) \sum_{i = 1}^{\lfloor \frac{n}{dk} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{dk} \rfloor} ijk^2 \\ =& \sum_{d = 1}^n d^3 \sum_{k = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(k) k^2 S(\lfloor \dfrac{n}{dk} \rfloor)^2 \end{align} \]

其中 \(S(n) = \frac{n(n + 1)}{2}\) ，令 \(T = dk\) 则：

\[\begin{align} =& \sum_{T = 1}^n T^2 S(\lfloor \dfrac{n}{T} \rfloor)^2 \sum_{d | T} d \mu(\dfrac{T}{d}) \\ =& \sum_{T = 1}^n S(\lfloor \dfrac{n}{T} \rfloor)^2 T^2 \varphi(T) \end{align} \]

杜教筛处理 \(T^2 \varphi(T)\) 的前缀和即可，具体为构造 \(g(n) = n^2\) 。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;

map<ll, int> mp;

int pri[N], phi[N], sum[N];
bool isp[N];

ll n;
int Mod, pcnt, inv6;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, phi[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, phi[i] = i - 1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                phi[i * pri[j]] = phi[i] * phi[pri[j]];
            else {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = add(sum[i - 1], 1ll * i * i % Mod * phi[i] % Mod);
}

inline int Sum1(ll n) {
    n %= Mod;
    return n * (n + 1) / 2 % Mod;
}

inline int Sum2(ll n) {
    n %= Mod;
    return n * (n + 1) % Mod * (n * 2 + 1) % Mod * inv6 % Mod;
}

int Sum(ll n) {
    if (n < N)
        return sum[n];

    if (mp.find(n) != mp.end())
        return mp[n];

    int res = 1ll * Sum1(n) * Sum1(n) % Mod;

    for (ll l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res = dec(res, 1ll * dec(Sum2(r), Sum2(l - 1)) * Sum(n / l) % Mod);
    }

    return mp[n] = res;
}

signed main() {
    scanf("%d%lld", &Mod, &n);
    sieve();
    inv6 = mi(6, Mod - 2);
    int ans = 0;

    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans = add(ans, 1ll * Sum1(n / l) * Sum1(n / l) % Mod * dec(Sum(r), Sum(l - 1)) % Mod);
    }

    printf("%d", ans);
    return 0;
}

P3172 [CQOI2015] 选数

求：

\[\sum_{i_1 = L}^R \sum_{i_2 = L}^R \cdots \sum_{i_n = L}^R [\gcd(i_1, i_2, \cdots, i_n) = k] \pmod{10^9 + 7} \]

\(n, k \le 10^9\) ，\(L, R \le 10^9\) ，\(R - L \le 10^5\)

令 \(l = \lceil \dfrac{L}{k} \rceil, r = \lfloor \dfrac{R}{k} \rfloor\) ，则：

\[\begin{align} & \sum_{i_1 = L}^R \sum_{i_2 = L}^R \cdots \sum_{i_n = L}^R [\gcd(i_1, i_2, \cdots, i_n) = k] \\ =& \sum_{i_1 = l}^r \sum_{i_2 = l}^r \cdots \sum_{i_n = l}^r [\gcd(i_1, i_2, \cdots, i_n) = 1] \\ =& \sum_{i_1 = l}^r \sum_{i_2 = l}^r \cdots \sum_{i_n = l}^r \sum_{d | i_1, d | i_2, \cdots, d | i_n} \mu(d) \\ =& \sum_{d = 1}^r \mu(d) (\lfloor \dfrac{r}{d} \rfloor - \lfloor \dfrac{l - 1}{d} \rfloor)^n \end{align} \]

#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1e7 + 7;

map<int, int> mp;

int pri[N], mu[N], sum[N];
bool isp[N];

int n, k, L, R, pcnt;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, mu[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, mu[i] = Mod - 1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                mu[i * pri[j]] = Mod - mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = add(sum[i - 1], mu[i]);
}

inline int Sum(int n) {
    if (n < N)
        return sum[n];

    if (mp.find(n) != mp.end())
        return mp[n];

    int res = 1;

    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res = dec(res, 1ll * (r - l + 1) * Sum(n / l) % Mod);
    }

    return mp[n] = res;
}

signed main() {
    sieve();
    scanf("%d%d%d%d", &n, &k, &L, &R);
    L = (L - 1) / k, R /= k;
    int ans = 0;

    for (int l = 1, r; l <= R; l = r + 1) {
        r = R / (R / l);

        if (L / l)
            r = min(r, L / (L / l));

        ans = add(ans, 1ll * mi(R / l - L / l, n) * dec(Sum(r), Sum(l - 1)) % Mod);
    }

    printf("%d", ans);
    return 0;
}

P4619 [SDOI2018] 旧试题

给定 \(A, B, C\) ，求：

\[\sum_{i = 1}^A \sum_{j = 1}^B \sum_{k = 1}^C d(ijk) \]

\(A, B, C \le 2 \times 10^5\)

首先套路地推莫反：

\[\begin{aligned} & \sum_{i = 1}^A \sum_{j = 1}^B \sum_{k = 1}^C d(ijk) \\ =& \sum_{i = 1}^A \sum_{j = 1}^B \sum_{k = 1}^C \sum_{a \mid i} \sum_{b \mid k} \sum_{c \mid k} [\gcd(a, b) = 1] [\gcd(a, c) = 1] [\gcd(b, c) = 1] \\ =& \sum_{a = 1}^A \sum_{b = 1}^B \sum_{c = 1}^C [\gcd(a, b) = 1] [\gcd(a, c) = 1] [\gcd(b, c) = 1] \lfloor \frac{A}{a} \rfloor \lfloor \frac{B}{b} \rfloor \lfloor \frac{C}{c} \rfloor \\ =& \sum_{x = 1}^{\min(A, B)} \sum_{y = 1}^{\min(A, C)} \sum_{z = 1}^{\min(B, C)} \mu(x) \mu(y) \mu(z) \sum_{\mathrm{lcm}(x, y) \mid a} \lfloor \frac{A}{a} \rfloor \sum_{\mathrm{lcm}(x, z) \mid b} \lfloor \frac{B}{b} \rfloor \sum_{\mathrm{lcm}(y, z) \mid c} \lfloor \frac{C}{c} \rfloor \\ =& \sum_{x = 1}^{\min(A, B)} \sum_{y = 1}^{\min(A, C)} \sum_{z = 1}^{\min(B, C)} \mu(x) \mu(y) \mu(z) \sum_{a = 1}^{\lfloor \frac{A}{\mathrm{lcm}(x, y)} \rfloor} \lfloor \frac{\lfloor \frac{A}{\mathrm{lcm}(x, y)} \rfloor}{a} \rfloor \sum_{b = 1}^{\lfloor \frac{B}{\mathrm{lcm}(x, z)} \rfloor} \lfloor \frac{\lfloor \frac{B}{\mathrm{lcm}(x, z)} \rfloor}{b} \rfloor \sum_{c = 1}^{\lfloor \frac{C}{\mathrm{lcm}(y, z)} \rfloor} \lfloor \frac{\lfloor \frac{C}{\mathrm{lcm}(y, z)} \rfloor}{c} \rfloor \end{aligned} \]

记 \(f(n) = \sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor\) ，这可以 \(O(V \sqrt{V})\) 数论分块预处理，继续推式子得到：

\[\sum_{x = 1}^{\min(A, B)} \sum_{y = 1}^{\min(A, C)} \sum_{z = 1}^{\min(B, C)} \mu(x) \mu(y) \mu(z) F \left( \lfloor \frac{A}{\mathrm{lcm}(x, y)} \rfloor \right) F \left( \lfloor \frac{B}{\mathrm{lcm}(x, z)} \rfloor \right) F \left( \lfloor \frac{C}{\mathrm{lcm}(y, z)} \rfloor \right) \]

事实上很多 \((x, y, z)\) 的贡献都是 \(0\) ，有贡献的三元组一定满足：

  • \(\mu(x) \ne 0\) ，\(\mu(y) \ne 0\) ，\(\mu(z) \ne 0\) 。
  • \(\max \{ \mathrm{lcm}(x, y), \mathrm{lcm}(x, z), \mathrm{lcm}(y, z) \} \le \max \{ A, B, C \}\) 。

考虑将可能有贡献的点之间连边，枚举 \(\gcd = i\) ，然后拿出 \(i\) 的倍数 \(j, k\) 满足 \(\mu(j) \ne 0\) 且 \(\mu(k) \ne 0\) 连边。问题转化为三元环计数，可能需要额外处理一下有数字相等的情况。

事实上这张图的边数是很少的 \(m = 821535\) ，因此直接跑三元环计数即可做到 \(O(m \sqrt{m})\) 。

#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1e5 + 7, S = 1e6 + 7;

vector<int> g[N];

struct Edge {
    int u, v, w;
} e[S];

int pri[N], mu[N], deg[N], f[N], tag[N];
bool isp[N];

int A, B, C, n, tot, pcnt;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, mu[1] = 1;
    
    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, mu[i] = -1;
        
        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;
            
            if (i % pri[j])
                mu[i * pri[j]] = -mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }
}

inline int solve() {
    for (int i = 1; i <= n; ++i)
        g[i].clear();
    
    for (int i = 1; i <= tot; ++i) {
        int &u = e[i].u, &v = e[i].v;
        
        if (!(deg[u] == deg[v] ? u < v : deg[u] < deg[v]))
            swap(u, v);
        
        g[u].emplace_back(i);
    }
    
    int ans = 0;
    
    for (int i = 1; i <= n; ++i) {
        for (int j : g[i])
            tag[e[j].v] = e[j].w;
        
        for (int j : g[i])
            for (int k : g[e[j].v]) {
                if (!tag[e[k].v])
                    continue;

                int u = i, v = e[j].v, w = e[k].v, lcma = tag[e[k].v], lcmb = e[j].w, lcmc = e[k].w, res = 0;

                auto calc = [&](int a, int b, int c) {
                    return 1ll * f[a / lcma] * f[b / lcmb] % Mod * f[c / lcmc] % Mod;
                };
                
                if (u == v && v == w)
                    res = add(res, calc(A, B, C));
                else if (u == v || u == w || v == w)
                    res = add(res, add(calc(A, B, C), add(calc(B, C, A), calc(C, A, B))));
                else {
                    res = add(res, add(calc(A, B, C), calc(A, C, B)));
                    res = add(res, add(calc(B, A, C), calc(B, C, A)));
                    res = add(res, add(calc(C, A, B), calc(C, B, A)));
                }
                
                if (mu[u] * mu[v] * mu[w] == 1)
                    ans = add(ans, res);
                else
                    ans = dec(ans, res);
            }
        
        for (int j : g[i])
            tag[e[j].v] = 0;
    }
    
    return ans;
}

signed main() {
    sieve();
    
    for (int i = 1; i < N; ++i)
        for (int l = 1, r; l <= i; l = r + 1)
            r = i / (i / l), f[i] = add(f[i], 1ll * (r - l + 1) * (i / l) % Mod);
    
    int T;
    scanf("%d", &T);
    
    while (T--) {
        scanf("%d%d%d", &A, &B, &C);
        memset(deg + 1, 0, sizeof(int) * (n = max(max(A, B), C))), tot = 0;
        
        for (int i = 1; i <= n; ++i)
            for (int j = i; j <= n; j += i) {
                if (!mu[j])
                    continue;

                for (int k = j; 1ll * j / i * k <= n; k += i)
                    if (mu[k] && __gcd(j, k) == i)
                        e[++tot] = (Edge){j, k, j / i * k}, ++deg[j], ++deg[k];
            }
        
        printf("%d\n", solve());
    }
    
    return 0;
}

单位根反演

公式：

\[[n \mid k] = \frac{1}{n} \sum_{i = 0}^{n - 1} \omega_n^{ki} \]

证明可以等比数列求和，也可以用向量加法。

推论：

\[[x \equiv y \pmod{n}] = \frac{1}{n} \sum_{i = 0}^{n - 1} \omega_n^{(x - y)i} \]

LOJ6485. LJJ 学二项式定理

给定 \(n, s, a_{0 \sim 3}\) ，求：

\[\left( \sum_{i = 0}^n \binom{n}{i} \times s^i \times a_{i \bmod 4} \right) \pmod{998244353} \]

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

\[\begin{aligned} & \sum_{i = 0}^n \binom{n}{i} \times s^i \times a_{i \bmod 4} \\ =& \sum_{i = 0}^n \binom{n}{i} \times s^i \times \sum_{j = 0}^3 a_j [i \bmod 4 = j] \\ =& \sum_{i = 0}^n \binom{n}{i} \times s^i \times \sum_{j = 0}^3 a_j \times \frac{1}{4}\sum_{k = 0}^3 \omega_4^{(i - j) k} \\ =& \frac{1}{4} \sum_{j = 0}^3 a_j \sum_{k = 0}^3 \omega_4^{-jk} \sum_{i = 0}^n \binom{n}{i} (s \omega_4^k)^i \\ =& \frac{1}{4} \sum_{j = 0}^3 a_j \sum_{k = 0}^3 \omega_4^{-jk} (s \omega_4^k + 1)^n \\ \end{aligned} \]

时间复杂度 \(O(\log n)\) 。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 998244353, buf = 911660635, inv4 = 1ll * ((Mod + 1) / 2) * ((Mod + 2) / 2) % Mod;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

signed main() {
    int T;
    scanf("%d", &T);

    while (T--) {
        ll n;
        int s, a[4], ans = 0;
        scanf("%lld%d%d%d%d%d", &n, &s, a + 0, a + 1, a + 2, a + 3);

        for (int i = 0; i <= 3; ++i)
            for (int j = 0, pw = 1; j <= 3; ++j, pw = 1ll * pw * buf % Mod)
                ans = add(ans, 1ll * a[i] * mi(buf, Mod - 1 - i * j) % Mod * 
                    mi(add(1ll * s * pw % Mod, 1), n % (Mod - 1)) % Mod);

        printf("%d\n", 1ll * inv4 * ans % Mod);
    }

    return 0;
}

P5591 小猪佩奇学数学

给定 \(n, p, k\) ，求：

\[\sum_{i = 0}^n \binom{n}{i} \times p^i \times \lfloor \frac{i}{K} \rfloor \pmod{998244353} \]

\(1 \le n, p < 998244353\) ，\(k \in \{ 2^w \mid 0 \le w \le 20 \}\)

用 \(\lfloor \frac{i}{k} \rfloor = \frac{i - i \bmod k}{k}\) 展开得到：

\[\begin{aligned} & \sum_{i = 0}^n \binom{n}{i} \times p^i \times \lfloor \frac{i}{k} \rfloor \\ =& \frac{1}{k} \sum_{i = 0}^n \binom{n}{i} \times p^i \times (i - i \bmod k) \\ =& \frac{1}{k} \left( \sum_{i = 0}^n \binom{n}{i} p^i i - \sum_{i = 0}^n \binom{n}{i} p^i (i \bmod k) \right) \end{aligned} \]

由吸收恒等式得到：

\[\sum_{i = 0}^n \binom{n}{i} p^i i = n p \sum_{i = 0}^n \binom{n - 1}{i - 1} p^{i - 1} = n p (p + 1)^{n - 1} \]

继续考虑后面一项：

\[\begin{aligned} &\sum_{i = 0}^n \binom{n}{i} p^i (i \bmod k) \\ =& \frac{1}{k} \sum_{i =0 }^n \binom{n}{i} p^i \sum_{j = 0}^{k - 1} j [i \bmod k = j] \\ =& \frac{1}{k} \sum_{i =0 }^n \binom{n}{i} p^i \sum_{j = 0}^{k - 1} j \sum_{t = 0}^{k - 1} \omega_k^{t(i - j)} \\ =& \frac{1}{k} \sum_{j = 0}^{k - 1} j \sum_{t = 0}^{k - 1} \omega_k^{-t j} \sum_{i = 0}^n \binom{n}{i} p^i \omega_k^{t i} \\ =& \frac{1}{k} \sum_{j = 0}^{k - 1} j \sum_{t = 0}^{k - 1} \omega_k^{-t j} (p \omega_k^t + 1)^n \\ =& \frac{1}{k} \sum_{t = 0}^{k - 1} (p \omega_k^t + 1)^n \sum_{j = 0}^{k - 1} j \omega_k^{-t j} \end{aligned} \]

考虑后面的 \(\sum_{j = 0}^{k - 1} j \omega_n^{-t j}\) ，设 \(S(n, k) = \sum_{i = 0}^{n - 1} i k^i\) ，则答案为：

\[\frac{1}{k} (n p (p + 1)^{n - 1} - \frac{1}{k} \sum_{t = 0}^{k - 1} (p \omega_k^t + 1)^n S(k - 1, \omega_k^{-t})) \]

下面考虑算 \(S\) 。当 \(k = 1\) 时有 \(S(n, 1) = \sum_{i = 0}^n i = \frac{(n - 1)n}{2}\) ，当 \(k \not = 1\) 时有：

\[\begin{aligned} k S(n, k) - S(n, k) &= \sum_{i = 0}^{n - 1} i k^{i + 1} - \sum_{i = 0}^{n - 1} i k^i \\ &= \sum_{i = 1}^n (i - 1) k^i - \sum_{i = 0}^{n - 1} i k^i \\ &= (n - 1) k^n - \sum_{i = 1}^{n - 1} k^i \\ &= (n - 1) k^n + \frac{k - k^{n + 1}}{k - 1} + 1 \\ &= (n - 1) k^n + \frac{1 - k^n}{k - 1} + 1 \end{aligned} \]

可得 \(S(n, k) = \frac{(n - 1) k^n + \frac{1 - k^n}{k - 1} + 1}{k - 1}\) 。注意到带入的 \(k\) 总是 \(n\) 次单位根，故 \(k^n = 1\) ，则 \(S(n, k) = \frac{n}{k - 1}\) 。

时间复杂度 \(O(k \log n)\) 。

#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353, rt = 3;

int n, p, k;

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline int S(int n, int k) {
    return k == 1 ? (1ll * n * (n - 1) / 2) % Mod : 1ll * n * mi(k - 1, Mod - 2) % Mod;
}

signed main() {
    scanf("%d%d%d", &n, &p, &k);
    vector<int> pw(k + 1);
    pw[0] = 1, pw[1] = mi(rt, (Mod - 1) / k);
    
    for (int i = 2; i <= k; ++i)
        pw[i] = 1ll * pw[i - 1] * pw[1] % Mod;

    int ans = 0;
    
    for (int i = 0; i < k; ++i)
        ans = add(ans, 1ll * mi(add(1ll * p * pw[i] % Mod, 1), n) * S(k, pw[k - i]) % Mod);
    
    printf("%d", 1ll * mi(k, Mod - 2) * 
        dec(1ll * n * p % Mod * mi(p + 1, n - 1) % Mod, 1ll * ans * mi(k, Mod - 2) % Mod) % Mod);
    return 0;
}

杜教筛

杜教筛被用于处理一类数论函数的前缀和问题，可以在低于线性时间复杂度内计算某些数论函数的前缀和。

考虑构造一个 \(S(n)\) 关于 \(S(\lfloor \frac{n}{i} \rfloor)\) 的递推式。对于任意数论函数 \(g\) ，必有：

\[\sum_{i = 1}^n (f * g)(i) = \sum_{i = 1}^n \sum_{d | i} g(d) f(\frac{i}{d}) = \sum_{i = 1}^n g(i) S(\lfloor \frac{n}{i} \rfloor) \]

那么可以得到递推式：

\[\begin{align} g(1) S(n) &= \sum_{i = 1}^n g(i) S(\lfloor \frac{n}{i} \rfloor) - \sum_{i = 2}^n g(i) S(\lfloor \dfrac{n}{i} \rfloor) \\ &= \sum_{i = 1}^n (f * g)(i) - \sum_{i = 2}^n g(i) S(\lfloor \frac{n}{i} \rfloor) \end{align} \]

假如可以构造恰当的数论函数 \(g\) 使得：

  • 可以快速计算 \(\sum_{i = 1}^n (f * g)(i)\) 。
  • 可以快速计算 \(g\) 的单点值，以便于用整除分块求解 \(\sum_{i = 2}^n g(i) S(\lfloor \frac{n}{i} \rfloor)\) 。

则可以用杜教筛求出 \(S(n)\) 。预处理 \(S(1 \sim n^{\frac{2}{3}})\) 的值即可做到 \(O(n^{\frac{2}{3}})\) 的复杂度，注意需要记忆化保证复杂度正确性。

P4213【模板】杜教筛（Sum）

求欧拉函数、莫比乌斯函数的前 \(n\) 项和。

\(n < 2^{31}\)

对于欧拉函数前缀和的求解，由 \(\varphi * 1 = id\) ，从而：

\[S(n) = \dfrac{n (n + 1)}{2} - \sum_{i = 2}^n S(\lfloor \frac{n}{i} \rfloor) \]

对于莫比乌斯函数前缀和的求解，由 \(\mu * 1 = \epsilon\) ，从而：

\[S(n) = 1 - \sum_{i = 2}^n S(\lfloor \frac{n}{i} \rfloor) \]

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;

int pri[N];
bool isp[N];

int pcnt;

namespace Phi {
map<ll, ll> mp;

ll s[N];
int f[N];

inline void prework() {
    for (int i = 1; i < N; ++i)
        s[i] = s[i - 1] + f[i];
}

inline ll Sum(ll n) {
    if (n < N)
        return s[n];

    if (mp.find(n) != mp.end())
        return mp[n];

    ll res = 1ll * n * (n + 1) / 2;

    for (ll l = 2, r; l <= n; l = r + 1)
        r = n / (n / l), res -= 1ll * (r - l + 1) * Sum(n / l);

    return mp[n] = res;
}
} // namespace Phi

namespace Mu {
map<ll, ll> mp;

ll s[N];
int f[N];

inline void prework() {
    for (int i = 1; i < N; ++i)
        s[i] = s[i - 1] + f[i];
}

inline ll Sum(ll n) {
    if (n < N)
        return s[n];

    if (mp.find(n) != mp.end())
        return mp[n];

    ll res = 1;

    for (ll l = 2, r; l <= n; l = r + 1)
        r = n / (n / l), res -= 1ll * (r - l + 1) * Sum(n / l);

    return mp[n] = res;
}
} // namespace Mu

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, Phi::f[1] = Mu::f[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, Phi::f[i] = i - 1, Mu::f[i] = -1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                Phi::f[i * pri[j]] = Phi::f[i] * (pri[j] - 1), Mu::f[i * pri[j]] = -Mu::f[i];
            else {
                Phi::f[i * pri[j]] = Phi::f[i] * pri[j], Mu::f[i * pri[j]] = 0;
                break;
            }
        }
    }

    Phi::prework(), Mu::prework();
}

signed main() {
    sieve();
    int T;
    scanf("%d", &T);

    while (T--) {
        ll n;
        scanf("%lld", &n);
        printf("%lld %lld\n", Phi::Sum(n), Mu::Sum(n));
    }

    return 0;
}

P3768 简单的数学题

给定 \(n, p\) ，求：

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

\(n \le 10^{10}\) ，\(5 \times 10^8 \le p \le 1.1 \times 10^9\) 且 \(p\) 为质数

考虑将 \(\gcd\) 用 \(\varphi\) 反演掉：

\[\begin{align} & \sum_{i = 1}^n \sum_{j = 1}^n ij \gcd(i, j) \\ =& \sum_{i = 1}^n \sum_{j = 1}^n ij \sum_{d | i, d | j} \varphi(d) \\ =& \sum_{d = 1}^n \varphi(d) \sum_{d | i} \sum_{d | j} ij \\ =& \sum_{d = 1}^n \varphi(d) d^2 (\sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} i)^2 \\ =& \sum_{d = 1}^n \varphi(d) d^2 \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} i^3 \\ \end{align} \]

直接杜教筛 \(\varphi(d) d^2\) ，为了消掉 \(d^2\) 取 \(g = id^2\) 即可。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;

map<ll, ll> mp;

ll sum[N];
int pri[N], phi[N];
bool isp[N];

ll n, Mod, inv6;
int pcnt;

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline void sieve() {
    memset(isp, true, sizeof(isp));
    isp[1] = false, phi[1] = 1;

    for (int i = 2; i < N; ++i) {
        if (isp[i])
            pri[++pcnt] = i, phi[i] = i - 1;

        for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
            isp[i * pri[j]] = false;

            if (i % pri[j])
                phi[i * pri[j]] = phi[i] * (pri[j] - 1);
            else {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
        }
    }

    for (int i = 1; i < N; ++i)
        sum[i] = (sum[i - 1] + 1ll * phi[i] * i % Mod * i % Mod) % Mod;
}

inline ll S1(ll n) {
    n %= Mod;
    return n * (n + 1) / 2 % Mod;
}

inline ll S2(ll n) {
    n %= Mod;
    return n * (n + 1) % Mod * (n * 2 + 1) % Mod * inv6 % Mod;
}

inline ll S3(ll n) {
    return S1(n) * S1(n) % Mod;
}

ll Sum(ll n) {
    if (n < N)
        return sum[n];

    if (mp.find(n) != mp.end())
        return mp[n];

    ll res = S3(n);

    for (ll l = 2, r; l <= n; l = r + 1)
        r = n / (n / l), res = (res - 1ll * (S2(r) - S2(l - 1) + Mod) % Mod * Sum(n / l) % Mod + Mod) % Mod;

    return mp[n] = res;
}

signed main() {
    scanf("%lld%lld", &Mod, &n);
    sieve(), inv6 = mi(6, Mod - 2);
    ll ans = 0;

    for (ll l = 1, r; l <= n; l = r + 1)
        r = n / (n / l), ans = (ans + 1ll * (Sum(r) - Sum(l - 1) + Mod) % Mod * S3(n / l) % Mod) % Mod;

    printf("%lld", ans);
    return 0;
}
posted @ 2025-06-26 13:43  wshcl  阅读(44)  评论(0)    收藏  举报