P1009 加强版——原来这个东西可以分治啊？？？
gitxiaozheng
·
·
算法·理论
前言
作者在发布一篇介绍 Binary Splitting 的文章后，有人这么吐槽：
确实麻烦，那能不能把 Binary Splitting 扩展运用到只与高精度整数有关的级数运算呢？
引入
::::info[P1009 加强版]{open}
原题见 U660888。
给定正整数 n，求 S=\sum^{n}_{i=1}i! 的值。数据范围为 1\le n \le 1.5\times10^5。
时间限制：2500ms。
::::
大家都知道，加强之前仅仅是个简单的高精度基础题。而现在，别看这题时间限制好像很宽松，其实直接用压 9 位的高精度乘法实现暴力算法连 n\le 5\times 10^4 的数据都过不去。
于是我们尝试分治。分治计算单个阶乘再相加是好想的，但仍然是不够快，甚至劣于暴力算法。其实，我们可以用一个分治算法来解决它。
前置知识
当然，你可能需要通过前言里的文章认识一下 Binary Splitting，这会有助于你理解算法的推导过程。当然你愣记结论也行。
算法介绍
设 n 为正整数，有函数 P(x) 和 R(x)，则如果我们遇到形如下列形式的式子：
\sum^{n}_{i=1} {P(i)}\prod^{i-1}_{j=1}{R(j)}
我们就可以分治计算。
定义部分和 S(l,r) 为：
S(l,r)=\sum^{r}_{i=l} {P(i)}\prod^{i-1}_{j=l}{R(j)}
然后定义辅助函数 R(l,r) 为：
R(l,r)=\prod^{r}_{i=l}R(i)
分治时取 m=\lfloor \frac{l+r}{2}\rfloor，然后使用下列分治递推式：
\begin{cases}
S(l,r)=S(l,m)+R(l,m) \cdot S(m+1,r)\\
R(l,r)=R(l,m) \cdot R(m+1,r)
\end{cases}
当 l=r-1 时有边界条件（以下记为边界条件 1）：
\begin{cases}
S(r-1,r)=S(r-1) \\
R(r-1,r)=R(r-1)
\end{cases}
那么，S(1,n+1) 即为答案。
::::info[注]{open}
边界条件的选取多种多样，视具体情况而定，以下提供一种可能的形式（记为边界条件 2）：
\begin{cases}
S(r-1,r)=S(r) \\
R(r-1,r)=R(r)
\end{cases}
如果使用边界条件 2，答案则为 S(0,n)。
在前言的文章中还提到了另外一种边界条件，读者可以到原文查看。
::::
分治计算即可。
分析与证明
::::info[提示]{open}
这一部分你如果熟练掌握了 Binary Splitting 就是一眼的结论。你也可以用其他的方法推导，都不难。可以视情况跳过。
::::
乍一看可能并不能很快地理解这是怎么来的。
但是我们已经知道，这个分治算法有一个父亲——Binary Splitting。
既然这个算法的正确性已经得到证明，我们不妨使用这个来推出本文的分治算法。
Binary Splitting 所求级数的基本形式是这个：
\sum^{n}_{i=1} \frac{P(i)}{R(i)} \prod^{i}_{j=1}\frac{R(j)}{Q(j)}
把 $Q(x)=1$ 代入：
$$
\begin{aligned}
&\sum^{n}_{i=1} \frac{P(i)}{R(i)} \prod^{i}_{j=1}\frac{R(j)}{1}\\
&=\sum^{n}_{i=1} {P(i)}\prod^{i-1}_{j=1}{R(j)}
\end{aligned}
$$
就能得到我们分治算法所求级数的基本形式了。
再看一下 Binary Splitting 计算过程中用到的 $P(l,r)$ 和 $Q(l,r)$：
$$
P(l,r)=\sum^{r}_{i=l} P(i) \cdot (\prod^{i-1}_{j=l}R(j)) \cdot (\prod^{r}_{j=i+1}Q(j))
$$
$$
Q(l,r)=\prod^{r}_{i=l}Q(i)
$$
代入 $Q(x)=1$，发现 $P(l,r)$ 与本文的 $S(l,r)$ 等价。
再看 Binary Splitting 分治过程用的递推式子：
$$
\begin{cases}
P(l,r)=P(l,m) \cdot Q(m+1,r)+R(l,m) \cdot P(m+1,r)\\
Q(l,r)=Q(l,m) \cdot Q(m+1,r)\\
R(l,r)=R(l,m) \cdot R(m+1,r)
\end{cases}
$$
$Q(x)=1$ 那么 $Q(l,r)$ 也等于 $1$，代入 $Q(l,r)=1$，再将 $P(l,r)$ 换成 $S(l,r)$ 即可得到本文分治算法的分治递推式。
于是我们发现并证明了本文的分治算法其实就是 Binary Splitting 的一个扩展。
## 时间复杂度
根据那篇文章的结论，这个算法时间复杂度应为 $O(n \log^3 n)$。
然而其实有一只 $\log$ 是 $\log_{10}$，因为将普通数值转为高精度数值或有一定开销。由于 Binary Splitting 的应用场景是求一类无理数小数点后 $n$ 位，在这个场景下涉及到较复杂的高精度小数计算，开销较大。因而这只 $\log$ 无法省略。
然而本文的模型是面向整数的运算，没有很大的开销，所以只要适当压位，这只 $\log$ 完全可以被优化至忽略不计，因此适当压位后，这个算法复杂度便可降为 $O(n\log^2 n)$。
## 优点和缺点
这个算法实质上是 Binary Splitting 在高精度整数运算方面的具体实例，连高精度小数类都不用实现，结论相对比较好记，所以挺好上手，是暴力的一个很好的替代方式。
那么缺点是什么呢？
**常数特别特别大！！！**
虽然应该不及 [P5050](https://www.luogu.com.cn/problem/P5050)，但是如果使用的 FFT 除了是迭代倍增的 FFT 以外一点优化没有，并且还不压位，想通过例题也是比较吃力的，下文题解部分将会列出一些卡常的方法以供参考。
::::info[常数较大的原因]{open}
其一是因为每次递归需要两次高精度乘法和一次高精度加法。
另外一个重要的原因是答案所处的量级也是很大的。别看 $n$ 上界只有 $1.5 \times 10^5$，本例题最大点的答案却大概处于 $10^{700000}$（$10$ 的 $7\times 10^5$ 次方）量级，$7\times 10^5$ 远比 $1.5 \times 10^5$ 大。
::::
## 例题题解
现在回看例题，想必思路已十分明晰！
注意到：
$$
\begin{aligned}
S&=\sum^{m}_{i=1}i! \\
&=\sum^{m}_{i=1}i\prod^{i-1}_{j=1}j
\end{aligned}
$$
于是可以得出：
$$
\begin{cases}
P(x)=x\\
R(x)=x
\end{cases}
$$
我们就可以使用上述算法分治计算了。
::::info[一些卡常小 tips]
- **压位**。要想干掉 $\log_{10}$，压位必不可少。你可以直接压 $8$ 位。但是为了保留精度，在做 FFT 时要将压 $8$ 位转换成压 $4$ 位，乘完后再转回去。可以参考代码实现。
- **三次变两次优化**。（注：三次变两次优化会让代码在相乘的数值域相差过大的情况下爆精度，平时写高精度题目时慎用。本题作者亲测可行。）
- **预处理单位根**。如果不做这一步的话极有可能炸精度，只能开 `long double`，[尽管可以通过](https://www.luogu.com.cn/record/263020936)，跑出来的结果却不甚令人满意。
- 对于比较小的高精度整数，我们直接**使用小常数 $O(n^2)$ 算法处理**。
做完上述优化之后已经足够了，最大点 340ms 左右，虽然仍然有较大优化空间，比如搭上指令集，隐式位逆序置换等，这里不给出了，因为~~我不会~~，并且没有太大必要。
::::
::::success[AC code & 提交记录]
代码有点长，本文只给出大致的代码框架。想看完整代码[点这里](https://www.luogu.me/paste/bqw46lco)，完整代码包含上面说到的大部分卡常技巧的一个参考实现。
（完整代码码风实在是太猎奇了，不过至少能读。）
```cpp
//edit by devc++
//by gitxiaozheng/dakada
#include<bits/stdc++.h>
using namespace std;
#define L long long
#define int unsigned int
#define bigint vector<int>
#define space putchar(' ')
#define enter putchar('\n')
inline L qcin();
void qcout(L);
const int base=1e8,FFTbase=1e4;
const double pi=acos(-1);
bigint P,R;
//实现FFT多项式乘法+高精度运算的代码省略
void bs(int l,int r,bigint& PT,bigint &RT){
if(l==r-1){
//这里使用边界条件 1。
PT.push_back(l);
RT.push_back(l);
return;
}
/*
P1:S(l,m)
P2:S(m+1,r)
R1:R(l,m)
R2:R(m+1,r)
*/
bigint P1,P2,R1,R2;
int mid=(l+r)>>1;
bs(l,mid,P1,R1);
bs(mid,r,P2,R2);
RT=bigintMul(R1,R2);
PT=bigintadd(bigintMul(P2,R1),P1);
return;
}
#undef int
signed main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout.tie(nullptr);
int n;
cin>>n;
bs(1,n+1,P,R);
for(int i=P.size()-1;i>=0;i--){
if(i==P.size()-1) cout<<P[i];
else cout<<setw(8)<<setfill('0')<<P[i];
}
return 0;
}
//前面的代码有快读快写函数的声明，然而这篇代码没有用到快读快写。。。
```
这是[提交记录](https://www.luogu.com.cn/record/263468065)。
::::
## 后记
类似于 [P10967](https://www.luogu.com.cn/problem/P10967)，这道题有加强版必然有更强的加强版：把 $n$ 开到 $10^9$ 但是要对答案取模，详情请看[这里](https://www.luogu.com.cn/article/z9cnl0pw)。
同时，还有一种加强方式：$n\le 10^6$，只要卡常卡到能在 500ms 以内（这是**所有测试点的时间总和**而不是单个测试点所需时间）通过 [P1919](https://www.luogu.com.cn/problem/P1919)，那么以本题的时限是完全可以通过的。