【算法讲12:杜教筛入门】
前置知识
- 积性函数 与 狄利克雷卷积 【算法讲7:积性函数(下)】
- 数论分块 【算法讲6:数论分块(整除分块)】
- 莫比乌斯反演 与 欧拉筛 【算法讲8:莫比乌斯函数及其反演(理论部分) | 欧拉筛】
- 记忆化搜索。应该学过搜索的人都会的吧…
引入
- 【问题描述】【模板】杜教筛 | 洛谷 P4213
给定 n n n ,求
∑ i = 1 n φ ( i ) 与 ∑ i = 1 n μ ( i ) \sum_{i=1}^n\varphi(i)\qquad 与\qquad \sum_{i=1}^n\mu(i) i=1∑nφ(i)与i=1∑nμ(i) - 【数据范围】
样例组数 T ≤ 10 T\le 10 T≤10
1 ≤ n ≤ 2 31 1\le n\le 2^{31} 1≤n≤231
思路
- 【暴力的做法】如果
n
n
n 是
1
e
7
1e7
1e7 级别的,我们可以
O
(
n
)
O(n)
O(n) 直接用 欧拉筛 筛出每一个函数值,然后用前缀求和即可。
但是这里 n n n 达到了 2 e 9 2e9 2e9 ,明显无论是 时间 还是 空间 都无法完成。于是我们有了 杜教筛 。 - 首先,设
f
f
f 是我们需要求的前缀和的积性函数。我们需要找到函数
g
、
h
g、h
g、h,满足:
h = f ∗ g h=f*g h=f∗g
我们设
S ( n ) = ∑ i = 1 n f ( i ) S(n)=\sum_{i=1}^nf(i) S(n)=i=1∑nf(i)
表示我们需要求的前缀和。接下来,我们推一会儿狮子:
∑ i = 1 n h ( i ) = ∑ i = 1 n ∑ d ∣ i g ( d ) f ( i d ) = ∑ d = 1 n g ( d ) ∑ i = 1 ⌊ n d ⌋ f ( i ) = ∑ d = 1 n g ( d ) S ( ⌊ n d ⌋ ) = g ( 1 ) S ( n ) + ∑ d = 2 n g ( d ) S ( ⌊ n d ⌋ ) \begin{aligned} \sum_{i=1}^n h(i) &= \sum_{i=1}^n \sum_{d|i}g(d)f(\frac{i}{d})\\ &= \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} f(i)\\ &= \sum_{d=1}^ng(d)S(\lfloor\frac{n}{d}\rfloor)\\ &=g(1)S(n)+\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)\\ \end{aligned} i=1∑nh(i)=i=1∑nd∣i∑g(d)f(di)=d=1∑ng(d)i=1∑⌊dn⌋f(i)=d=1∑ng(d)S(⌊dn⌋)=g(1)S(n)+d=2∑ng(d)S(⌊dn⌋)
于是,我们得到了杜教筛的式子:
g ( 1 ) S ( n ) = ∑ i = 1 n h ( i ) − ∑ d = 2 n g ( d ) S ( ⌊ n d ⌋ ) g(1)S(n)=\sum_{i=1}^n h(i)-\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)\\ g(1)S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋) - 那么杜教筛有什么注意事项?
我们需要寻找合适的 g 、 h g、h g、h,以至于我们可以很方便求出 ∑ n i = 1 h ( i ) \underset{i=1}{\overset{n}{\sum}}h(i) i=1∑nh(i)
后面那一项,可以用 整除分块 快速算出。
如果我们直接记忆化计算 S ( n ) S(n) S(n) ,时间复杂度大概为 O ( n 3 4 ) O(n^{\frac{3}{4}}) O(n43)
如果我们预先算出前 n \sqrt n n 项 f ( i ) f(i) f(i),那么时间复杂度大概为 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)
对于 φ \varphi φ 的杜教筛
- 关键是怎么找合适的
h
=
f
∗
g
h=f*g
h=f∗g
我们之前有 φ ∗ 1 = i d \varphi*1=id φ∗1=id,这是一个很好的卷积公式。直接带入即可得到:
S ( n ) = ∑ i = 1 n i − ∑ d = 2 n S ( ⌊ n d ⌋ ) S(n)=\sum_{i=1}^n i-\sum_{d=2}^nS(\lfloor\frac{n}{d}\rfloor) S(n)=i=1∑ni−d=2∑nS(⌊dn⌋)
h ( n ) h(n) h(n) 的前缀和非常好算,就是 n ( n + 1 ) 2 \cfrac{n(n+1)}{2} 2n(n+1),于是完成了。
对于 μ \mu μ 的杜教筛
- 关键是怎么找合适的
h
=
f
∗
g
h=f*g
h=f∗g
我们之前有 μ ∗ 1 = ε \mu*1=\varepsilon μ∗1=ε,这是一个很好的卷积公式。直接带入即可得到:
S ( n ) = ∑ i = 1 n ε ( i ) − ∑ d = 2 n S ( ⌊ n d ⌋ ) S(n)=\sum_{i=1}^n \varepsilon(i)-\sum_{d=2}^nS(\lfloor\frac{n}{d}\rfloor) S(n)=i=1∑nε(i)−d=2∑nS(⌊dn⌋)
h ( n ) h(n) h(n) 的前缀和非常好算,就是 1 1 1,于是又完成了。
核心代码
- 时间复杂度: O ( n + n 2 3 ) O(\sqrt n+n^{\frac{2}{3}}) O(n+n32)
/*
_ __ __ _ _
| | \ \ / / | | (_)
| |__ _ _ \ V /__ _ _ __ | | ___ _
| '_ \| | | | \ // _` | '_ \| | / _ \ |
| |_) | |_| | | | (_| | | | | |___| __/ |
|_.__/ \__, | \_/\__,_|_| |_\_____/\___|_|
__/ |
|___/
*/
const int MAX = 1e7+50;
const int TMAX = MAX - 50;
int cnt;
ll phi[MAX];
int mu[MAX];
int vis[MAX],prime[MAX];
void shai(int n){
phi[1] = mu[1] = 1;
for(int i = 2;i <= n;++i){
if(!vis[i]){
prime[++cnt] = i;
phi[i] = i-1;
mu[i] = -1;
}
for(int j = 1;j <= cnt && i * prime[j] <= n;++j){
vis[i * prime[j]] = 1;
if(i % prime[j]){
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
mu[i * prime[j]] = -mu[i];
}else{
phi[i * prime[j]] = phi[i] * prime[j];
mu[i * prime[j]] = 0;
break;
}
}
}
for(int i = 1;i <= n;++i){
phi[i] += phi[i-1];
mu[i] += mu[i-1];
}
}
unordered_map<int,ll>PHI;
ll fd_phi(int n){
if(n <= TMAX)return phi[n];
if(PHI[n])return PHI[n];
ll L = 2,R = 0;
ll res = 0;
while(L <= n){
R = n / (n / L);
res += (R - L + 1) * fd_phi(n / L);
L = R + 1;
}
res = (1LL+n)*n/2 - res;
PHI[n] = res;
return res;
}
unordered_map<int,int>MU;
ll fd_mu(int n){
if(n <= TMAX)return mu[n];
if(MU[n])return MU[n];
ll L = 2,R = 0;
ll res = 0;
while(L <= n){
R = n / (n / L);
res += (R - L + 1) * fd_mu(n / L);
L = R + 1;
}
res = 1 - res;
MU[n] = res;
return res;
}
int main()
{
shai(TMAX);
int T;scanf("%d",&T);
while(T--){
int n;scanf("%d",&n);
printf("%lld %d\n",fd_phi(n),fd_mu(n));
}
return 0;
}
例子
来看一个比较综合的例子。
- 【题目描述】function | HDU5608
给定一个数论函数 f f f,它满足以下性质:
N 2 − 3 N + 2 = ∑ d ∣ N f ( d ) N^2-3N+2=\sum_{d|N}f(d) N2−3N+2=d∣N∑f(d)
你需要求出:
∑ i = 1 n f ( i ) \sum_{i=1}^nf(i) i=1∑nf(i) - 【数据范围】
限时 3000 M s 3000Ms 3000Ms
内存 65536 K b 65536Kb 65536Kb
样例组数 T ≤ 500 T\le 500 T≤500
N ≤ 1 0 9 N\le 10^9 N≤109 - 【思路】
我们设 F ( n ) = ∑ d ∣ n f ( i ) F(n)=\underset{d|n}{\sum}f(i) F(n)=d∣n∑f(i)。根据莫比乌斯反演,我们有: f ( i ) = ∑ d ∣ n μ ( d ) F ( n d ) f(i)=\underset{d|n}{\sum}\mu(d)F(\frac{n}{d}) f(i)=d∣n∑μ(d)F(dn)
我们容易在 O ( n log n ) O(n\log n) O(nlogn) 的时间内算出前 1 e 6 1e6 1e6 个 f ( i ) f(i) f(i)
接下来,我们就要使用杜教筛了。因为题目已知 f ∗ 1 = F f*1=F f∗1=F,我们使用这个即可:
S ( n ) = ∑ i = 1 n F ( i ) − ∑ i = 2 n S ( ⌊ n i ⌋ ) S(n)=\sum_{i=1}^nF(i)-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor) S(n)=i=1∑nF(i)−i=2∑nS(⌊in⌋)
前面的我们很好得到,因为 F ( n ) = n 2 − 3 n + 2 F(n)=n^2-3n+2 F(n)=n2−3n+2。后面的数论分块即可。 - 【一些优化】
(1)你可能算了半天某个值 f ( i ) = 0 f(i)=0 f(i)=0,但是你如果直接调用 i f ( f ( i ) ! = 0 ) if(f(i)!=0) if(f(i)!=0) 来判断这个值有没有算过不可取。应该用 ( i t = f i n d ( i ) ) ! = F . e n d ( ) (it=find(i))!=F.end() (it=find(i))!=F.end()来证明之前确实没有储存过该值。
(2)因为用到 f i n d find find ,所以我们选择 m a p map map 更加稳妥,因为哈希表的储存和查找,在数比较密集的情况下还是可能会 O ( N ) O(N) O(N) 的。
(3)因为哈希表是像数组一样 O ( N ) O(N) O(N) 内存的,但是 m a p map map 内部为红黑树,如果你 1 ∼ 1 e 6 1\sim 1e6 1∼1e6 的值全部存在 m a p map map 中,之后很多新加进来的数字再一存,会导致 M L E MLE MLE ,所以小值存数组里还是必要的。
(4)取模操作,一些加减法后取模操作尽量变成 i f ( x ≥ M O D ) x − = M O D if(x\ge MOD)x-=MOD if(x≥MOD)x−=MOD 之类的语句
(5)可以化简的式子直接化简,不然乘法多了取模多了又容易 T L E TLE TLE
这几条优化,把 T L E 、 M L E TLE、MLE TLE、MLE 的代码直接纠正了回来。
本来我没有用以上几个优化,预处理只能处理了 5 e 5 5e5 5e5 个数字,勉强 1800 M s 1800Ms 1800Ms 和 24600 K 24600K 24600K。
但是看到他们的提交,时间不会卡的这么死,就仔细研究了一下可以应用到的优化。
核心代码
- 时间:
780
M
s
/
3000
M
s
780Ms/3000Ms
780Ms/3000Ms
内存: 14612 K / 65536 K 14612K/65536K 14612K/65536K
时间复杂度: O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)
/*
_ __ __ _ _
| | \ \ / / | | (_)
| |__ _ _ \ V /__ _ _ __ | | ___ _
| '_ \| | | | \ // _` | '_ \| | / _ \ |
| |_) | |_| | | | (_| | | | | |___| __/ |
|_.__/ \__, | \_/\__,_|_| |_\_____/\___|_|
__/ |
|___/
*/
const int MAX = 1e6+50;
const ll MOD = 1e9+7;
ll qpow(ll a,ll n){/* */ll res = 1LL;while(n){if(n&1)res=res*a%MOD;a=a*a%MOD;n>>=1;}return res;}
ll inv(ll a){/* */return qpow(a,MOD-2);}
const int TMAX = MAX - 50;
bool vis[MAX];
int prime[MAX];
int mu[MAX];
int cnt;
ll f[MAX]; /// 存小值
map<int,ll>F; /// 存大值
map<int,ll>::iterator it;
const ll iv3 = inv(3);
ll qiu(ll x){
return x * x - 3 * x + 2;
}
void shai(int n){
mu[1] = 1;
for(int i = 2;i <= n;++i){
if(!vis[i]){
prime[++cnt] = i;
mu[i] = -1;
}
for(int j = 1;j <= cnt && i * prime[j] <= n;++j){
vis[i*prime[j]] = 1;
if(i % prime[j])mu[i*prime[j]] = -mu[i];
else {
mu[i*prime[j]] = 0;
break;
}
}
}
for(int i = 1; i <= n; i++){
for(int j = i; j <= n; j += i){
f[j] = (f[j] + 1LL * mu[i] * qiu(j / i)) % MOD;
}
f[i] += f[i - 1];
f[i] = (f[i] % MOD + MOD) % MOD;
}
}
ll fd_func(ll n){
if(n <= TMAX)return f[n];
if((it = F.find(n)) != F.end())return it->second;
ll L = 2,R = 0;
ll res = n * (n + 1) % MOD * (n - 4) % MOD * iv3 % MOD + 2 * n;
res = (res % MOD + MOD) % MOD;
while(L <= n){
R = n / (n / L);
res -= (R - L + 1) * fd_func(n / L) % MOD;
if(res < 0)res += MOD;
L = R + 1;
}
return F[n] = res;
}
int main()
{
shai(TMAX);
int T;scanf("%d",&T);
while(T--){
ll a;scanf("%lld",&a);
printf("%lld\n",fd_func(a));
}
return 0;
}