【数论】卷积反演大集合

不知道为啥脑抽要学数论，骂声一片中发现数论还没入门（悲）。

1. 狄利克雷卷积与数论函数

1.1 数论函数

定义：数论函数为值域为整数的函数。

简单数论函数：

• \(I(n)\)，恒等函数，恒等为 \(1\)。
• \(e(n)\)，元函数，卷积中的单位元，若 \(n=1\)，\(e(n)=1\)。否则为 \(e(n)=0\)。
• \(id(n)\)，单位函数，\(id(n)=n\)。
• \(id^x(n)\)，幂函数，\(id^x(n)=n^x\)。

在数论函数的基础上，有以下。

定义：完全积性函数，\(f(1)=1\)，且对于任意的整数 \(a\) 和 \(b\) 有 \(f(ab)=f(a)f(b)\)。

\(e(n),I(n),id(n)\) 是完全积性函数。

定义：积性函数，\(f(1)=1\)，且对于任意互质的 \(a\) 和 \(b\) 有 \(f(ab)=f(a)f(b)\)。

欧拉函数与莫比乌斯函数皆为积性函数。证明显然。

1.2 狄利克雷卷积

定义：\((f * g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d})\)。

或 \((f * g)(n)=\sum\limits_{xy=n}f(x)g(y)\)。（优美的对称结构）

交换律：\(f*g=g*f\)

结合律：\((f*g)*h=f*(g*h)\)

会用就行。

2. 莫比乌斯反演

2.1 莫比乌斯函数

函数 \(f\) 与 \(F\) 满足如下关系：

\(F(n)=\sum\limits_{d|n}f(d)\)

由狄利克雷卷积得：\(F(n)=\sum\limits_{d|n}(1*f(d))\)

于是乎就有：\(F=I*f\)

若已知 \(f\) 求 \(F\)，则需反演。

变形：\(f=F*I^{-1}\)

这里的 \(I^{-1}\) 即为莫比乌斯函数，记为 \(\mu\)。

积性函数的逆还是积性函数，所以 \(\mu\) 为积性函数。

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

2.2 莫比乌斯反演

• 嵌套式反演公式：由 \(\mu*I=e\) 得 \(\sum\limits_{d|n}\mu(d)=[n=1]\)。

好了，你现在已经数论入门了！！。

3. 欧拉反演

3.1 欧拉函数

记 \(\varphi(n)\) 表示小于 \(n\) 的正整数中与 \(n\) 互质的数的个数。

计算方法：将总数减去与 \(n\) 不互质的数的个数。将 \(n\) 质因数分解后，记 \(p_1,p_2,p_3,\dots,p_k\) 为 \(n\) 的质因子，则有：

\[\varphi(n)=n-\sum\limits_{i=1}^{n}\frac{n}{p_i}+\sum\limits_{1\le i<j \le n}\frac{n}{p_ip_j}-\sum\limits_{1\le i<j<k \le n}\frac{n}{p_ip_jp_k}+\dots+(-1)^{k+1}\frac{n}{p_1p_2p_3\dots p_k} \]

因数分解可得：

\[\varphi(n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})(1-\frac{1}{p_3})\dots(1-\frac{1}{p_k}) \]

计算欧拉函数的步骤中运用到了容斥原理，可见容斥原理的应用范围非常广。

欧拉函数的性质：

1. 若 \(\gcd(a,b)=1\)，则有 \(\varphi(a)·\varphi(b)=\varphi(ab)\)，即欧拉函数为积性函数。

   证明：

   设 \(A\) 为 \(a\) 的质因子集合，集合元素记为 \(p_1,p_2,p_3,\dots,p_k\)，\(B\) 为 \(b\) 的质因子集合，集合元素记为 \(P_1,P_2,P_3,\dots,P_m\)。

   \(\because\gcd(a,b)=1\)，\(\therefore A \cap B=\varnothing\).

   \(\because\varphi(a)=a(1-\frac{1}{p_1})(1-\frac{1}{p_2})(1-\frac{1}{p_3})\dots(1-\frac{1}{p_k}),\varphi(b)=b(1-\frac{1}{P_1})(1-\frac{1}{P_2})(1-\frac{1}{P_3})\dots(1-\frac{1}{P_m})\)

   \(\therefore \varphi(ab)=ab(1-\frac{1}{p_1})(1-\frac{1}{p_2})(1-\frac{1}{p_3})\dots(1-\frac{1}{p_k})(1-\frac{1}{P_1})(1-\frac{1}{P_2})(1-\frac{1}{P_3})\dots(1-\frac{1}{P_m})\)

   \(\therefore\varphi(a)·\varphi(b)=\varphi(ab)\)

2. 若 \(p\) 为质数，\(\varphi(p)=p-1\).

   证明：

   \(\because p\) 为质数，\(\therefore\) 小于 \(p\) 的正整数中与 \(p\) 互质的数的个数为 \(p-1\).

3. 若 \(p\) 为质数，\(\varphi(p)=p^k-p^{k-1}\).

   证明：

   \(\because p\) 为质数，\(\therefore\) 与 \(p\) 不互质的数只能为 \(p\) 的倍数，共有 \(p^{k-1}\) 个。

   \(\therefore\) \(\varphi(p)=p^k-p^{k-1}\).

不用怀疑了，全从数学大礼包上贺的。

3.2 欧拉反演

• 嵌套式欧拉反演：由 \(\varphi * I = id\) 得 \(\sum\limits_{d|n}\varphi(d)=n\)

4. 反演习题

P3455 [POI2007] ZAP-Queries

Problem

给定 \(a, b, c\)，求 \(\sum\limits_{i=1}^a\sum\limits_{j=1}^b[(i,j)=c]\).

Solve

莫反好题。

转化形式：

\[\sum\limits_{i=1}^{\lfloor\frac{a}{c}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{c}\rfloor}[(i,j)=1] \]

代入嵌入式莫比乌斯反演：

\[\sum\limits_{i=1}^{\lfloor\frac{a}{c}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{c}\rfloor}\sum\limits_{d|(i,j)}\mu(d) \]

由 \(d|(i,j)\Rightarrow d|i,d|j\) 可知：

\[\sum\limits_{i=1}^{\lfloor\frac{a}{c}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{c}\rfloor}\sum\limits_{d|i,d|j}\mu(d) \]

交换枚举顺序：

\[\sum\limits_{d=1}^{\min(\lfloor\frac{a}{c}\rfloor,\lfloor\frac{b}{c}\rfloor)}\mu(d)\sum\limits_{d|i}^{\lfloor\frac{a}{c}\rfloor}\sum\limits_{d|j}^{\lfloor\frac{b}{c}\rfloor}1 \]

化简：

\[\sum\limits_{d=1}^{\min(\lfloor\frac{a}{c}\rfloor,\lfloor\frac{b}{c}\rfloor)}\mu(d)\lfloor\frac{\frac{a}{c}}{d}\rfloor\lfloor\frac{\frac{b}{c}}{d}\rfloor \]

\(\mu(d)\) 可以用前缀和处理，\(\lfloor\frac{\frac{a}{c}}{d}\rfloor\lfloor\frac{\frac{b}{c}}{d}\rfloor\) 整除分块即可。

Code

点击查看代码

#include <bits/stdc++.h>
#define int long long
#define H 19260817
#define rint register int
#define For(i,l,r) for(rint i=l;i<=r;++i)
#define FOR(i,r,l) for(rint i=r;i>=l;--i)
#define MOD 1000003
#define mod 1000000007

using namespace std;

inline int read() {
  rint x=0,f=1;char ch=getchar();
  while(ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
  while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
  return x*f;
}

void print(int x){
  if(x<0){putchar('-');x=-x;}
  if(x>9){print(x/10);putchar(x%10+'0');}
  else putchar(x+'0');
  return;
}

const int N = 5e4 + 10;

int T, a, b, d, prime[N], tot, mu[N], sum[N];

bool vis[N];

void Euler() {
  vis[1] = mu[1] = 1;
  for (int i = 2; i <= N-10; i++) {
    if(!vis[i]) prime[++tot] = i, mu[i] = -1;
    for (int j = 1; i * prime[j] <= N-10 && j <= tot; j++) {
      vis[i * prime[j]] = 1;
      mu[i * prime[j]] = (i % prime[j] == 0 ? 0 : -mu[i]);
      if(i % prime[j] == 0) break;
    } 
  }
}

signed main() {
  T = read();
  Euler();
  For(i,1,N) {
    sum[i] = sum[i-1] + mu[i];
  }
  while(T--) {
    a = read(), b = read(), d = read();
    int ans = 0, n = a / d, m = b / d;
    for (int l = 1, r = 0; l <= min(n, m); l = r + 1) {
      r = min(n / (n / l), m / (m / l));
      ans += (sum[r] - sum[l-1]) * (n / l) * (m / l);
    }
    cout << ans << '\n';
  }
  return 0;
}

P2522 [HAOI2011] Problem b

Problem

给定 \(a, b, c, d, k\)，求 \(\sum\limits_{i=a}^b\sum\limits_{i=c}^d[(i,j)=k]\).

Solve

根据上一题得到启发：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[(i,j)=k]=\sum\limits_{d=1}^{\min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}\mu(d)\lfloor\frac{\frac{n}{k}}{d}\rfloor\lfloor\frac{\frac{m}{k}}{d}\rfloor \]

若上题求得是一个 \(n \times m\) 的 \((i,j)\) 矩阵，则此题求以 \((a,c),(b,d)\) 为顶点的子矩阵。

在上一题的基础上二维前缀和求解即可。

Code

点击查看代码

#include <bits/stdc++.h>
#define int long long
#define H 19260817
#define rint register int
#define For(i,l,r) for(rint i=l;i<=r;++i)
#define FOR(i,r,l) for(rint i=r;i>=l;--i)
#define MOD 1000003
#define mod 1000000007

using namespace std;

inline int read() {
  rint x=0,f=1;char ch=getchar();
  while(ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
  while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
  return x*f;
}

void print(int x){
  if(x<0){putchar('-');x=-x;}
  if(x>9){print(x/10);putchar(x%10+'0');}
  else putchar(x+'0');
  return;
}

const int N = 5e4 + 10;

int T, a, b, c, d, k, prime[N], tot, mu[N], sum[N];

bool vis[N];

void Euler() {
  vis[1] = mu[1] = 1;
  for (int i = 2; i <= N-10; i++) {
    if(!vis[i]) prime[++tot] = i, mu[i] = -1;
    for (int j = 1; i * prime[j] <= N-10 && j <= tot; j++) {
      vis[i * prime[j]] = 1;
      mu[i * prime[j]] = (i % prime[j] == 0 ? 0 : -mu[i]);
      if(i % prime[j] == 0) break;
    } 
  }
}

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

signed main() {
  T = read();
  Euler();
  For(i,1,N) {
    sum[i] = sum[i-1] + mu[i];
  }
  while(T--) {
    a = read(), b = read(), c = read(), d = read(), k = read();
    cout << calc(b, d) - calc(a-1, d) - calc(b, c-1) + calc(a-1, c-1) << '\n';
  }
  return 0;
}

P2158 [SDOI2008] 仪仗队

Problem

给定一个 \(n\times n\) 的点矩阵，求位于左下角能看见多少个不被遮挡的点。

Solve

若 \((x,y)\) 能被看见，则满足 \((i,j)=1\)。

若存在点 \((kx,ky)\) 则它一定会被 \((x,y)\) 挡住。

所求即为 \(\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}[(i,j)=1]\).

莫比乌斯反演求解。

注意细节：计算时按 \(2+\sum\limits_{i=1}^{n-1}\sum\limits_{j=1}^{n-1}[(i,j)=1]\) 计算。\(n=1\) 需要特判。

Code

点击查看代码

#include <bits/stdc++.h>
#define int long long
#define H 19260817
#define rint register int
#define For(i,l,r) for(rint i=l;i<=r;++i)
#define FOR(i,r,l) for(rint i=r;i>=l;--i)
#define MOD 1000003
#define mod 1000000007

using namespace std;

inline int read() {
  rint x=0,f=1;char ch=getchar();
  while(ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
  while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
  return x*f;
}

void print(int x){
  if(x<0){putchar('-');x=-x;}
  if(x>9){print(x/10);putchar(x%10+'0');}
  else putchar(x+'0');
  return;
}

const int N = 4e4 + 10;

int n, prime[N], tot, mu[N], sum[N];

bool vis[N];

void Eular() {
  vis[1] = mu[1] = 1;
  for (int i = 2; i <= 4e4; i++) {
    if(!vis[i]) prime[++tot] = i, mu[i] = -1;
    for (int j = 1; i * prime[j] <= 4e4 && j <= tot; j++) {
      vis[i * prime[j]] = 1;
      mu[i * prime[j]] = (i % prime[j] == 0 ? 0 : -mu[i]); 
      if(i % prime[j] == 0) break;
    }
  }
}

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

signed main() {
  n = read();
  Eular();
  For(i,1,n) sum[i] = sum[i-1] + mu[i];
  cout << (n == 1 ? 0 : calc(n)) << '\n';
  return 0;
}

P2398 GCD SUM

Problem

给定正整数 \(n\)，求 \(\sum\limits_{i=1}^n \sum\limits_{j=1}^n \gcd(i, j)\).

Solve

欧拉反演模板题。

原式：

\[\sum\limits_{i=1}^n \sum\limits_{j=1}^n (i, j) \]

嵌入式欧拉反演：

\[\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sum\limits_{d|(i,j)}\varphi(d) \]

由 \(d|(i,j)\Rightarrow d|i,d|j\) 可知：

\[\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sum\limits_{d|i,d|j}\varphi(d) \]

转化形式：

\[\sum\limits_{d=1}^n\varphi(d)\sum\limits_{d|i}^n\sum\limits_{d|j}^n 1 \]

化简：

\[\sum\limits_{d=1}^n\varphi(d)\big(\lfloor\frac{n}{d}\rfloor\big)^2 \]

前缀和 + 整除分块即可求解。注意欧拉筛欧拉函数时，是以质数为讨论根据。

Code

点击查看代码

#include <bits/stdc++.h>
#define int long long
#define H 19260817
#define rint register int
#define For(i,l,r) for(rint i=l;i<=r;++i)
#define FOR(i,r,l) for(rint i=r;i>=l;--i)
#define MOD 1000003
#define mod 1000000007

using namespace std;

inline int read() {
  rint x=0,f=1;char ch=getchar();
  while(ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
  while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
  return x*f;
}

void print(int x){
  if(x<0){putchar('-');x=-x;}
  if(x>9){print(x/10);putchar(x%10+'0');}
  else putchar(x+'0');
  return;
}

const int N = 1e5 + 10;

int n, prime[N], tot, phi[N], sum[N];

bool vis[N];

void Eular() {
  vis[1] = phi[1] = 1;
  for (int i = 1; i <= N - 1; i++) {
    if(!vis[i]) prime[++tot] = i, phi[i] = i-1;
    for (int j = 1; j <= tot && i * prime[j] <= N - 1; j++) {
      vis[i * prime[j]] = 1;
      phi[i * prime[j]] = (i % prime[j] == 0 ? phi[i] * prime[j] : phi[i] * (prime[j]-1));
      if(i % prime[j] == 0) break;
    }
  }
}

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

signed main() {
  n = read();
  Eular();
  For(i,1,n) sum[i] = sum[i-1] + phi[i];
  cout << calc(n) << '\n';
  return 0;
}