杜教筛瞎推[学习笔记]
O, 前言
对于 bzoj3944 来说, 和莫比乌斯反演等其他知识关系不大, 但是 \(\mu\) 函数在自变量较大情况下的前缀和在反演题中也是会被用到的.
接下来通过 bzoj3944 Sum 一题引入杜教筛的思想.
参考资料
铃悬的数学小讲堂 -- 杜教筛 by 铃悬 https://lx-2003.blog.luogu.org/dujiao-sieve
一, 例题
给定一个正整数 \(N(0\le N\le 2^{31}-1)\), 求 \(\sum_{i=1}^n\varphi(i)\) 和 \(\sum_{i=1}^n\mu(i)\). 多组询问.
1. 分析
对于多组询问求积性函数前缀和, 我们可以 \(O(n)\) 线筛预处理,\(O(1)\) 询问. 其中空间复杂度当然也是 \(O(n)\) 的, 所以我们无法直接线筛, 但是之后的部分还是会用到的.
2. 推导
其实地上本没有路, 走的人多了, 也便成了路.-- 鲁迅
对于一般的问题, 我们要求 \(S(n)=\sum_{i=1}^nf(i)?\).
构造 \(g(n)\),\((f*g)(n)=\sum_{d|n}f(d)g\left(\frac nd\right)\), 我们要求的是 \(S(n)\), 式子中可以从右边分离出 \(f(n)\), 再进一步推导, 求出方程左边的前缀和, 有
\[ \begin{aligned} \sum_{i=1}^n(f*g)(i)&=\sum_{i=1}^n\sum_{d|i}f(d)g\left(\frac id\right)\&=\sum_{i=1}^n\sum_{x=1}^i\sum_{xy=i}f(x)g(y)\\end{aligned} \]
因为每对数 \((x,y)\) 最多只会做一次贡献, 所以我们可以直接枚举 \(x\cdot y\) 的结果来统计贡献.
\[ \begin{aligned} \sum_{i=1}^n(f*g)(i)&=\sum_{y=1}^n\sum_{xy\le n}f(x)g(y)\&=\sum_{y=1}^ng(y)\sum_{xy\le n}f(x)\&=\sum_{y=1}^ng(y)\sum_{i=1}^{\left\lfloor\frac ny\right\rfloor}f(i) \end{aligned} \]
我们在上面定义了 \(S(n)\), 所以式子继续化成了
\[ \sum_{i=1}^n(f*g)(i)=\sum_{y=1}^ng(y)S\left(\left\lfloor\frac ny \right\rfloor\right) \]
这时我们可以整出 \(S(n)\) 来了, 当右侧式子中 \(y=1\) 时, 变为 \(g(1)S(n)\).
移项得到
\[ g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{y=2}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right) \]
这时我们看到了可以数论分块的 \(\left\lfloor\frac ny\right\rfloor\), 看上去像个递归过程. 假定 \(g(n)\) 的前缀和可以在 \(O(1)\) 的复杂度内求出, 且 \((f*g)(i)\) 的前缀和 -- 即 \(\sum_{i=1}^n(f*g)(i)\) 可以在 \(O(1)\sim O(\sqrt n)\) 的时间内求出, 则可以保证复杂度.
因此构造 \(g?\) 是一个重点.
3. 复杂度
在构造之前, 先考虑一下复杂度.
假定在理想状态下, 求出 \(g(n)\) 的前缀和和 \((f*g)(i)\) 的前缀和都是 \(O(1)\) 的.
复杂度分析为
\[ T(n)=\sum_{i=1}^\sqrt n \sqrt{\left\lfloor\frac ni\right\rfloor}=\int_1^\sqrt n\sqrt{\frac nx}\mathrm dx \]
式子可以化为 \(\sqrt n\times x^{-\frac 12} \mathrm dx\), 找到函数 \(\sqrt n\times 2x^{\frac 12}=2\sqrt {nx}\)? 的导数是它, 然后带入
\[ \begin{aligned} T(n)&=\int_1^\sqrt n\sqrt n\times x^{-\frac 12}\mathrm dx\&=2\sqrt{n\sqrt n}-2\sqrt n\&=O\left(n^{\frac 34}\right) \end{aligned} \]
同时, 我们可以线筛预处理出前面一段的答案. 此外, 我们从 \(S(n)\) 进入整个递归过程, 这之后也只会调用自变量为 \(\left\lfloor\frac ny\right\rfloor\) 的 \(S\) 函数. 接着, 由于 \(\left\lfloor\frac{\left\lfloor\frac ny\right\rfloor}{x}\right\rfloor=\left\lfloor\frac{n}{xy}\right\rfloor\), 自变量只可能是 \(\left\lfloor\frac ni\right\rfloor\) 的 \(O(\sqrt n)?\) 种取值中的一种, 因此可以记忆化.
假设我们预处理出了前 \(n^c\) 的答案, 那么后面一部分积分的上界就是 \(n^{1-c}\) 了, 因为只有当 \(k<c\) 时需要计算 \(S\left(\frac{n}{n^{k}}\right)=S\left(n^{1-k}\right)\).
则
\[ \begin{aligned} T(n)&=O(n^c)+\int_1^{n^{1-c}}\sqrt n\times x^{-\frac 12}\mathrm dx\&=O(n^c)+2\sqrt{n\cdot n^{1-c}}-2\sqrt n\&=O(n^c)+O\left(n^{1-\frac c2}\right) \end{aligned} \]
由基本不等式得二者取等时 \(T(n)\) 有最小值, 此时 \(c=1-\frac c2,\ c=\frac 23\),\(T(n)_\min=O\left(n^\frac 23\right)\).
4. 构造
对于一些常见的函数如 \(id(n),1(n),\epsilon(n)\),\(g(n)\) 的前缀和比较好求, 因此在杜教筛时优先考虑这些函数.
对于 \(\varphi(n)\), 由于 \((\varphi*1)(n)=n\),\((\varphi*1)(n)\) 和 \(1(n)\) 的前缀和都比较好求, 所以可以利用它.
把它带入移项后的式子:
\[ \begin{aligned} g(1)S(n)&=\sum_{i=1}^n(f*g)(i)-\sum_{y=2}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right)\1\times S(n)&=\sum_{i=1}^n(\varphi*1)(i)-\sum_{y=2}^n1\times S\left(\left\lfloor\frac ny\right\rfloor\right)\S(n)&=\frac{n(n+1)}2-\sum_{y=2}^nS\left(\left\lfloor\frac ny\right\rfloor\right) \end{aligned} \]
然后递归就可以了.
对于 \(S(x)\) 的计算(\(n\) 是常数)
当 \(x\le \left(2^{31}-1\right)^{\frac 23}\) 时, 可以直接返回已经存下的值, 其余情况根据关于 \(n\) 的分母, 即 \(\left\lfloor\frac nx\right\rfloor\) 的值查询是否记忆过这个答案, 此时只需要存 \(O\left(n^\frac 13\right)\) 种状态.
5. 代码
- #include<cstdio>
- #include<cstring>
- #define ll long long
- const int mxm=2000000;
- int pri[1<<20],cnt=0,n;
- ll phi[2000100],mu[2000100];
- bool is[2000100];
- ll f[2000100],g[2000100];
- bool used[2000];
- ll calc(int x)
- {
- if(x<=mxm)
- return phi[x];
- int t=n/x;
- if(used[t])
- return g[t];
- used[t]=1;
- ll ans=(ll)(1ll+x)*x/2;
- for(int i=2;i<=x;++i)
- {
- int nxt=x/(x/i);
- ans-=calc(x/i)*(nxt-i+1);
- i=nxt;
- if(i==2147483647)
- break;
- }
- return g[t]=ans;
- }
- ll calc1(int x)
- {
- if(x<=mxm)
- return mu[x];
- int t=n/x;
- if(used[t])
- return g[t];
- used[t]=1;
- ll ans=1;
- for(int i=2;i<=x;++i)
- {
- int nxt=x/(x/i);
- ans-=calc1(x/i)*(nxt-i+1);
- i=nxt;
- if(i==2147483647)
- break;
- }
- return g[t]=ans;
- }
- int main()
- {
- is[0]=is[1]=1;
- phi[1]=1;
- mu[1]=1;
- for(int i=2;i<=mxm;++i)
- {
- if(!is[i])
- {
- pri[++cnt]=i;
- mu[i]=-1;
- phi[i]=i-1;
- }
- for(int j=1;j<=cnt&&i*pri[j]<=mxm;++j)
- {
- is[i*pri[j]]=1;
- if(i%pri[j]==0)
- {
- mu[i*pri[j]]=0;
- phi[i*pri[j]]=pri[j]*phi[i];
- break;
- }
- else
- {
- mu[i*pri[j]]=-mu[i];
- phi[i*pri[j]]=(pri[j]-1)*phi[i];
- }
- }
- mu[i]+=mu[i-1];
- phi[i]+=phi[i-1];
- }
- int T;
- scanf("%d",&T);
- while(T--)
- {
- scanf("%d",&n);
- memset(used,0,sizeof(used));
- printf("%lld",calc(n));
- memset(used,0,sizeof(used));
- printf("%lld\n",calc1(n));
- }
- return 0;
- }
在数论题中最好通过带入最大的数 (在本题中是 \(2147483647\)) 来检验程序的运行时间和数据规模. 在本题中, 当循环到 \(i=2147483647\) 时,++i 就会爆 int, 在这里需要特判.
数组不要开小了!!
二, 练习
感觉把反演给忘了, 先复习反演.
??????
杜教筛瞎推 学习笔记[杜教筛] [数学]
来源: http://www.bubuko.com/infodetail-3019028.html