[学习笔记]初等数论
最大公约数 \(gcd\)
欧几里得算法（辗转相除法）：
\[\gcd(a, b) = \gcd(b, a \bmod b) \]代码：
int gcd(int a, int b){
return b ? gcd(b, a%b) : a;
}
或者直接使用 __gcd(a, b)。
辗转相减法：
\[\gcd(a, b) = \gcd(a, b-a) \]推广到 \(n\) 项：
\[\gcd(a_1, a_2, \dots, a_n) = \gcd(a_1, a_2-a_1, a_3-a_2,\dots, a_n-a_{n-1}) \]注意：因为 \(gcd\) 不可能是负数，所以当权值为负时答案应该取绝对值。
线性筛
欧拉筛。不多讲了。
int n;
int prime[N], cnt;
bool is_prime[N];
void getprime(){
for(int i=2; i<=n; i++) is_prime[i] = 1;
for(int i=2; i<=n; i++){
if(is_prime[i]){
prime[++cnt] = i;
}
for(int j=1; j<=cnt && i*prime[j]<=n; j++){
is_prime[i*prime[j]] = 0;
if(i%prime[j] == 0) break;
}
}
}
积性函数
数论函数：
数论函数指定义域为正整数的函数。数论函数也可以视作一个数列。
积性函数：
若一个数论函数 \(f\) 满足，对所有互质的 \(n, m\)，都满足 \(f(nm) = f(n)f(m)\)，则称 \(f\) 是一个积性函数。
此外，如果对所有 \(n, m\)，都满足 \(f(nm) = f(n)f(m)\)，则称 \(f\) 为完全积性函数。
性质：都可以通过筛法来快速求解。
常见的积性函数：
欧拉函数 \(\varphi(n)\)：
表示 \(\left[ 1, n \right]\) 中与 \(n\) 互质的数的数量。
欧拉公式：记 \(n = \prod p_i^{c_i}\)，则 \(\varphi(n) = \prod (p_i-1)p_i^{c_i-1}\)。
性质：
\(\sum\limits_{d \mid n} \varphi(d) = n\)
莫比乌斯函数 \(\mu(n)\)：
\(\mu(n) = \begin{cases} (-1)^k & n = p_1 \hspace{0.1em} p_2 \dots p_k \hspace{0.3em} (\hspace{0.1em}p_i\hspace{0.1em}是质数) \\ 0 & n \hspace{0.1em} 的某个质因子指数 \hspace{0.1em} c \ge 2 \end{cases}\)
性质：
\(\sum\limits_{d \mid n} \mu(d) = \left[n = 1\right]\)
单位函数 \(\epsilon(n)\)：
\(\epsilon(n) = \begin{cases} 1 & n = 1 \\ 0 & n \neq 1 \end{cases}\)
筛 \(\varphi(n)\)
从欧拉筛出发考虑。下面的 \(i\) 指当前的数，\(p\) 指质数表中的质数。
-
当 \(i\) 是质数时，显然 \(\varphi(i) = i-1\)。
-
当 \(i \not\equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(i, p\) 互质时，根据积性函数定义得：\(\varphi(ip) = \varphi(i)\varphi(p) = \varphi(i) \times (p-1)\)。
-
当 \(i \equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(p\) 是 \(i\) 的一个因子时。设 \(i\) 唯一分解中 \(p\) 的指数是 \(c\)，因为 \(\varphi(p^c) = (p-1)p^{c-1}\)，因此 \(\varphi(ip) = \varphi(i) \times p\)。
解释：
-
对于 \(\varphi(p^c) = (p-1)p^{c-1}\)，\(\left[ 1, p\right]\) 内有 \(p-1\) 个数与其互质。所以 \(\left[ 1, p^c\right]\) 内有 \((p-1)p^{c-1}\) 个数与其互质。
-
对于 \(\varphi(ip) = \varphi(i) \times p\)。因为 \(\varphi(ip) = \varphi(\frac{i}{p^c}) \cdot \varphi(p^{c+1})\)，又因为 \(\frac{i}{p^c}\) 与 \(p^c\) 互质，所以 \(\varphi(\frac{i}{p^c}) \varphi(p^c) = \varphi(i)\)。又因为 \(\varphi(p^{c+1}) = (p-1)p^{c}\)。所以将 \(\varphi(\frac{i}{p^c})\) 和 \(\varphi(p^{c+1})\) 两项拆开消元后即可得到 \(\varphi(ip) = \varphi(i) \times p\)。
也可通过欧拉公式（欧拉函数的定义）理解。
-
代码实现：
int n;
int prime[N], phi[N], cnt;
bool is_prime[N];
void getprime(){
phi[1] = 1;
for(int i=2; i<=n; i++) is_prime[i] = 1;
for(int i=2; i<=n; i++){
if(is_prime[i]){
phi[i] = i-1;
prime[++cnt] = i;
}
for(int j=1; j<=cnt && i*prime[j]<=n; j++){
is_prime[i*prime[j]] = 0;
if(i%prime[j] == 0){
phi[i*prime[j]] = prime[j]*phi[i];
break;
}
else phi[i*prime[j]] = (prime[j]-1)*phi[i];
}
}
}
筛 \(\mu(n)\)
还是考虑三种情况。
- 当 \(i\) 是质数时，显然 \(\mu(i) = -1\)。
- 当 \(i \not\equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(i, p\) 互质时，根据积性函数定义得：\(\mu(ip) = \mu(i)\mu(p) = -\mu(i)\)。
- 当 \(i \equiv 0 \hspace{0.3em} (\bmod \hspace{0.3em} p)\) 时，即 \(p\) 是 \(i\) 的一个因子时。根据该函数定义得 \(\mu(ip) = 0\)。
代码实现应该差不多。
P2568 GCD
首先我们枚举质数：
\[\sum_{p \in \rm prime} \sum_{i=1}^{n} \sum_{j=1}^{n} \left[ \gcd(i, j) = p \right] \]变形：
\[\sum_{p \in \rm prime} \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{p} \rfloor} \left[ \gcd(i, j) = 1 \right] \]接下来改变 \(j\) 的枚举上界（其中 \(-1\) 的原因是 \(i=j=1\) 时的答案会被重复统计，因此这里的 \(-1\) 是在第二层中的）：
\[\sum_{p \in \rm prime} \left( \sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \left( 2\sum_{j=1}^{i} \left[ \gcd(i, j) = 1 \right] \right) -1 \right) \]此时发现 \(\sum\limits_{j=1}^{i} \left[ \gcd(i, j) = 1 \right]\) 就是 \(\varphi(i)\)，故原式可化为：
\[\sum_{p \in \rm prime} \left( 2\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor} \varphi(i) -1 \right) \]可以前缀和求 \(\sum \varphi(i)\) 的值，最后枚举 \(p \in \rm prime\) \({\rm and}\) \(p \le n\) 并统计答案即可。
时间复杂度：\(O(n)\)。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e7 + 5;
int prime[N], phi[N], cnt;
bool is_prime[N];
ll ans, sum[N];
int main(){
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
int n; cin>>n;
phi[1] = 1;
for(int i=2; i<=n; i++) is_prime[i] = 1;
for(int i=2; i<=n; i++){
if(is_prime[i]){
phi[i] = i-1;
prime[++cnt] = i;
}
for(int j=1; j<=cnt && i*prime[j]<=n; j++){
is_prime[i*prime[j]] = 0;
if(i%prime[j] == 0){
phi[i*prime[j]] = prime[j]*phi[i];
break;
}
else phi[i*prime[j]] = (prime[j]-1)*phi[i];
}
}
for(int i=1; i<=n; i++) sum[i] = sum[i-1]+phi[i];
for(int i=1; i<=cnt; i++) ans += 2*sum[n/prime[i]]-1;
cout<<ans;
return 0;
}
P10463 Interval GCD
辗转相减法（前置知识）：
\[\gcd(a, b) = \gcd(a, b-a) \]推广到 \(n\) 项：
\[\gcd(a_1, a_2, \dots, a_n) = \gcd(a_1, a_2-a_1, a_3-a_2,\dots, a_n-a_{n-1}) \]看到题意中把 \(A_l, A_{l+1}, \dots, A_r\) 都加上 \(d\)，自然想到差分。
令差分完之后的数组为 \(B\)，那么 \(\gcd(B_{l+1}, B_{l+2}, \dots, B_r) = \gcd(A_{l+1}, A_{l+2}, \dots A_r)\)，而 \(\gcd(B_l,B_{l+1}) \neq \gcd(A_l, A_{l+1})\)，所以 \(A_l\) 需要单独维护。
以下用线段树维护差分数组的 \(\gcd\)，用树状数组维护差分数组，最后通过前缀和求得 \(A_l\)。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 5e5+5;
#define lowbit(x) x&(-x)
#define ls(x) ((x)<<1)
#define rs(x) ((x)<<1|1)
#define gcd(x, y) abs(__gcd((x), (y)))
int n, m;
ll a[N], cha[N], c[N];
struct node{
int l, r;
ll v;
}tr[N<<2];
void add(int x, ll y){
while(x <= n){
c[x] += y;
x += lowbit(x);
}
}
ll sum(int x){
ll ret = 0;
while(x){
ret += c[x];
x -= lowbit(x);
}
return ret;
}
void pushup(int x){
tr[x].v = gcd(tr[ls(x)].v, tr[rs(x)].v);
}
void buildtr(int x, int l, int r){
tr[x].l = l, tr[x].r = r;
if(l == r){
tr[x].v = cha[l];
return;
}
int mid = (l+r) >> 1;
buildtr(ls(x), l, mid);
buildtr(rs(x), mid+1, r);
pushup(x);
}
void upd(int x, int tar, ll d){
if(tr[x].l > tar || tr[x].r < tar)
return;
if(tr[x].l == tr[x].r && tr[x].l == tar){
tr[x].v += d;
return;
}
int mid = (tr[x].l+tr[x].r) >> 1;
if(tar<=mid) upd(ls(x), tar, d);
upd(rs(x), tar, d);
pushup(x);
}
ll query(int x, int l, int r){
if(tr[x].l > r || tr[x].r < l)
return 0;
if(tr[x].l >= l && tr[x].r <= r)
return tr[x].v;
return gcd(query(ls(x), l, r), query(rs(x), l, r));
}
int main(){
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
cin>>n>>m;
for(int i=1; i<=n; i++){
cin>>a[i];
cha[i] = a[i] - a[i-1];
add(i, cha[i]);
}
buildtr(1, 1, n);
while(m--){
char op; int l, r; cin>>op>>l>>r;
if(op == 'C'){
ll d; cin>>d;
upd(1, l, d);
upd(1, r+1, -d);
add(l, d);
add(r+1, -d);
} else if(op == 'Q'){
cout<<gcd(sum(l), query(1, l+1, r))<<"\n";
}
}
return 0;
}
USACO08DEC Patting Heads S
从枚举倍数的角度思考：对于一个数 \(i\)，若它在原数组中出现了 \(cnt_i\) 次，那么 \(i\) 这个数会对它的每一个倍数产生 \(cnt_i\) 的贡献。这样就可以通过查询这样产生的贡献的总和来计算答案了。
注意最后要减去一个对自己的贡献。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e6+5;
int a[N], w[N], c[N];
int main(){
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
int m = 0; // 值域
int n; cin>>n;
for(int i=1; i<=n; i++){
cin>>a[i];
c[a[i]]++;
m = max(m, a[i]);
}
for(int i=1; i<=m; i++){
for(int j=i; j<=m; j+=i){
w[j] += c[i];
}
}
for(int i=1; i<=n; i++)
cout<<w[a[i]]-1<<"\n"; // 减去自己对自己的贡献
return 0;
}
Divan and Kostomuksha (easy version)
考虑 \(dp_i\) 表示目前前缀 \(\gcd = i\)，前缀 \(\gcd\) 和的最大值。
首先计算出序列 \(a\) 中是 \(i\) 的倍数的数量，记为 \(cnt_i\)。（这里处理方法和上一题相似）
接下来考虑转移，假设是从 \(j\) 转移，转移到 \(i\)，相当于在 \(j+1\) 至 \(f_j\) 所在位置全部填上 \(a\) 中是 \(i\) 的倍数但不是 \(j\) 的倍数的数，容易发现多的是 \(i\) 的倍数的数量为 \(cnt_i - cnt_j\)。所以转移式为：
\[dp_i = \max_{i \mid j}(dp_j + (cnt_i - cnt_j) \times i) \]#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 5e6;
ll cnt[N+5], a[N+5];
ll dp[N+5];
//cnt 表示 a 中能被 i 整除的数的个数
int main(){
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll n, m = 0; cin>>n; // m 是值域
for(int i=1; i<=n; i++){
cin>>a[i];
cnt[a[i]]++;
m = max(m, a[i]);
}
for(int i=1; i<=m; i++){
for(int j=i*2; j<=m; j+=i){
cnt[i] += cnt[j];
}
}
ll ans = 0;
for(int i=m; i>=1; i--){
dp[i] = cnt[i]*i;
for(int j=i*2; j<=m; j+=i){
dp[i] = max(dp[i], dp[j]+1ll*(cnt[i]-cnt[j])*i);
}
ans = max(ans, dp[i]);
}
cout<<ans;
return 0;
}