[演算法入門]杜教篩
#1.0 基礎知識
#1.1 積性函式
#1.1.1 定義
設 \(f\) 是數論函式,若對任意互質的正整數 \(a,b\),都有 \(f(ab)=f(a)f(b)\),則稱 \(f\) 是積性函式。若對任意的正整數 \(a,b\),都有 \(f(ab)=f(a)f(b)\),則稱 \(f\) 是完全積性的。
#1.1.2 性質
若 \(f\) 是積性函式,且 \(n=p_1^{\alpha_1}p_2^{\alpha_2}···p_s^{\alpha_s}\) 是 \(n\) 的標準分解,則有 \(f(n)=f(p_1^{\alpha_1})f(p_2^{\alpha_2})···f(p_s^{\alpha_s})\)
#1.1.3 一些常見的積性函式
- 單位函式
單位函式 \(\epsilon(n)\) 定義為:
\[\epsilon(n)=[n=1] =\begin{cases}1,\quad n=1,\\0,\quad n\ne1.\end{cases} \]其中 \([\text{condition}]\) 表示當 \(\text{condition}\) 為真時取值為 \(1\),否則為 \(0\) 的函式。單位函式是完全積性函式。
- 除數函式
除數函式 \(\sigma_k(n)\)
約數個數 \(\sigma_0(n)\) 常記為 \(d(n)\),約數和 \(\sigma_1(n)\) 常記為 \(\sigma(n).\)
除數函式都是積性函式。
- \(\text{Euler}\) 函式
\(\text{Euler}\) 函式 \(\varphi(n)\) 表示不超過 \(n\) 且與 \(n\) 互質的正整數的個數。由 \(n\) 的標準分解並結合容斥原理,我們可以得到 \(\text{Euler}\) 函式的顯示錶達式:
\[\varphi(n)=n\cdot\prod\limits_{i=1}^s(1−\dfrac{1}{p_i}), \]其中 \(n=p_1^{\alpha_1}p_2^{\alpha_2}···p_s^{\alpha_s}\)
對於任意 \(n\),\(\text{Euler}\) 函式有如下性質:
\[n=\sum\limits_{d|n}\varphi(d). \]證明:
若 \(\gcd(n,i)=d\),那麼 \(\gcd(\dfrac{n}{d},\dfrac{i}{d})=1\)。而 \(\dfrac{i}{d}\leq\dfrac{n}{d},\dfrac{i}{d}\in\Z\),故這樣的 \(i\) 有 \(φ(\dfrac{n}{d})\) 個。考慮所有 \(d|n\),我們也就考慮到了所有 \(1\) 到 \(n\) 之間的 \(n\) 個整數,因此有
\[n=\sum\limits_{d|n}\varphi(\dfrac{n}{d})=\sum\limits_{d|n}\varphi(d). \]證畢.-
冪函式 \(id_k(n)\), \(id_k(n)=n^k,id=id_1.\)
-
恆等函式 \(I(n)\),這個函式的值恆為 \(1.\)
-
\(\text{Mobius}\) 函式
其中 \(p_1,\cdots,p_s\) 是不同素數。可以看出,\(\mu(n)\) 恰在 \(n\) 無平方因子時非零。易見 \(\mu\) 是積性函式。
\(\text{Mobius}\) 函式具有以下性質:
\[\sum\limits_{d|n}\mu(d)=\epsilon(n). \]證明:
當 \(n = 1\) 時顯然成立,下面證 \(n>1\) 時結論成立。
當根據 \(\text{Mobius}\) 函式的定義,當且僅當 \(n\) 無平方因子時非零,故能對答案造成貢獻的 \(d\) 中每個素因子的次數只能為 \(0\) 或 \(1\),設 \(n\) 的唯一分解式為 \(n=p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_s^{\alpha_s}\),於是有:
\[\sum\limits_{d|n}\mu(d)=\sum\limits_{k=0}^s(-1)^k\dbinom{s}{k}=(1-1)^s=0. \]證畢.#1.2 \(\text{Dirichlet}\) 卷積
#1.2.1 定義
設 \(f,g\) 是數論函式,考慮數論函式 \(h\) 滿足
\[h(n)=\sum\limits_{d|n}f(d)g(\dfrac{n}{d}), \]則稱 \(h\) 為 \(f\) 和 \(g\) 的 \(\text{Dirichlet}\) 卷積,記作 \(h=f*g.\)
#1.2.2 性質
單位函式 \(\epsilon\) 是 \(\text{Dirichlet}\) 卷積的單位元,即對於任意函式 \(f\),有 \(\epsilon∗f=f∗\epsilon=f\)。\(\text{Dirichlet}\) 卷積滿足交換律和結合律。
如果 \(f,g\) 都是積性函式,那麼 \(f∗g\) 也是積性函式。
#1.2.3 一些關係
通過觀察,不難發現有些積性函式可以用 \(\text{Dirichlet}\) 卷積的形式聯絡起來,比如:
- \(\mu*I=\epsilon.\)
- \(\varphi*I=id.\)
以上兩個是較為常用的轉換關係。
#1.3 整除分塊
定理:\(\left\lfloor\tfrac{n}{d}\right\rfloor\) 的值不超過 \(2\sqrt n\) 個。很顯然,不證。
假設我們要求 \(\sum_{i=1}^n\left\lfloor\frac{n}{i}\right\rfloor.\)
通過打表(以 \(n=10\) 為例),我們得到:
不難發現,對於每一個值相同的塊,它的起始點為 \(i\),它的最後一個數就是 \(n/(n/i).\)
for (int l = 1,r;l <= n;l = r + 1){
r = n / (n / l);
ans += (r - l + 1) * (n / l);
}
再根據我們一開始得到的定理可知,上面程式的時間複雜度為 \(O(\sqrt n).\)
#2.0 杜教篩(SUM)
#2.1 有啥用?
杜教篩是一種利用 \(\text{Dirichlet}\) 卷積來構造遞推式,從而對積性函式進行快速求和的方法。
比如,我要求 \(\sum\limits_{i=1}^n\varphi(i)\),當 \(n\) 的規模在 \(10^7\) 時,線性篩還可以勉強一戰,再高一些,就需要更優的方法了。至於杜教篩為啥更優,後面再說。
#2.2 怎麼辦?
假定我們要求 \(S(n)=\sum\limits_{i=1}^nf(i)\),我們需要構造出這樣兩個積性函式 \(h,g\),滿足 \(h=f*g\)。
\[\begin{aligned} \sum\limits_{i=1}^nh(i)&=\sum\limits_{i=1}^n\sum\limits_{d|i}g(d)\cdot f(\dfrac{i}{d})\\ &=\sum\limits_{d=1}^ng(d)\cdot\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f(i)\\ &=\sum\limits_{d=1}^ng(d)\cdot S(\left\lfloor\dfrac{n}{d}\right\rfloor)\\ &=g(1)\cdot S(n)+\sum\limits_{d=2}^ng(d)\cdot S(\left\lfloor\dfrac{n}{d}\right\rfloor)\\ \end{aligned} \]於是有
\[g(1)\cdot S(n)=\sum\limits_{i=1}^nh(i)-\sum\limits_{d=2}^ng(d)\cdot S(\left\lfloor\dfrac{n}{d}\right\rfloor).\tag1 \]式 \((1)\) 便是杜教篩的一個常用形式。顯然我們希望構造出的 \(h\) 的字首和要好求一些。
拋開右式前半部分不談,我們來看後半部分 \(\sum_{d=2}^xg(d)\cdot S(\left\lfloor\tfrac{n}{d}\right\rfloor).\)
看到這個熟悉的形式,不由得讓我們想起了 整除分塊,然後我們可以遞迴地解決,並採用記憶化搜尋。
當然,我們可以與線性篩結合,進一步降低時間複雜度。
#2.3 時間複雜度分析
演算法的總時間複雜度就是計算所有特殊點處的函式值的時間複雜度。
回憶特殊點的結構,時間複雜度 \(T(n)\) 可以估計為
\[T(n)=\sum\limits_{i=1}^\sqrt nO\left(\sqrt i\right)+\sum\limits_{i=1}^\sqrt nO\left(\sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor}\right). \]顯然式中第一項漸進意義上小於第二項。
而對於式中第二項我們可以利用積分估計:
\[\sum\limits_{i=1}^\sqrt nO\left(\sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor}\right)=O\left(\int_1^{\sqrt n}\sqrt{\dfrac{n}{x}}dx\right)=O(n^{\frac{1}{2}}\cdot n^{\frac{1}{4}})=O(n^{\frac{3}{4}}). \]於是演算法的時間複雜度為 \(O(n^{\frac{3}{4}})\)。
假設我們使用 \(\text{Euler}\) 篩預先求出了 \(\varphi\) 的前 \(S\) 項,那麼遞迴部分的時間複雜度變為:
\[\sum\limits_{i=1}^{\frac{n}{S}}O\left(\sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor}\right)=O\left(\int_1^{\frac{n}{S}}\sqrt{\dfrac{n}{x}}dx\right)=O\left(n^{\frac{1}{2}}\cdot\sqrt{\dfrac{n}{S}}\right)=O\left(\dfrac{n}{\sqrt{S}}\right). \]結合 \(\text{Euler}\) 篩的時間複雜度為 \(O(S)\),於是總體時間複雜度為 \(O\left(S+\frac{n}{\sqrt{S}}\right).\)
如果取 \(S=n^{\frac{2}{3}}\),於是總體時間複雜度為 \(O(n^{\frac{2}{3}})\).
#2.4 例題
練習是必不可少的。
#2.4.1 \(\text{Euler}\) 函式
求 \(\sum_{i=1}^n\varphi(i).\)
熟知
\[\varphi*I=id. \]令 \(\phi(n)=\sum_{i=1}^n\varphi(i)\),於是有
\[\begin{aligned} \phi(n)=\sum\limits_{i=1}^ni-\sum\limits_{d=2}^n\phi(\left\lfloor\dfrac{n}{d}\right\rfloor). \end{aligned} \]#2.4.2 \(\text{Mobius}\) 函式
求 \(\sum_{i=1}^n\mu(i).\)
同樣考慮
\[\mu*I=\epsilon. \]令 \(\Mu=\sum_{i=1}^n\mu(i)\),於是有
\[\begin{aligned} \Mu(n)=1-\sum\limits_{d=2}^n\Mu(\left\lfloor\dfrac{n}{d}\right\rfloor). \end{aligned} \]#2.5 程式碼實現
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <tr1/unordered_map> //unordered_map 在 c++98 中需要該標頭檔案
#define ll long long
#define mset(l,x) memset(l,x,sizeof(l))
using namespace std;
using namespace tr1; //需要該名稱空間
const int N = 6000010;
const int INF = 0x3fffffff;
const int L = 5000000;
bool vis[N];
int mu[N],phi[N],sum1[N],n,prim[N],cnt,t;
ll sum2[N];
unordered_map <int,int> w1;
unordered_map <ll,ll> w2;
inline void prework(int x){
phi[1] = mu[1] = 1;
for (int i = 2;i <= x;i ++){
if (!vis[i]){
prim[++ cnt] = i;
phi[i] = i - 1;mu[i] = -1;
}
for (int j = 1;j <= cnt && prim[j] * i <= x;j ++){
vis[prim[j] * i] = 1;
if (i % prim[j] == 0){
phi[prim[j] * i] = phi[i] * prim[j];
break;
}
else {
mu[i * prim[j]] = -mu[i];
phi[i * prim[j]] = phi[i] * (prim[j] - 1);
}
}
}
for (int i = 1;i <= x;i ++)
sum1[i] = sum1[i - 1] + mu[i],sum2[i] = sum2[i - 1] + phi[i];
}
inline int djsmu(int x){
if (x <= L) return sum1[x];
if (w1[x]) return w1[x];
int ans = 1;
for (ll l = 2,r;l <= x;l = r + 1){
r = x / (x / l);
ans -= (r - l + 1) * djsmu(x / l);
}
return w1[x] = ans;
}
inline ll djsphi(ll x){
if (x <= L) return sum2[x];
if (w2[x]) return w2[x];
ll ans = x * (1 + x) / 2;
for (ll l = 2,r;l <= x;l = r + 1){
r = x / (x / l);
ans -= (r - l + 1) * djsphi(x / l);
}
return w2[x] = ans;
}
int main(){
scanf("%d",&t);
prework(L);
while (t --){
scanf("%d",&n);
printf("%lld %d\n",djsphi(n),djsmu(n));
}
return 0;
}
上面使用的 unordered_map
本身是衝著比 map
少個 \(\log\) 才用的,但發現實際用時差的不多(而且似乎還容易被卡)。