[学习笔记]程设作业的一些YY，高精度除法/FFT/Karastuba/牛顿迭代
高精度除法
最近程序设计的作业有个高精度除法，正好我抽到写这题的解题报告，借此机会稍微研究了一下相关的做法，于是就有了这些东西。
这篇博客大概是围绕高精乘法/除法展开，从朴素做法出发，简单介绍一下分治乘法（Karatsuba算法）及一些优化、经典\(O(n\log n)\)复杂度的FFT、牛顿迭代法实现高精度除法等等。FFT更多地是整理一个大概的框架和自己的口胡，关于FFT网上有许多优秀的博客（比如我就是看@miskcoo的那篇《从多项式乘法到快速傅里叶变换》学的），我也就没有更多地记录。牛顿迭代法则是参考了倪泽堃2011年的课件《理性愉悦——高精度数值计算》里提到的方法。
一些朴素做法
要计算\(n\)位大数\(A\)和\(m\)位大数\(B\)的商\(A/B\)，一个很暴力的做法是用减法模拟除法，如果\(A<B\)则结果直接返回0，否则如果\(A\geq B\)，每次选取\(A\)的前\(m\)位\(A'\)，和\(B\)做减法，直到不能做减法为止，注意这样没对\(A',B\)至多做\(9\)次减法，每次的复杂度是\(O(m)\)，一共做\(O(n-m)\)次，\(n,m\)接近的话复杂度依然是\(O(n^2)\)级别的。
另一种牺牲时间的暴力做法是直接考虑二分答案，假设\(\lfloor A/B \rfloor\)的答案为\(t\)，则一定有不等式：\(Bt\leq A<B(t+1)\)，这里的\(t\)就出现了单调性，每次乘法的复杂度是\(O(n^2)\)，而二分的代价是\(O(\log 值域)\)级别的，注意这里值域其实是\(O(10^n)\)而不是\(O(n)\)，所以总复杂度会达到\(O(n^3)\)，对于范围不大的高精度除法可以考虑直接用这种三方的做法实现（比如这次程设作业，\(n\)只有500，\(O(n^3)\)足够通过），实现上会比上一种除法来得方便。
下面提到的Karatsuba和FFT可以优化乘法的效率，不过最快也只能把这种方法优化到\(O(n^2\log n)\)，还不如暴力的\(O(n^2)\)。但这并不失为一个思路——后面我们会顺着这个思路给出一个复杂度\(O(n\log n)\)的牛顿迭代的做法。
分治乘法——Karatsuba
原本是在Miskcoo的博客学FFT来着，刚好看见这个东西，看起来挺简单的就试着实现了一下：
原理很简单，我们要算\(x\times y\)的结果，假设它们中最大的数位长度是\(n\)，我们选取一个\(m\)（一般就是\(n/2\)，不过可能涉及到取整的问题，我就再拿一个变量）成\(x=A\times 10^{m}+B,y=C\times 10^{m}+D\)，接着\(x\times y=AC\times 10^{2m}+(AD+BC)\times 10^m +BD\)，这样看起来是要四次乘法，复杂度\(T(n)=4T(\frac{n}{2})+O(n)\)，对应主定理的一条：存在一个\(\xi>0:O(n)=O(n^{\log_2 4-\xi})\)，则\(T(n)=\Theta(n^{\log_2 4})\)，这样子做总复杂度还是\(O(n^2)\)，和暴力相比没什么变化。
不过如果我们改写\((AD+BC)=(A+B)(C+D)-(AC+BD)\)的话，就可以把乘法次数减少一次，于是新的复杂度\(T_2(n)=3T_2(\frac{n}{2})+O(n)\)，复杂度\(O(n^{\log _2 3})\)，大约是\(O(n^{1.58})\)。
实现部分大致如下：split是把数拿去拆分，pow10(A,d)就是把数字\(A\)乘个\(10^{POWER*d}\)，不过我写的常数非常大…跑的也很慢。
Big kara(Big x,Big y,int d){
if(x<y)swap(x,y);
if(x.len<=1&&y.len<=1){
Big res=Big(0);
res.a[0]=x.a[0]*y.a[0];
if(res.a[0]>=BASE){
res.a[++res.len-1]=res.a[0]/(BASE+1);
res.a[0]%=BASE+1;
}
return res;
}
int n=x.len,m;m=n/2;
Big A,B,C,D,AC,BD,H;
B=split(x,0,m-1);A=split(x,m,x.len-1);
D=split(y,0,m-1);C=split(y,m,y.len-1);
AC=kara(A,C,d+1);BD=kara(B,D,d+1);
H=kara(A+B,C+D,d+1)-(AC+BD);
AC=pow10(AC,2*m);H=pow10(H,m);
Big res=AC+H+BD;
while(res.len>1&&res.a[res.len-1]==0)res.len--;
return res;
}
之后问了一下@Moebius Meow学姐为什么会这么慢…发现这段代码各种问题…
一个很大的问题是每层递归都需要申请一段新的内存，内存大小是\(O(n)\)的而不是\(O(1)\)的（这里的每个Big就是一个数组），这就会让效率很低。常见的优化（经典的分治做矩阵乘法也可以这么优化）就是开个全局变量\(A\)当内存池用，每层递归把声明新变量改成直接在内存池上操作，算完了一层就直接把这块内存全释放，伪代码：
kara{
申请res的空间
记录内存池起点
申请临时的各种内存以及进行递归
计算res
重置内存池起点（弹出子递归树的内存）
返回res
}
这样操作最大需要的内存对应着的就是递归树的一条从根到叶子的链，也就是\(M(n)=M(\frac{n}{2})+kn\)，解得\(M(n)\approx2kn\)。
另一个优化是在规模比较小的时候考虑\(O(n^2)\)的暴力，大的时候再分治，效率会快许多。
不过由于今天的重点不在于这个分治乘法，加上这个作者时间不够（可恶，作业好多，学不完了.jpg），就没有具体实现了（x）
多项式乘法——FFT
具体的知识还是看miskcoo的博客：
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
大致思路就是我们考虑计算两个多项式\(\begin{aligned}F(x)=\sum_{i=0}^na_ix^i\end{aligned}\)和\(\begin{aligned}G(x)=\sum_{i=0}^m b_ix^i\end{aligned}\)的乘积\(R(x)\)，进位之后把\(x=10\)代入\(R(x)\)就是我们要的高精度乘法的答案。我们常见的是上面这种系数表示，直接计算每一项的系数：\(R(x)=F(x)G(x)\)，对应的第\(k\)项的系数\(\begin{aligned}c_k=\sum_{i+j=k}a_ib_j\end{aligned}\)，这样会达到\(O(n^2)\)的复杂度。
但是如果我们手上的是两个点值表示的多项式：\(n\)次多项式用\(n+1\)个点\((x_k,f(x_k))(k=0,\dots,n)\)确定一个\(n\)次多项式，就能在\(O(n)\)的时间内实现点值表示的乘法（直接\(R(x_k)=F(x_k)G(x_k)\)）。朴素的系数表示到点值表示的算法需要\(O(n^2)\)的时间（如果没什么优化，这些\(x\)互相没什么关系，我们需要对\(O(n)\)个\(x\)分别求值），点值表示转化到系数表示方法就很多，比如看成解方程\(Xa=Y\)，用高斯消元在\(O(n^3)\)的时间解决，以及还有\(O(n^2)\)的拉格朗日插值，不过如果只是这样就和朴素的做法没什么差别了。
而FFT要做的就是这件事情：离散傅里叶变换DFT实现系数表示到点值表示的转换，离散傅里叶逆变换IDFT实现点值到系数的转换，两个转换都能做到\(O(n\log n)\)的复杂度。
DFT是先把\(n\)进行适当地扩大（先扩大成比原来\(n\)大的最小的2的次幂，接着再扩大一倍），然后再巧妙地选取\(n\)个自变量——\(n\)个\(n\)次单位根\(\omega_n^k=\cos(\frac{2k\pi}{n})+i\sin(\frac{2k\pi}{n}),(k=0,\dots,n-1)\)，来代入一个\(n-1\)次多项式，利用单位根的性质实现分治求解：考虑我们需要计算一个\(\begin{aligned}A(\omega_n^k)=\sum_{i=0}^{n-1} a_i\omega_n^{ki}\end{aligned}\)，\(\omega_n^k\)的次数按照奇偶性拆开：
\(\begin{aligned}A(\omega_n^k)=\sum_{i=0}^{n/2-1}a_{2i}\omega_n^{2ki}+\omega_n^k\sum_{i=0}^{n/2-1}a_{2i+1}\omega_n^{2ki}\end{aligned}\)，接着是利用\(\omega_n^{2k}=(e^{2k\pi/n})^2=(e^{2k\pi/(n/2)})=\omega_{n/2}^k\)，于是有\(\begin{aligned}A(\omega_n^k)=\sum_{i=0}^{n/2-1}a_{2i}\omega_{n/2}^{ki}+\omega_n^k\sum_{i=0}^{n/2-1}a_{2i+1}\omega_{n/2}^{ki}\end{aligned}\)，这两个\(\sum\)的形式又和\(A(\omega_n^k)\)很相似了，不过这样看起来问题规模好像没有怎么缩小（还是要对每个\(k\)递归计算），但是如果再利用\(\omega_{n}^{k+\pi/2}=\omega_{n}^k\cdot\omega_n^{\pi/2}=-\omega_{n}^k\)，就能用\(k=0,\dots,n/2-1\)的答案计算\(k=n/2,\dots,n-1\)的答案，这样复杂度就是\(T(n)=2T(n/2)+O(n),T(n)=O(n\log n)\)了。
而IDFT则和FFT类似，数学上可以证明，只要对点值表示得到的那\(n\)个点再做一次FFT，并把FFT里的\(\omega_n^k\)换成\(\omega_n^{-k}\)，就能得到对应的系数了。
这里就先贴个代码。
void init_eps(int p){
double ang=2.0*PI/p;
rep(i,0,p-1)eps[i]=cmd(cos(ang*i),sin(ang*i));
rep(i,0,p-1)inv_eps[i]=conj(eps[i]);
}
void FFT(int n,cmd *x,cmd *w){
for(int i=0,j=0;i<n;++i){
if(i>j)swap(x[i],x[j]);
for(int l=(n>>1);(j^=l)<l;l>>=1);
}
for(int i=2;i<=n;i<<=1){
int m=i>>1;
for(int j=0;j<n;j+=i)rep(k,0,m-1){
cmd z=x[j+m+k]*w[n/i*k];
x[j+m+k]=x[j+k]-z;
x[j+k]+=z;
}
}
}
Big operator*(const Big &x,const Big &y){
Big res=Big(0);
int p=1;
while(p<max(x.len,y.len))p<<=1;p<<=1;
init_eps(p);
rep(i,0,p-1)f[i]=g[i]=cmd(0,0);//initialize
rep(i,0,p-1)f[i]=x.a[i];
rep(i,0,p-1)g[i]=y.a[i];
FFT(p,f,eps);FFT(p,g,eps);
rep(i,0,p-1)f[i]*=g[i];
FFT(p,f,inv_eps);
res.len=p;
rep(i,0,p-1)res.a[i]=(int)(0.5+f[i].real()/p);
rep(i,0,p-1)if(res.a[i]>BASE){
res.a[i+1]+=res.a[i]/(BASE+1);
res.a[i]=res.a[i]-res.a[i]/(BASE+1)*(BASE+1);
}
while(res.len>1&&res.a[res.len-1]==0)res.len--;
return res;
}
注意迭代实现的部分
for(int i=0,j=0;i<n;++i){
if(i>j)swap(x[i],x[j]);
for(int l=(n>>1);(j^=l)<l;l>>=1);
}
for(int i=2;i<=n;i<<=1){
int m=i>>1;
for(int j=0;j<n;j+=i)rep(k,0,m-1){
cmd z=x[j+m+k]*w[n/i*k];
x[j+m+k]=x[j+k]-z;
x[j+k]+=z;
}
}
前面一个for循环是在把下标分组，接着第一层i的循环是在从下往上迭代，m则是这一层子问题的长度，接着j的那层循环则是遍历这一层的所有子问题，k的一层则是对每一个子问题扫一遍，里面的东西则是像滚动数组那样合并答案。
以及我在用FFT写高精乘法的时候企图顺便写一个压位，交上去疯狂WA…于是又去问了@Moebius Meow学姐，学姐一下子就看出来问题，考虑一个长度为\(L\)的数字压\(x\)位，做FFT之后的一位最大就可以达到\(\frac{L}{x}10^{2x}\)（每一项对应的都是一个卷积），double有效位数大概是14位，一般FFT处理的高精乘法的长度可以达到\(10^6\)，就算我压了2位，最大也有\(10^9\)多了，精度就已经很危险了。似乎一个比较好用的做法是FNTT（快速数论变换）之类的东西，不过作者还不会。
牛顿迭代法
上面的二分的一个瓶颈在于收敛得还是不够快，要说什么收敛速度很快的迭代方式，那肯定就是牛顿迭代啦…
考虑任意一个在定义域\(D\)上连续且\(n\)阶可导的函数\(f\)，在\(x_0(x_0\in D)\)一点的Taylor展开：
\(\begin{aligned}f(x)=\sum_{i=0}^{\infty} f^{(i)}(x_0)\frac{(x-x_0)^i}{i!}\end{aligned}\)，考虑令\(f(x)=0\)，任意函数的根就转化成了求一个多项式的根，但毕竟右边是个无穷项的多项式，求精确解是别想了（何况高次多项式也没有解析解）。但其实上我们只要取Taylor展开的前几项做近似，来逼近方程的根，这样往往也能达到很好的效果：牛顿迭代就是取前两项来逼近。
考虑近似\(f(x)\approx f(x_0)+f'(x_0)(x-x_0)=0\)，得\(\begin{aligned}x=x_0-\frac{f(x_0)}{f'(x_0)}\end{aligned}\)，对应的迭代方式是：\(\begin{aligned}x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}\end{aligned}\)，放在这里，我们如果要计算一个数\(B\)的倒数\(x\)：\(\frac{1}{x}=B\)，就可以考虑\(f(x)=\frac{1}{x}-B\)，对应的牛顿迭代\(\begin{aligned}x_{k+1}=x_k-\frac{\frac{1}{x_k}-B}{-\frac{1}{x_k^2}}=x_k+x_k-Bx_k^2=x_k(2-Bx_k)\end{aligned}\)，容易验证这样可以收敛到我们想要的结果：
考虑构造数列\(\{a_n\}:0<a_1<\frac{1}{B},a_{n+1}=a_n(2-Aa_n)\)，定义函数\(g(x)=2x-Ax^2(0\leq x\leq \frac{1}{A})\)，\(g'(x)=2(1-Ax)>0\)，\(g(x)\)在定义域上递增，易得对应数列递增，同时\(a_{n+1}=a_n(2-Aa_n)>a_n(2-1)>0\)，\(a_{n+1}=a_n(2-Aa_n)<\frac{2}{A}\)，有界，单调有界准则告诉我们\(a_n\)收敛，设极限是\(X》0\)，则\(X=2X-AX^2,X=\frac{1}{A}\)。
那么这个东西收敛得有多快呢？考虑再构造一个\(b_n=\frac{1}{B}-a_n\)来描述迭代\(n\)次之后，我们用来迭代的\(a_n\)距离它的极限还差多少，有\(a_n=\frac{1}{B}-b_n\)，可以推得\(b_{n+1}=\frac{1}{B}-\frac{2}{B}+2b_n+B(\frac{1}{B^2}+b_n^2-\frac{2b_n}{B})=Bb_n^2=B(Bb_1)^{2^n-1}\)。
\(Bb_1\)很显然小于1，注意到它是二阶收敛的（次数是\(2^n\)级别！）这是什么概念？每多迭代一次，有效的数位就会多一倍，要得到一个比较令人满意的结果，大约只需要迭代\(O(\log\log n)\)次，例如值域如果是\(10^{10^6}\)，二分需要大约\(10^6\)次，而牛顿迭代大约只要十几次就可以了！
回到原来的问题，我们先用牛顿迭代计算出\(1/B\)的结果，再计算\(A\times(1/B)\)即可，求\(1/B\)的时候为了避免做实数的高精度运算，可以考虑类似前面FFT那样，把要计算的结果看成是一个关于\(x\)的多项式，只不过这里是\(x=0.1\)的情形，最后取小数点后一定长度的结果作为答案，再去和\(A\)乘。同时注意我们前面对迭代的分析，数列\(a_n\)其实是严格小于我们的极限\(1/B\)的，所以这会导致我们在遇到刚刚好能整除的时候，结果可能比正确答案小了1，所以对于算出来的一个答案\(t\)，我们还要另外判断一下正确答案是\(t\)还是\(t+1\)。
同时如果一开始开足够长的数来迭代，需要迭代\(O(\log \log 10^n)=O(\log n)\)次，复杂度是\(O(n\log^2n)\)的，但是如果我们一边迭代一边扩展长度，每次把长度扩展一倍，这样时间\(T(n)\approx T(n/2)+O(n\log n)\)，这里\(O(n\log n)\)占主导，对应主定理的存在\(\epsilon>0:O(n\log n)=\Omega(n^{\log_2 1+\epsilon})=\Omega(n^\epsilon)\)，\(T(n)=O(n\log n)\)，这样又降了一个\(\log\)，最后一次乘法也用FFT实现。
就这样，我们就获得了做高精度除法的\(O(n\log n)\)的算法！！
Poly Inverse(Poly A){
Poly B,C;B.resize(2);B[1]=100/(A[0]*10+(A.size()>1?A[1]:0));
for(int s=1;s<=10;++s){//迭代次数
C.resize(1<<s),B.resize(1<<s);
for(int i=0;i<min(1ll<<s,(int)A.size());++i) C[i]=A[i];
for(int i=min(1ll<<s,(int)A.size());i<1<<s;++i) C[i]=0;
C=B*C;
for(int i=1;i<C.size();++i) C[i]=-C[i];
C[0]=2-C[0];
for(int i=C.size()-1;i;--i){
C[i-1]+=(C[i]-9)/10;
C[i]=(C[i]+10)%10;
}
B=B*C;
B.resize(1<<s);
}
return B;
}
结束啦
感谢前队友/室友lzx花了一晚上的时间教会了我FFT的原理，以及@Moebius Meow细心地解答我的各种问题。
