C++ 大整数运算 高精度除法
前言
这篇文章主要是对于大整数类的设计过程中,如何实现并改进长除法(模拟竖式法)的一个总结。
高精度除法的分类和比较
虽然有些文章在讨论大整数的除法运算时,喜欢分成高精度除以高精度和高精度除以低精度(短除法)两种情况并分开实现,但在本文中将不做这种区分。事实上,短除法就是长除法的一个缩略。
按照wiki百科,除法大致可以分为以下两类:
- 慢速算法:每步只能确定一位的算法。长除法属于此类,时间复杂度为O(N2)。
- 快速算法:每步可以确定多位的算法。牛顿迭代法属于此类,配合用FFT实现的乘法,时间复杂度为O(NlogN)。
一些讨论
这里讨论的是为了实现大整数类所必需的一些基础。厘清它们对提升算法的效率是有意义的。
1. 关于底层
标准C++的语法中基本的数据类型都带有基本的四则运算,这里的底层主要指硬件或者编译器。
要实现大整数的除法运算,至少需要用到加、减两种基本运算,乘、除和取余虽然不是必须的,但我们仍然假定底层(硬件或者软件模拟)已经实现,因为后面的优化将用到它们。
2. 关于数据结构
数据结构有多种实现方式,互有优劣。数组最简单,但长度固定,需要考虑溢出。STL的vector不用担心溢出,但不是线程安全的。链表线程安全,但操作复杂且效率低下。
本文中我们将使用数组,其中数组的低位存储数值的低位,数组的最高位存储数值的符号。
3. 关于数据类型
数据类型的选择,取决于存储效率和处理能力。用1个int型存储1个十进制有效数字显然是低效的,但用long型存储9个十进制有效数字也未必最佳(假如系统是16位的,则long型的四则运算本身就需要软件模拟)。有些类型的大小也依赖于编译器的实现,这也是在设计过程中需要考虑的。
本文先从1个char型存储1位十进制有效数字开始设计,暂将优化(压位)放到后面讨论。这样的优点是可以使用字符串进行存储,在输入输出上有所便宜。
长除法的实现
1. 除法的基本原理
本质上,慢速算法都是通过减法来实现除法。不考虑效率的话,下面的代码就可以实现除法。
Integer Integer::operator/(const Integer& itg)const{//...//prod=*this;//div=itg;result=zero;while(!(prod<div)){prod=prod-div;result=result+one;}//...
}
显然,这个实现的时间复杂度是指数级的。现实中大概只有幼儿园的小朋友会使用。
2. 优化策略
如何高效率地完成减法,是优化的一个方向。直观的办法是通过减去除数的某个已知的倍数,来达到提高效率的目的。下面的代码是一个最简单的实现。
Integer Integer::operator/(const Integer& itg)const{//...//prod=*this;//div=itg;result=zero;quot=one;base=div;while(!(prod<div)){if(prod<base){quot=one;base=div;}else{prod=prod-base;result=result+quot;quot=quot+quot;base=base+base;}}//...
}
这个实现的时间复杂度已经达到O(N2),但无奈常数太大。因为每次采用的是翻倍,所以这里计算N是以2为底取的对数,而常用的笔算方法是以10为底取的对数。而且每次计算都要调用大整数的加减法,这也是一笔不小的开销。
3. 长除法的直减版本
下面的代码是采用移位策略的实现,已是相当接近于列竖式笔算的版本。
Integer Integer::operator/(const Integer& itg)const{Integer result=0,prod,div,base;int loop,offset,quot,reg,trans=0;/*处理除0异常*/prod=*this;if(*(prod.p+Length-1)=='-') *(prod.p+Length-1)='+';div=itg;if(*(div.p+Length-1)=='-') *(div.p+Length-1)='+';while(!(prod<div)){base=div;quot=1;//确定移位的具体值loop=prod.index-1;offset=base.index-1;while(*(prod.p+loop)==*(base.p+offset)){if(offset>0){loop--;offset--;}else break;}if(char2int(*(prod.p+loop))>=char2int(*(base.p+offset)))offset=prod.index-base.index;elseoffset=prod.index-base.index-1;//移位if(offset>0){for(loop=base.index-1;loop>=0;loop--)*(base.p+loop+offset)=*(base.p+loop);for(loop=0;loop<offset;loop++)*(base.p+loop)=int2char(0);base.index=base.index+offset;quot=quot+offset;}if(prod.index>base.index) *(base.p+prod.index-1)=int2char(0);//结果初始化,仅执行一次if(quot>result.index){for(loop=0;loop<quot;loop++)*(result.p+loop)=int2char(0);result.index=quot;}while(!(prod<base)){//执行减法for(loop=offset;loop<prod.index;loop++){reg=char2int(*(prod.p+loop))-char2int(*(base.p+loop));if(trans){reg=reg-1;trans=0;}if(reg<0){reg=reg+Radix;trans=1;}*(prod.p+loop)=int2char(reg);}//清除前导0loop=prod.index-1;while(loop>=0){if(*(prod.p+loop)==int2char(0))loop--;elsebreak;}prod.index=loop+1;if(prod.index==0) *(prod.p+Length-1)='0';//商加1*(result.p+quot-1)=int2char(char2int(*(result.p+quot-1))+1);}}if(result.index!=0) *(result.p+Length-1)=sign2c(sign2i(*(p+Length-1))*sign2i(*(itg.p+Length-1)));return result;
}
移位的本质是一种乘法,只是利用了数据结构的特点而得到了极大的便宜。
长除法的优化
仔细分析长除法的直减版本,可以发现以下两点:
- 存储效率太低,每个char至少能表达256个数值,但我们只使用了其中10个,这直接导致了处理能力的浪费。
- 为了确定每一位的数值,我们平均执行了5次减法(假定数值是随机的),但现实中我们在列竖式的时候绝对没有做这么多次减法。
第一个问题似乎容易解决,我们只要在选定的数据类型中使用尽可能大的进制(比如对char我们可以使用百进制)即可。对加、减、乘三种运算而言,这是正确的。但对除法,情况有所不同。
仍以char为例,当进制由10变成100的时候,存储长度可以缩短一半,在其他不变的情况下,时间复杂度的常数可以变为原本的1/4。但是为了确定每一位,平均执行减法的次数从5次变成了50次!实际的常数变化率为1/4*(50/5)=2.5!除法的运行时间反而变长了。
所以不管是为了压位,还是为了进一步减少减法提升效率,我们都需要试商。
1. 试商
如上文提及,我们在列竖式的时候总是能够快速的猜到每一位的最终值(或附近),我们是怎么做到的呢?下面的故事引自华罗庚的《天才与锻炼》:
有一天教授要给学生们出一道计算题。一位助手取来了题目,是一个871位数开97方,要求答案有9位有效数字。教授开始在黑板上抄这个数:456,378,192,765,431,892,634,578,932,246,653,811,594,667,891,992,354……当抄到二百多位后,教授的手已经发酸了。“唉!”他叹了一口气,把举着的手放下甩了一下。这时一位学生噗嗤一声笑了起来,对教授说,当您写出八位数字后,我已把答案算出来了,它是588,415,036。那位助手也跟着笑了,他说,本来后面这些数字是随便写的,它们并不影响答数。这时教授恍然大悟,“哈哈,我常给你们讲有效数字,现在我却把这个概念忘了。”
如引文所述,假如我们把大整数都用科学记数法表示,为了得到结果的最开始一位或几位有效数字,被除数和除数只需要几位有效数字就足够了。那么实际计算中到底需要几位呢?下面的不等式可以给出答案:
图中,X的整数部分即为当前位的最终值,P为进制。
由图可知,如果对被除数进行截断,对除数进行凑整(注:我们希望试商的结果不大于最终值,毕竟减过头就麻烦大了。 ),只要当除数的有效数字大于进制,就可以保证试商结果小于等于最终值,并且试商结果的误差就不大于1。
2. 长除法的试商版本
下面的代码是采用试商法的实现,相比于上面的直减版本,更接近列竖式的笔算。
Integer Integer::operator/(const Integer& itg)const{Integer result(false),prod,div,base;int loop,offset,quot,mult,reg,trans=0;long lp,lb;/*处理除0异常*/prod=*this;if(*(prod.p+Length-1)=='-') *(prod.p+Length-1)='+';div=itg;if(*(div.p+Length-1)=='-') *(div.p+Length-1)='+';while(!(prod<div)){base=div;quot=1;//确定移位的具体值loop=prod.index-1;offset=base.index-1;while(*(prod.p+loop)==*(base.p+offset)){if(offset>0){loop--;offset--;}else break;}if(char2int(*(prod.p+loop))>=char2int(*(base.p+offset)))offset=prod.index-base.index;elseoffset=prod.index-base.index-1;//移位if(offset>0){for(loop=base.index-1;loop>=0;loop--)*(base.p+loop+offset)=*(base.p+loop);for(loop=0;loop<offset;loop++)*(base.p+loop)=int2char(0);base.index=base.index+offset;quot=quot+offset;}if(prod.index>base.index) *(base.p+prod.index-1)=int2char(0);//结果初始化,仅执行一次if(quot>result.index){for(loop=0;loop<quot;loop++)*(result.p+loop)=int2char(0);result.index=quot;}//试商lp=char2int(*(prod.p+prod.index-1));lb=char2int(*(base.p+prod.index-1));for(loop=prod.index-2;loop>=(base.index-2)&&loop>=0;loop--){lp=lp*Radix+char2int(*(prod.p+loop));lb=lb*Radix+char2int(*(base.p+loop));}if(div.index>2) lb=lb+1;mult=static_cast<int>(lp/lb);while(!(prod<base)){if(mult){//执行试商后的减法for(loop=offset;loop<prod.index;loop++){reg=char2int(*(prod.p+loop))-mult*char2int(*(base.p+loop));if(trans){reg=reg-trans;trans=0;}if(reg<0){if(reg%Radix){trans=1-reg/Radix;reg=reg%Radix+Radix;}else{trans=-reg/Radix;reg=0;}}*(prod.p+loop)=int2char(reg);}}else{//执行试商后的修正for(loop=offset;loop<prod.index;loop++){reg=char2int(*(prod.p+loop))-char2int(*(base.p+loop));if(trans){reg=reg-1;trans=0;}if(reg<0){trans=1;reg=reg+Radix;}*(prod.p+loop)=int2char(reg);}}//清除前导0loop=prod.index-1;while(loop>=0){if(*(prod.p+loop)==int2char(0))loop--;elsebreak;}prod.index=loop+1;if(prod.index==0) *(prod.p+Length-1)='0';//商加上试商结果和修正if(mult){*(result.p+quot-1)=int2char(char2int(*(result.p+quot-1))+mult);mult=0;}else *(result.p+quot-1)=int2char(char2int(*(result.p+quot-1))+1);}}if(result.index!=0) *(result.p+Length-1)=sign2c(sign2i(*(p+Length-1))*sign2i(*(itg.p+Length-1)));return result;
}
可以看到,如果配合压位,短除法可以被长除法近乎完全涵盖,两者的效率差可以控制在一个相当有限的范围内。
3. 效率比较
下面的截图是直减版本(1.0.7)和试商版本(1.0.8)的效率比较。测试方法为对梅森数M61使用1亿以内的质数进行取余运算,看是否能整除(当然不可以,因为M61是质数。取这个数的原因是它在64位以内,可以用来和系统自身的除法进行效率比较 )。
可以看到试商法比直减法快了近2.5倍(=32.698/13.166)。
4. 压位
在有了试商以后,char就可以采用百进制,这将进一步提升大整数的效率。但更进一步使用更大的数据类型时就需要注意硬件和编译器的支持问题了。
事实上,在上面的试商版本中,在最差情况下被除数需要3个字节,除数需要2个字节以保证试商结果,故试商时采用long进行除法运算,虽然符合C++语言的规范,但这对于16位的系统并不友好。
另一方面,针对不同的硬件和编译器,int在2字节和4字节之间摇摆,long在4字节和8字节之间徘徊,这都会影响代码的使用(不过对于专业码农,这可能不算个问题吧 )。
其他
关于试商法的论述,还可以参考:
- kedixa的C++大整数运算,里面的除法部分给出了试商的另一个来源。
- 百度百科词条高位试商法,不过笔者认为写得相当一般。
关于大整数类的完整实现,可以从C++大整数类高精度运算库下载,该版本使用short类型并采用万进制,32位(64位)系统中耗时约为long long的6(10)倍。
C++ 大整数运算 高精度除法相关推荐
- 大整数运算(高精度运算)C/C++
前言 这种类型,在做题过程中多为观察所给数据可能造成的大小来选择是否使用. 属于模板类型,学习者理解其格式并记住大致框架即可熟练应用. 一.什么是大整数(高精度) 想知道什么是大整数,不如换一个解释的 ...
- java 大整数编程_Java编程--RSA算法中的大整数运算
Java编程–RSA算法中的大整数运算 RSA原理浅析 RSA是利用陷门单向函数实现的,其安全基础依赖于大整数的分解问题的难解性 算法过程 为了加深对RSA算法的了解,接下来通过简单的一个例子来分析一 ...
- 九度OJ 1037:Powerful Calculator(强大的计算器) (大整数运算)
时间限制:1 秒 内存限制:32 兆 特殊判题:否 提交:1821 解决:528 题目描述: Today, facing the rapid development of business, SJTU ...
- 随手记——大整数运算模板(进化史)
大整数相加.相乘 2019年2月23日09:46:21 基本上没有空间浪费.关键是思路清晰,实现起来方便,字符串倒过来放到vector里(倒过来 方便进位运算),然后做完运算再逆序回来. 关于乘法运算 ...
- C语言 大整数运算(加、减、乘)
题目:大整数计算 背景介绍: 大整数一般指超过十尾的十进制整数,假定不超过五十位.这类大整数在C语言系统中因超界溢出而不能直接表达或计算. 实现方法: 以字符串形式输入.输出和存放大整数,计算时可以将 ...
- matlab将数扩大为整数,MATLAB如何完成大整数运算问题?
Forcal+HugeCalc可以计算下面的数(只有指数部分不能是大整数): (1000!)^1000%(2000!) 代码: !using["HugeCalc"]; mvar: ...
- Go语言的big包实现大整数运算
程序虽然写出来了,但是不知道如何用一个大数(例如100位的大数)去初始化一个大数变量,比较遗憾! Go语言程序: // bigint project main.go package mainimpor ...
- A1136 | 字符串处理、大整数运算
题目链接: https://www.patest.cn/contests/pat-a-practise/1136 今天是12月17号.最近这几天都有点不在状态.已经整整一周没有练算法了,自从12.3考 ...
- 大整数运算之 大整数加法、减法、乘法
其实大整数的问题都是在像我们打草稿的时候列竖式一样的,不要告诉我你不知道什么叫竖式~!其实我开始也不知道它叫这个名字: 所谓竖式,就是你打草稿算算术的方法,小学知识:比如你写 11+9: 11 + ...
最新文章
- sdr 软件_SDR 软件定义的无线电
- 理解java中的两种接口
- 【Spark Summit EU 2016】在在线学习中使用Structured Streaming流数据处理引擎
- 高级C语言教程-关键字和运算符
- Qracle学习:排序
- python四舍五入round_四舍五入就用round( )?Python四舍五入的正确打开方式!
- .NET Core性能测试组件BenchmarkDotNet 支持.NET Framework Mono
- 空白世界地图打印版_考研准考证打印什么时候_中国研究生招生信息网官网
- 金融时间序列计算分析题1
- C语言动态存储分配函数
- 让小黑人360度旋转的制作技巧
- 中国天然石墨行业市场供需与战略研究报告
- 故宫网站遭“围攻”!
- JavaScript 原型精髓 #一篇就够系列
- 在线制作车牌效果图_在线快速生成,苹果设备在线样机
- 清华大学计算机学院2021拟录取,【盛世清北】2020年清华大学(清华)计算机系考研复试拟录取信息...
- 作战管理系统:现代化作战体系核心
- Revit 和 ArchiCAD 在软件设计理念方面的对比
- OA办公系统能帮助企业做些什么?
- 记录一下Material Dialogs的使用