[学习笔记]程设作业的一些YY,高精度除法/FFT/Karastuba/牛顿迭代

高精度除法

最近程序设计的作业有个高精度除法,正好我抽到写这题的解题报告,借此机会稍微研究了一下相关的做法,于是就有了这些东西。

这篇博客大概是围绕高精乘法/除法展开,从朴素做法出发,简单介绍一下分治乘法(Karatsuba算法)及一些优化、经典O(nlogn)复杂度的FFT、牛顿迭代法实现高精度除法等等。FFT更多地是整理一个大概的框架和自己的口胡,关于FFT网上有许多优秀的博客(比如我就是看@miskcoo的那篇《从多项式乘法到快速傅里叶变换》学的),我也就没有更多地记录。牛顿迭代法则是参考了倪泽堃2011年的课件《理性愉悦——高精度数值计算》里提到的方法。

一些朴素做法

要计算n位大数Am位大数B的商A/B,一个很暴力的做法是用减法模拟除法,如果A<B则结果直接返回0,否则如果AB,每次选取A的前mA,和B做减法,直到不能做减法为止,注意这样没对A,B至多做9次减法,每次的复杂度是O(m),一共做O(nm)次,n,m接近的话复杂度依然是O(n2)级别的。

另一种牺牲时间的暴力做法是直接考虑二分答案,假设A/B的答案为t,则一定有不等式:BtA<B(t+1),这里的t就出现了单调性,每次乘法的复杂度是O(n2),而二分的代价是O(log)级别的,注意这里值域其实是O(10n)而不是O(n),所以总复杂度会达到O(n3),对于范围不大的高精度除法可以考虑直接用这种三方的做法实现(比如这次程设作业,n只有500,O(n3)足够通过),实现上会比上一种除法来得方便。

下面提到的Karatsuba和FFT可以优化乘法的效率,不过最快也只能把这种方法优化到O(n2logn),还不如暴力的O(n2)。但这并不失为一个思路——后面我们会顺着这个思路给出一个复杂度O(nlogn)的牛顿迭代的做法。

分治乘法——Karatsuba

原本是在Miskcoo的博客学FFT来着,刚好看见这个东西,看起来挺简单的就试着实现了一下:

原理很简单,我们要算x×y的结果,假设它们中最大的数位长度是n,我们选取一个m(一般就是n/2,不过可能涉及到取整的问题,我就再拿一个变量)成x=A×10m+B,y=C×10m+D,接着x×y=AC×102m+(AD+BC)×10m+BD,这样看起来是要四次乘法,复杂度T(n)=4T(n2)+O(n),对应主定理的一条:存在一个ξ>0:O(n)=O(nlog24ξ),则T(n)=Θ(nlog24),这样子做总复杂度还是O(n2),和暴力相比没什么变化。

不过如果我们改写(AD+BC)=(A+B)(C+D)(AC+BD)的话,就可以把乘法次数减少一次,于是新的复杂度T2(n)=3T2(n2)+O(n),复杂度O(nlog23),大约是O(n1.58)

实现部分大致如下:split是把数拿去拆分,pow10(A,d)就是把数字A乘个10POWERd,不过我写的常数非常大…跑的也很慢。

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(n2)+kn,解得M(n)2kn

另一个优化是在规模比较小的时候考虑O(n2)的暴力,大的时候再分治,效率会快许多。

不过由于今天的重点不在于这个分治乘法,加上这个作者时间不够(可恶,作业好多,学不完了.jpg),就没有具体实现了(x)

多项式乘法——FFT

具体的知识还是看miskcoo的博客:

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

大致思路就是我们考虑计算两个多项式F(x)=i=0naixiG(x)=i=0mbixi的乘积R(x),进位之后把x=10代入R(x)就是我们要的高精度乘法的答案。我们常见的是上面这种系数表示,直接计算每一项的系数:R(x)=F(x)G(x),对应的第k项的系数ck=i+j=kaibj,这样会达到O(n2)的复杂度。

但是如果我们手上的是两个点值表示的多项式:n次多项式用n+1个点(xk,f(xk))(k=0,,n)确定一个n次多项式,就能在O(n)的时间内实现点值表示的乘法(直接R(xk)=F(xk)G(xk))。朴素的系数表示到点值表示的算法需要O(n2)的时间(如果没什么优化,这些x互相没什么关系,我们需要对O(n)x分别求值),点值表示转化到系数表示方法就很多,比如看成解方程Xa=Y,用高斯消元在O(n3)的时间解决,以及还有O(n2)的拉格朗日插值,不过如果只是这样就和朴素的做法没什么差别了。

而FFT要做的就是这件事情:离散傅里叶变换DFT实现系数表示到点值表示的转换,离散傅里叶逆变换IDFT实现点值到系数的转换,两个转换都能做到O(nlogn)的复杂度。

DFT是先把n进行适当地扩大(先扩大成比原来n大的最小的2的次幂,接着再扩大一倍),然后再巧妙地选取n个自变量——nn次单位根ωnk=cos(2kπn)+isin(2kπn),(k=0,,n1),来代入一个n1次多项式,利用单位根的性质实现分治求解:考虑我们需要计算一个A(ωnk)=i=0n1aiωnkiωnk的次数按照奇偶性拆开:

A(ωnk)=i=0n/21a2iωn2ki+ωnki=0n/21a2i+1ωn2ki,接着是利用ωn2k=(e2kπ/n)2=(e2kπ/(n/2))=ωn/2k,于是有A(ωnk)=i=0n/21a2iωn/2ki+ωnki=0n/21a2i+1ωn/2ki,这两个的形式又和A(ωnk)很相似了,不过这样看起来问题规模好像没有怎么缩小(还是要对每个k递归计算),但是如果再利用ωnk+π/2=ωnkωnπ/2=ωnk,就能用k=0,,n/21的答案计算k=n/2,,n1的答案,这样复杂度就是T(n)=2T(n/2)+O(n),T(n)=O(nlogn)了。

而IDFT则和FFT类似,数学上可以证明,只要对点值表示得到的那n个点再做一次FFT,并把FFT里的ωnk换成ωnk,就能得到对应的系数了。

这里就先贴个代码。

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之后的一位最大就可以达到Lx102x(每一项对应的都是一个卷积),double有效位数大概是14位,一般FFT处理的高精乘法的长度可以达到106,就算我压了2位,最大也有109多了,精度就已经很危险了。似乎一个比较好用的做法是FNTT(快速数论变换)之类的东西,不过作者还不会。

牛顿迭代法

上面的二分的一个瓶颈在于收敛得还是不够快,要说什么收敛速度很快的迭代方式,那肯定就是牛顿迭代啦…

考虑任意一个在定义域D上连续且n阶可导的函数f,在x0(x0D)一点的Taylor展开:
f(x)=i=0f(i)(x0)(xx0)ii!,考虑令f(x)=0,任意函数的根就转化成了求一个多项式的根,但毕竟右边是个无穷项的多项式,求精确解是别想了(何况高次多项式也没有解析解)。但其实上我们只要取Taylor展开的前几项做近似,来逼近方程的根,这样往往也能达到很好的效果:牛顿迭代就是取前两项来逼近。

考虑近似f(x)f(x0)+f(x0)(xx0)=0,得x=x0f(x0)f(x0),对应的迭代方式是:xk+1=xkf(xk)f(xk),放在这里,我们如果要计算一个数B的倒数x1x=B,就可以考虑f(x)=1xB,对应的牛顿迭代xk+1=xk1xkB1xk2=xk+xkBxk2=xk(2Bxk),容易验证这样可以收敛到我们想要的结果:

考虑构造数列{an}:0<a1<1B,an+1=an(2Aan),定义函数g(x)=2xAx2(0x1A)g(x)=2(1Ax)>0g(x)在定义域上递增,易得对应数列递增,同时an+1=an(2Aan)>an(21)>0an+1=an(2Aan)<2A,有界,单调有界准则告诉我们an收敛,设极限是X0,则X=2XAX2,X=1A

那么这个东西收敛得有多快呢?考虑再构造一个bn=1Ban来描述迭代n次之后,我们用来迭代的an距离它的极限还差多少,有an=1Bbn,可以推得bn+1=1B2B+2bn+B(1B2+bn22bnB)=Bbn2=B(Bb1)2n1

Bb1很显然小于1,注意到它是二阶收敛的(次数是2n级别!)这是什么概念?每多迭代一次,有效的数位就会多一倍,要得到一个比较令人满意的结果,大约只需要迭代O(loglogn)次,例如值域如果是10106,二分需要大约106次,而牛顿迭代大约只要十几次就可以了!

回到原来的问题,我们先用牛顿迭代计算出1/B的结果,再计算A×(1/B)即可,求1/B的时候为了避免做实数的高精度运算,可以考虑类似前面FFT那样,把要计算的结果看成是一个关于x的多项式,只不过这里是x=0.1的情形,最后取小数点后一定长度的结果作为答案,再去和A乘。同时注意我们前面对迭代的分析,数列an其实是严格小于我们的极限1/B的,所以这会导致我们在遇到刚刚好能整除的时候,结果可能比正确答案小了1,所以对于算出来的一个答案t,我们还要另外判断一下正确答案是t还是t+1

同时如果一开始开足够长的数来迭代,需要迭代O(loglog10n)=O(logn)次,复杂度是O(nlog2n)的,但是如果我们一边迭代一边扩展长度,每次把长度扩展一倍,这样时间T(n)T(n/2)+O(nlogn),这里O(nlogn)占主导,对应主定理的存在ϵ>0:O(nlogn)=Ω(nlog21+ϵ)=Ω(nϵ)T(n)=O(nlogn),这样又降了一个log,最后一次乘法也用FFT实现。

就这样,我们就获得了做高精度除法的O(nlogn)的算法!!

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细心地解答我的各种问题。

posted @   yoshinow2001  阅读(1180)  评论(1编辑  收藏  举报
编辑推荐:
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
· Linux系列:如何用heaptrack跟踪.NET程序的非托管内存泄露
· 开发者必知的日志记录最佳实践
· SQL Server 2025 AI相关能力初探
阅读排行:
· 震惊!C++程序真的从main开始吗?99%的程序员都答错了
· 【硬核科普】Trae如何「偷看」你的代码?零基础破解AI编程运行原理
· 单元测试从入门到精通
· 上周热点回顾(3.3-3.9)
· winform 绘制太阳,地球,月球 运作规律
点击右上角即可分享
微信分享提示