[学习笔记]程设作业的一些YY,高精度除法/FFT/Karastuba/牛顿迭代
高精度除法
最近程序设计的作业有个高精度除法,正好我抽到写这题的解题报告,借此机会稍微研究了一下相关的做法,于是就有了这些东西。
这篇博客大概是围绕高精乘法/除法展开,从朴素做法出发,简单介绍一下分治乘法(Karatsuba算法)及一些优化、经典
一些朴素做法
要计算
另一种牺牲时间的暴力做法是直接考虑二分答案,假设
下面提到的Karatsuba和FFT可以优化乘法的效率,不过最快也只能把这种方法优化到
分治乘法——Karatsuba
原本是在Miskcoo的博客学FFT来着,刚好看见这个东西,看起来挺简单的就试着实现了一下:
原理很简单,我们要算
不过如果我们改写
实现部分大致如下:split是把数拿去拆分,pow10(A,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学姐为什么会这么慢…发现这段代码各种问题…
一个很大的问题是每层递归都需要申请一段新的内存,内存大小是
kara{ 申请res的空间 记录内存池起点 申请临时的各种内存以及进行递归 计算res 重置内存池起点(弹出子递归树的内存) 返回res }
这样操作最大需要的内存对应着的就是递归树的一条从根到叶子的链,也就是
另一个优化是在规模比较小的时候考虑
不过由于今天的重点不在于这个分治乘法,加上这个作者时间不够(可恶,作业好多,学不完了.jpg),就没有具体实现了(x)
多项式乘法——FFT
具体的知识还是看miskcoo的博客:
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
大致思路就是我们考虑计算两个多项式
但是如果我们手上的是两个点值表示的多项式:
而FFT要做的就是这件事情:离散傅里叶变换DFT实现系数表示到点值表示的转换,离散傅里叶逆变换IDFT实现点值到系数的转换,两个转换都能做到
DFT是先把
而IDFT则和FFT类似,数学上可以证明,只要对点值表示得到的那
这里就先贴个代码。
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学姐,学姐一下子就看出来问题,考虑一个长度为
牛顿迭代法
上面的二分的一个瓶颈在于收敛得还是不够快,要说什么收敛速度很快的迭代方式,那肯定就是牛顿迭代啦…
考虑任意一个在定义域
考虑近似
考虑构造数列
那么这个东西收敛得有多快呢?考虑再构造一个
回到原来的问题,我们先用牛顿迭代计算出
同时如果一开始开足够长的数来迭代,需要迭代
就这样,我们就获得了做高精度除法的
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细心地解答我的各种问题。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
· Linux系列:如何用heaptrack跟踪.NET程序的非托管内存泄露
· 开发者必知的日志记录最佳实践
· SQL Server 2025 AI相关能力初探
· 震惊!C++程序真的从main开始吗?99%的程序员都答错了
· 【硬核科普】Trae如何「偷看」你的代码?零基础破解AI编程运行原理
· 单元测试从入门到精通
· 上周热点回顾(3.3-3.9)
· winform 绘制太阳,地球,月球 运作规律