数论进阶 - james1 的博客

高级的数论知识更为困难，但也更有意思。本文将介绍省选以内的常见数论知识点，进一步学习模意义下的数论、数论函数和整除性相关问题，并通过题目介绍它们的应用。

整除性问题

与整除相关的内容是数论的根基，我们再来梳理一下这些内容。

算术基本定理

之前介绍过算术基本定理的标准分解式：

$n=\prod p_i^{\alpha_i}$

它还可以写成 $p$ 进赋值序列。记 $v_p(n)=\max\{k\in \mathbb{N}\mid p^k\mid n\}$ ，定义 $n$ 的 $p$ 进赋值序列为： $\{v_2(n),v_3(n),v_5(n),v_7(n),\cdots\}$ 。

这有什么用呢？ $p$ 进赋值序列是从数论角度刻画正整数的重要工具，在这个视角下，每个整数均可以定位在一个无穷维空间上的整点，每个维度都对应一个素数。

取模与整除

首先： $a\bmod b=a-b\times \left\lfloor\frac{a}{b}\right\rfloor$ ，实现了取模与整除的转化，需要根据具体情况分析使用哪种形式。

而 $\left\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\right\rfloor=\left\lfloor\frac{a}{bc}\right\rfloor$ 可以将一部分的除法问题转化为乘法问题。

数论分块

数论分块。 $\left\lfloor\frac{n}{i}\right\rfloor$ 只有 $O(\sqrt{n})$ 中不同的取值，并且每一种取值都是一个连续的区间。

比如光速幂是利用这一点来实现的，这个技巧也可以在 $O(\sqrt{n})$ 时间内枚举到所有区间。

最经典的应用就是快速计算 $\sum_{i=1}^n f(i)g\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$ ，需要 $O(1)$ 计算 $\sum_{i=l}^{r}f(i)$ ，然后将 $g$ 相同的打包计算。

使得 $\lfloor\frac{n}{i}\rfloor=\lfloor\frac{n}{j}\rfloor$ 成立的最大满足 $i\le j\le n$ 的块的右端点为 $\left\lfloor\cfrac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor$ 。为什么？令 $k=\lfloor\frac{n}{i}\rfloor$ ，可知 $k\le \frac{n}{i}$ ，那么 $\lfloor\frac{n}{k}\rfloor\ge i$ ，所以 $j=\frac{n}{i}$ 。时间复杂度 $O(\sqrt{n})$ 。

模板，代码如下：

查看代码

#include <iostream>
#include <cstdio>

using namespace std;
typedef long long i64;

i64 H(i64 n) {
    i64 res = 0; i64 l = 1, r;
    while (l <= n) {
        r = n / (n / l);
        res += (r - l + 1) * (n / l);
        l = r + 1;
    }
    return res;
}

int main(void) {
    int T; i64 x; cin >> T;
    while (T--) { cin >> x; cout << H(x) << "\n"; }
    return 0;
}

同样可以扩展到有多个形如 $\lfloor\frac{n}{i}\rfloor$ 的求和，只需要对 $r$ 取 $\min$ ，时间复杂度为 $O(\sum \sqrt{n})$ 。

扩展欧几里得算法

辗转相除法和更相减损法是一个很重要的结构，出现时一般都伴随着 $\gcd$ 相关的结论。

类欧几里得算法

大致意义是使用一个类似辗转相除法的方法进行函数求和。

设：

$f(a,b,c,n)=\sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor$

需要给出一个 $O(\log n)$ 的算法。

数论分块似乎应对不了这样的求和，因此我们考虑一个奇怪的方式。要将 $a\ge c$ 或 $b\ge$ 的情况转化为 $a,b<c$ 。

我们有：

$\begin{aligned} f(a,b,c,n) &=\sum_{i=0}^n \left\lfloor\frac{(a\bmod c)i+(a-a\bmod c)i+(b\bmod c)+(b-b\bmod c)}{c}\right\rfloor\\ &=\sum_{i=0}^n \left\lfloor\frac{(a\bmod c)i+(b\bmod c)}{c}\right\rfloor +\frac{(a-a\bmod c)i}{c}+\frac{b-b\bmod c}{c} \end{aligned}$

注意第二行中的后两个一定是可以整除的。

因此原式等于：

$\begin{aligned} &=\sum_{i=0}^n \left\lfloor\frac{(a\bmod c)i+(b\bmod c)}{c}\right\rfloor+\left\lfloor\frac{a}{c}\right\rfloor x+\left\lfloor\frac{b}{c}\right\rfloor \\ &=\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor+(n+1)\left\lfloor\frac{b}{c}\right\rfloor+ \sum_{i=0}^n\left\lfloor\frac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c} \right\rfloor\\ &=\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor +(n+1)\left\lfloor\frac{b}{c}\right\rfloor+f(a\bmod c,b\bmod c,c,n) \end{aligned}$

现在有 $a,b<c$ 了，然后呢？式子中只有 $i$ 一个变量。原式中的条件是 $0\le i\le n$ ，贡献是 $\left\lfloor\cfrac{ai+b}{c}\right\rfloor$ 。考虑将贡献合并计算，将式子改为：

$\sum_{i=0}^n\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1}1\\$

这个 $j$ 被 $i$ 限制着，不好做，我们交换 $i,j$ ，让 $j$ 被 $n$ 限制：

$\sum_{j=0}^{\left\lfloor \frac{an+b}{c} \right\rfloor-1}\sum_{i=0}^n\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor\right]$

万能欧几里得算法

万能欧几里得可以解决几乎全部类欧几里得算法能解决的问题。其优势在于，不管要求什么式子的值，代码都类似。

数论函数

我们一般研究数论函数中的积性函数。

常见积性函数

单位函数 $\epsilon(n)=[n=1]$ ，是完全积性函数。
常数函数 $1(n)=1$ ，是完全积性函数。
恒等函数 $\operatorname{id}_k(n)=n^k$ ，是完全积性函数，当 $k=1$ 时可以省略不写。
除数函数 $\sigma_k(n)=\sum_{d\mid n}d^k$ 。这样的话， $\sigma_0(n)=\tau(n)$ ， $\sigma_1(n)$ 代表约数和，有 $\sigma_k(n)=\prod_{i=1}^s\left(\sum_{j=0}^{\alpha_i}p_i^{jk}\right)=\prod_{i=1}^s\sigma_k(p_i^{\alpha_i})$ 。
欧拉函数， $\varphi(n)$ 代表 $1\sim n$ 中与 $n$ 互质的数的个数，你肯定不会忘记它。
本质不同质因子个数函数 $\omega(n)=\sum_{p}[p\mid n]$ 。

线性筛求积性函数

方法很简单，记录 $low$ 表示最小质因子的最高次幂，然后分类讨论一下即可。

for (int i = 2; i <= n; i++) {
    if (!vis[i]) prime[++tot] = i, f[i] = ..., low[i] = i;  // 单独算 f(p)
    for (int j = 1; j <= tot && i * prime[j] <= n; j++) {
        vis[i * prime[j]] = 1;
        if (i % prime[j] == 0) {  // i 与 p 不互质
            low[i * prime[j]] = low[i] * prime[j];
            if (i == low[i]) f[i * prime[j]] = ...;  // i = p^k，单独算 f(p^{k+1})
            else f[i * prime[j]] = f[i / low[i]] * f[low[i] * prime[j]];
            break;
        }
        low[i * prime[j]] = prime[j];
        f[i * prime[j]] = f[i] * f[prime[j]];  // i 与 p 互质，f(ip) = f(i)f(p)
    }
}

狄利克雷卷积

狄利克雷卷积是数论函数的基本运算。定义：

$h(n)=\sum_{d\mid n}f(d)g\left(\frac{n}{d}\right)$

简记为 $h=f*g$ 。

狄利克雷卷积具有交换律、结合律和分配律。

首先， $\epsilon*f=f$ 。因此单位函数 $\epsilon$ 为狄利克雷卷积的单位元，那么就可以定义数论函数的逆元 $f^{-1}$ ，满足 $f*f^{-1}=\epsilon$ 。

一个 $f$ 存在逆元，当且仅当 $f(1)\ne 0$ ，并且逆元唯一。 $f=g$ 的充要条件是 $f*h=g*h(h(1)\ne 1)$ 。这样一来有 $g(n) = -\cfrac{\sum\limits_{d \mid n, d \neq n} g(d)f\left(\dfrac n d\right)} {f(1)}$ 。

积性函数的狄利克雷卷积是积性函数，积性函数的逆元也是积性函数。

以下是一些常用的狄利克雷卷积：

$1*1=\tau$ ，
$\varphi *1=\operatorname{id}$ 。

狄利克雷前缀和

任意数论函数 $f$ 卷常数函数 $1$ 等价于做 $f$ 的狄利克雷前缀和，即： $g=f*1,g(n)=\sum_{d\mid n}f(d)$ 。含义是对每个 $n$ 计算给定数论函数在其所因数处的取值和。

将每个数写成无穷序列 $a_n=\{c_1,c_2,\cdots\}$ 表示 $n=\prod p_i^{c_i}$ 。由于 $x\mid y$ 的充要条件为 $a_x(c_i)\le a_y(c_i)$ ，因此 $f*1$ 可以看成对下标做其无穷序列的高维前缀和。

模板，采用高维前缀和实现即可。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef unsigned int uint; 

uint n, seed; 
inline uint getnext(void) {
	seed ^= seed << 13, seed ^= seed >> 17, seed ^= seed << 5;
	return seed;
}

bool v[20000005];
uint prime[10000005], tot; 
uint a[20000005];

int main(void) {
    cin >> n >> seed;
    for (int i = 1; i <= n; ++i) a[i] = getnext(); 
    for (int i = 2; i <= n; ++i) {
        if (!v[i]) prime[++tot] = i; 
        for (int j = 1; j <= tot && i * prime[j] <= n; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break;
        }
    }
    for (int i = 1; i <= tot; ++i)
        for (int j = 1; j * prime[i] <= n; ++j)
            a[j * prime[i]] += a[j]; 
    uint ans = 0; 
    for (int i = 1; i <= n; ++i) ans ^= a[i];
    cout << ans << "\n"; return 0;
}

莫比乌斯反演

我们定义莫比乌斯函数：

$\mu(n)=\begin{cases}1,&n=1,\\0,&\exists d>1,d^{2}\mid n, \\(-1)^{\omega(n)},&\mathrm{otherwise}. \end{cases}$

使用线性筛求解莫比乌斯函数

mu[1] = 1; 
for (int i = 2; i <= n; ++i) {
    if (!v[i]) prime[++tot] = i, mu[i] = -1; 
    for (int j = 1; j <= tot && i * prime[j] <= n; ++j) {
        v[i * prime[j]] = 1; 
        if (i % prime[j] == 0) {
            mu[i * prime[j]] = 0; 
            break; 
        }
        mu[i * prime[j]] = -mu[i];
    }
}

我们定义莫比乌斯函数：

$\mu(n)=\begin{cases}1,&n=1,\\0,&\exists d>1,d^{2}\mid n, \\(-1)^{\omega(n)},&\mathrm{otherwise}. \end{cases}$

为什么要定义这样一个奇怪的函数呢？实际上它是在对 $\mathbb{N}$ 做容斥。设 $g(n)=\sum_{n\mid d} f(d)$ ，已知 $g$ ，要求 $f(1)$ 。 $f(1)$ 等于 $g$ 在 $1$ 的倍数处取的取值和，减去质数处的取值和，但是多减了两个质数乘积的地方，因此还要加回来。这就是容斥原理，容斥系数是 $(-1)^{\omega(n)}$ 。比如求无平方因子的数的个数就是像这样容斥。

莫比乌斯函数不仅是积性函数，还满足：

$\sum_{d\mid n}\mu(d)=\begin{cases} 1,&n=1,\\ 0,&n\neq 1. \end{cases}$

上述式子描述的其实是 $\mu*1=\epsilon$ ，这是莫比乌斯函数最重要的性质，其引出了性质： $\sum_{d\mid n}\mu(d)=[n=1]$ ，我们可以将这个判断式转化为和式，这样更加方便计算。两者的转化称之为莫比乌斯反演。

莫比乌斯反演有以下结论：

若 $g=f*1$ ，则 $f=g*\mu$ 。
若 $g(n)=\sum_{n\mid d}f(d)$ ，则 $f(n)=\sum_{n\mid d} \mu\left(\dfrac d n\right) g(d)$ 。这也被称为莫比乌斯变换。
由于 $\varphi * 1 = \operatorname{id}$ ，因此 $\operatorname{id}*\mu =\varphi$ 。

一个常见的应用是， $[\gcd(i,j)=1]=\sum_{d\mid \gcd(i,j)}\mu(d)$ 。虽然看起来这像是废话，但这一步将 $i,j$ 互质转化为了枚举 $\gcd(i,j)$ 的约数 $d$ 。如果 $i,j$ 同样需要枚举，那么枚举 $d$ 并计算合法的 $i,j$ 个数（当且仅当 $d\mid i,d\mid j$ ）即可。具体来说：

$\begin{aligned} \sum\limits_{i = 1} ^ n \sum\limits_{j = 1} ^ m [\gcd(i, j) = 1] & = \sum\limits_{i = 1} ^ n \sum\limits_{j = 1} ^ m \sum\limits_{d\mid \gcd(i, j)} \mu(d) \\ & = \sum\limits_{d = 1} ^ {\min(n, m)} \mu(d) \sum\limits_{i = 1} ^ n \sum\limits_{j = 1} ^ m [d\mid i\land d\mid j] \\ & = \sum\limits_{d = 1} ^ {\min(n, m)} \mu(d) \left\lfloor \dfrac n d \right\rfloor \left\lfloor \dfrac m d \right\rfloor \\ \end{aligned}$

就可以直接使用数论分块在 $O(\sqrt{n}+\sqrt{m})$ 的时间内计算。

关于莫比乌斯反演的更多应用请在题目中寻找。

欧拉反演

欧拉函数有一条性质： $n = \sum_{d\mid n} \varphi(d)$ 。

用与莫比乌斯反演一样的套路，有：

$\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid \gcd(i,j)} \varphi(d)\\ &=\sum_{d=1}^{\min(n,m)}\varphi(d)\sum\limits_{i = 1} ^ n \sum\limits_{j = 1} ^ m [d\mid i\land d\mid j] \\ &=\sum\limits_{d = 1} ^ {\min(n, m)} \varphi(d) \left\lfloor \frac n d \right\rfloor \left\lfloor \frac m d \right\rfloor \end{aligned}$

这就是欧拉反演的公式。

模意义下的数论

数论算法的综合应用

Problemset

这里的题可能有一些难，也会与别的算法综合，但是都应该完成。

简单问题

不适合放在《数论初步》，但是又没用到本文介绍的知识点的简单题。

[SP5971] LCMSUM - LCM Sum

Portal.

$\begin{aligned} \sum\limits_{i = 1} ^ n \operatorname{lcm}(i, n) & = n \sum\limits_{i = 1} ^ n \frac{i}{\gcd(i, n)} \\ & = n \sum\limits_{d\mid n} \sum\limits_{i = 1} ^ n \frac{i}{d} [\gcd(i, n) = d] \\ & = n \sum\limits_{d\mid n} \sum\limits_{i = 1} ^ {\frac n d} i \left[\gcd\left(i, \frac n d\right) = 1\right] \end{aligned}$

后面这个东西是在 $d$ 以内所有与 $d$ 互质的数的和，是欧拉函数的结论，答案是 $\frac{n\times \varphi(n)}{2}$ ，注意 $d=1$ 要特判掉。

线性筛出 $1*\operatorname{id}*\ \varphi$ 即可。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 
const int N = 1000000; 

int tot, prime[500005], low[1000005]; 
i64 f[1000005];

int main(void) {
    f[1] = 1; 
    for (int i = 2; i <= N; ++i) {
        if (!f[i]) prime[++tot] = low[i] = i, f[i] = 1 + 1ll * i * (i - 1); 
        for (int j = 1; j <= tot && i * prime[j] <= N; ++j) {
            if (i % prime[j] == 0) {
                low[i * prime[j]] = low[i] * prime[j];
                if (i == low[i]) f[i * prime[j]] = f[i] + 1ll * (i * prime[j]) * i * (prime[j] - 1); 
                else f[i * prime[j]] = f[i / low[i]] * f[low[i] * prime[j]];
                break;
            }
            low[i * prime[j]] = prime[j];
            f[i * prime[j]] = f[i] * f[prime[j]];
        } 
    }
    int T, n; scanf("%d", &T); 
    while (T--) {
        scanf("%d", &n); 
        printf("%lld\n", (f[n] + 1) * n / 2);
    }
    return 0;
}

[六省联考 2017] 相逢是问候

Portal.

给定 $p,c$ ，维护一个正整数序列，支持区间替换为 $c^{a_i}$ ，区间求和（对 $p$ 取模）。

数论分块

这里的数论分块不会很难，因为其往往会和莫比乌斯反演和杜教筛综合。

[Luogu P2424] 约数和

Portal.

将 $f$ 做前缀和之后就可以直接使用数论分块计算。

查看代码

#include <iostream>
#include <cstdio>

using namespace std;
typedef long long i64;

i64 H(int n) {
    i64 res = 0, l = 1, r;
    while (l <= n) {
        r = n / (n / l);
        res += (l + r) * (r - l + 1) * (n / l) / 2;
        l = r + 1;
    }
    return res;
}

int main(void) {
    int x, y; cin >> x >> y;
    cout << H(y) - H(x - 1) << "\n";
    return 0;
}

[CQOI2007] 余数求和

Portal.

将取模拆开即可。

查看代码

#include <iostream>
#include <cstdio>

using namespace std;
typedef long long i64;

int main(void) {
    i64 n, k;
    scanf("%lld%lld", &n, &k);
    i64 ans = n * k; 
    for (i64 l = 1, r; l <= n; l = r + 1) {
        if (k / l != 0) r = min(n, k / (k / l));
        else break;
        ans -= (l + r) * (r - l + 1) * (k / l) / 2;
    }
    printf("%lld\n", ans);
    return 0;
}

[CTT2012] 模积和

Portal.

直接数论分块即可，需要将式子拆开。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 
const int P = 19940417; 

i64 n, m; 
i64 a, b, c;

i64 calc(int n) {
    return 3323403ll * n % P * (n + 1) % P * (n * 2 % P + 1) % P; 
}
i64 calc(int l, int r) {
    return calc(r) - calc(l - 1); 
}

int main(void) {
    cin >> n >> m; if (n > m) swap(n, m); 
    a = n * n % P; b = m * m % P; c = n * m % P * n % P; 
    for (int l = 1, r = 0; l <= n; l = r + 1) {
        r = n / (n / l); 
        a = (a - 1ll * (l + r) * (r - l + 1) / 2 % P * (n / l) % P) % P; 
    }
    for (int l = 1, r = 0; l <= m; l = r + 1) {
        r = m / (m / l); 
        b = (b - 1ll * (l + r) * (r - l + 1) / 2 % P * (m / l) % P) % P; 
    }
    for (int l = 1, r = 0; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l)); 
        i64 len = 1ll * (l + r) * (r - l + 1) / 2 % P; 
        c = (c + calc(l, r) * (n / l) % P * (m / l) % P - m * len % P * (n / l) % P - n * len % P * (m / l) % P) % P; 
    }
    cout << ((a * b % P - c) % P + P) % P << "\n";  
    return 0; 
}

莫比乌斯反演

带有特殊技巧的题目会标注。

莫比乌斯反演往往只能去掉艾弗森括号（虽然莫比乌斯函数自身还有容斥的性质），因此重要的反而是推式子。

[HAOI2011] Problem b

Portal.

利用二维差分进行计算，而且发现只有 $k$ 的倍数有用，因此可以将所有数除以 $k$ 转变为求互质。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 

int mu[50005], prime[50005], tot;
bool v[50005]; 

i64 calc(int x, int y) {
    i64 res = 0; 
    for (int l = 1, r; l <= x && l <= y; l = r + 1) {
        r = min(x / (x / l), y / (y / l));
        res += 1ll * (mu[r] - mu[l - 1]) * (x / l) * (y / l); 
    }
    return res; 
}

int main(void) {
    mu[1] = 1; 
    for (int i = 2; i <= 50000; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1; 
        for (int j = 1; j <= tot && i * prime[j] <= 50000; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break; 
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 2; i <= 50000; ++i) mu[i] += mu[i - 1]; 
    int T; scanf("%d", &T); while (T--) {
        int a, b, c, d, k; scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        a = (a - 1) / k, c = (c - 1) / k, b /= k, d /= k; 
        printf("%d\n", calc(a, c) - calc(a, d) - calc(b, c) + calc(b, d));
    }
    return 0;
}

欧拉反演 | [NOI2010] 能量采集

Portal.

要求的显然是：

$\sum_{i=1}^n\sum_{i=1}^m 2\times \gcd(i,j)-1$

然后欧拉反演直接拍上去。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 

int n, m; 
int prime[50005], tot;
i64 phi[100005];

int main(void) {
    scanf("%d%d", &n, &m);
    if (m > n) swap(n, m);
    phi[1] = 1; 
    for (int i = 2; i <= n; ++i) {
        if (!phi[i]) prime[++tot] = i, phi[i] = i - 1; 
        for (int j = 1; j <= tot && i * prime[j] <= n; ++j) {
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j]; 
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1); 
        }
    }
    for (int i = 1; i <= n; ++i) phi[i] += phi[i - 1]; 
    i64 ans = 0; 
    for (int l = 1, r; l <= m; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += 1ll * (phi[r] - phi[l - 1]) * (n / l) * (m / l); 
    }
    printf("%lld\n", 2 * ans - 1ll * n * m);
    return 0;
}

莫比乌斯函数与容斥原理 | [Luogu P4318] 完全平方数

Portal.

设 $f(n)$ 代表 $1\sim n$ 中无平方因子的数的个数，然后二分答案，找到 $f(r)\ge K$ 的最小 $r$ 。现在考虑如何计算 $f$ 。

首先去掉 $2^2,3^2,\cdots,p^2$ 的倍数，但是其中有的会算入两次，所以对于 $6^2$ 这种要减去。而像 $4^2$ 这种在 $2^2$ 时已经算完了，不能统计。发现这就是莫比乌斯函数，因此莫比乌斯函数是在对 $\mathbb{N}$ 做容斥。答案是：

$f(n)=\sum_{i}\mu(i)\left\lfloor\frac n {i^2}\right\rfloor$

查看代码

#include <bits/stdc++.h>
using namespace std;
const int N = 200000;

int mu[200005], prime[100005], tot; 
bool v[200005];

int calc(int n) {
    int res = 0; 
    for (int i = 1; i * i <= n; ++i)
        res += mu[i] * (n / (i * i));
    return res; 
}

int main(void) {
    mu[1] = 1; 
    for (int i = 2; i <= N; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1; 
        for (int j = 1; j <= tot && i * prime[j] <= N; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i]; 
        }
    }
    int T, n; scanf("%d", &T); 
    while (T--) {
        scanf("%d", &n); 
        int L = 0, R = 2000000001; 
        while (L + 1 != R) {
            int mid = L + (R - L) / 2; 
            if (calc(mid) >= n) R = mid;
            else L = mid;
        }
        printf("%d\n", R);
    }
    return 0;
}

[Luogu P2257] YY 的 GCD

Portal.

根之前的推法是一样的，不过外面要多枚举一个 $p$ 。

推完之后有这么一个东西：

$\sum_{p\mid k}\mu\left(\frac k p\right)$

这个东西可以直接在 $O(n\log \log n)$ 的时间预处理。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 
const int N = 10000000; 

int mu[10000005], prime[4000005], tot; 
i64 f[10000005];
bool v[10000005]; 

int main(void) {
    mu[1] = 1; 
    for (int i = 2; i <= N; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && i * prime[j] <= N; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= tot; ++i)
        for (int j = 1; j * prime[i] <= N; ++j)
            f[j * prime[i]] += mu[j];
    for (int i = 2; i <= N; ++i) f[i] += f[i - 1];
    int T; scanf("%d", &T); while (T--) {
        int n, m; scanf("%d%d", &n, &m);
        if (n > m) swap(n, m);
        i64 ans = 0; 
        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans += (f[r] - f[l - 1]) * (n / l) * (m / l);
        } printf("%lld\n", ans); 
    }
    return 0;
}

[SDOI2015] 约数个数和

Portal.

关键是：

$\tau(ij)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]$

然后随便推就可以了。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 
const int N = 50000; 

int prime[20005], tot; 
bool v[50005]; 
int mu[50005], g[50005]; 

void init(void) {
    mu[1] = 1; 
    for (int i = 2; i <= N; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1; 
        for (int j = 1; j <= tot && i * prime[j] <= N; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i]; 
        }
    }
    for (int i = 1; i <= N; ++i) mu[i] += mu[i - 1]; 
    for (int i = 1; i <= N; ++i) {
        i64 res = 0; 
        for (int l = 1, r; l <= i; l = r + 1) {
            r =  i / (i / l); 
            res += 1ll * (r - l + 1) * (i / l); 
        }
        g[i] = res; 
    }
}

int main(void) {
    int T, n, m; scanf("%d", &T); init(); 
    while (T--) {
        scanf("%d%d", &n, &m); 
        if (n > m) swap(n, m); 
        i64 ans = 0; 
        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l)); 
            ans += 1ll * (mu[r] - mu[l - 1]) * g[n / l] * g[m / l]; 
        } printf("%lld\n", ans); 
    }
    return 0;
}

[国家集训队] Crash 的数字表格

Portal.

来推式子！

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

设 $F(n,m)=\sum_{i=1}^n\sum_{j=1}^m ij[\gcd(i,j)=1]$ ，则原式为 $\sum_{d=1}^n d\times F(\lfloor \frac n d\rfloor, \lfloor \frac m d \rfloor)$ 。

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

后面那坨也可以通过数论分块计算。

查看代码

#include <bits/stdc++.h>
using namespace std;
const int P = 20101009; 
const int N = 10000000; 

int n, m, mu[10000005], prime[4000005], tot; 
bool v[10000005]; 

int F(int n) { return 1ll * n * (n + 1) / 2 % P; }
int sum(int n, int m) {
    int res = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l)); 
        res = (res + 1ll * (mu[r] - mu[l - 1]) * F(n / l) % P * F(m / l) % P) % P; 
    } 
    return res; 
}

int main(void) {
    mu[1] = 1; 
    for (int i = 2; i <= N; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1; 
        for (int j = 1; j <= tot && i * prime[j] <= N; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i]; 
        }
    }
    for (int i = 2; i <= N; ++i) mu[i] = (mu[i - 1] + 1ll * i * i % P * mu[i] % P) % P;
    scanf("%d%d", &n, &m); if (n > m) swap(n, m); 
    int ans = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l)); 
        ans = (ans + 1ll * (l + r) * (r - l + 1) / 2 % P * sum(n / l, m / l) % P) % P; 
    }
    printf("%d\n", (ans + P) % P); return 0;
}

[CF1780F] Three Chairs

Portal.

将 $a$ 排序，要求的就是：

$\begin{aligned} ans&=\sum_{1\le i<j\le n} (j-i+1)[\gcd(a_i, a_j)=1]\\ &=\sum_{1\le i<j\le n} (j-i+1)\sum_{d\mid a_i,d\mid a_j}\mu(d)\\ &=\sum_{d}\mu(d)\sum_{1\le i<j\le k} p_j-p_i-1 \end{aligned}$

其中 $p$ 表示是 $d$ 的倍数的数的下标位置。简单转化一下贡献，后面那一坨就是：

$-\frac{k(k-1)} 2+\sum_{i=1}^{k}p_i\times (2i-k-1)$

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 
const int N = 300000; 

int n; 
int a[300005], mu[300005];
int prime[100005], tot, pos[300005], idx[300005]; 
bool v[300005]; 

int main(void) {
    scanf("%d", &n); 
    for (int i = 1; i <= n; ++i) scanf("%d", a + i); 
    sort(a + 1, a + n + 1); mu[1] = 1; 
    for (int i = 1; i <= n; ++i) pos[a[i]] = i; 
    for (int i = 2; i <= N; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && i * prime[j] <= N; ++j) {
            v[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i]; 
        }
    }
    i64 ans = 0; 
    for (int d = 1; d <= N; ++d) if (mu[d]) {
        int k = 0; 
        for (int i = d; i <= N; i += d) 
            if (pos[i]) idx[++k] = pos[i];
        i64 res = -1ll * k * (k - 1) / 2; 
        for (int i = 1; i <= k; ++i)
            res += 1ll * idx[i] * (2 * i - k - 1);
        ans += mu[d] * res; 
    }
    printf("%lld\n", ans); 
    return 0;
}

[Luogu P4844] LJJ 爱数数

Portal.

方程两边同乘 $abc$ 可以得到 $ac+bc=ab\Rightarrow ab-ac-bc+c^2=c^2\Rightarrow (a-c)(b-c)=c^2$ 。

设 $a-c=x,b-c=y$ ，则 $xy=c^2$ ，可以证明 $a\perp b\perp c\iff x\perp y$ 。

需要满足 $\max\{x,y\}+c\le n$ ，设 $x>y$ ，则有 $x+c\le n$ 。最后让答案 $\times 2$ ，再加上一个 $a=b=2,c=1$ 即可。

由于 $x\perp y$ ，因此 $c^2$ 的每一种质因子只能都给 $x$ 或都给 $y$ ，这样有 $x=i^2,y=j^2(i>j)$ ，因此 $i^2+ij\le n\Rightarrow j\le \frac n i - i$ ，设 $R_i=\min\{i-1,\frac n i -1\}$ ，则答案为 $\sum_{i=1}^{n}\sum_{j=1}^{R_i}[\gcd(i,j)=1]$ ，莫比乌斯反演碾过去即可，时间复杂度为 $O(\sqrt{n}\log \sqrt{n})$ 。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef long long i64; 

i64 n, R[1000005]; int m; 
int mu[1000005], prime[400005], tot; 

bool v[1000005];

void init(int n) {
    mu[1] = 1; 
    for (int i = 2; i <= n; ++i) {
        if (!v[i]) prime[++tot] = i, mu[i] = -1; 
        for (int j = 1; j <= tot && i * prime[j] <= n; ++j) {
            v[i * prime[j]] = 1; 
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i]; 
        }
    }
}

int main(void) {
    cin >> n; if (n == 1) return puts("0"), 0; 
    for (int i = 2; i <= n; ++i) {
        R[i] = n / i - i; if (i - 1 < R[i]) R[i] = i - 1; 
        if (R[i] <= 0) break; 
        m = i; 
    } init(m);
    i64 ans = 0; 
    for (int d = 1; d <= m; ++d) if (mu[d]) {
        for (int i = d; i <= m; i += d)
            ans += mu[d] * (R[i] / d); 
    }
    printf("%lld\n", ans * 2 + 1); return 0;
}