高精度除法之估商法

落叶

我用VB完整的实现了估商法，支持正负，小数，科学计数法，之所以说完整实现，是因为这个程序是我正在写的（高精度表达式计算机）的子程序，能够和其它计算程序配合，实现高精度加减乘除，阶乘，乘方，开N次方，函数等混合运算.
例：能够运行万位精度，开10000次方根正确运算，开N次方根是用的迭代法，为什么举这个例子，因为参入运算的子程序如果有问题，是计算不出正确结果的，现在这个计算机还没完善:
先说优点，国内目前找到的表达式计算机(不支持高精度)，而且不完善，几乎无纠错功能,最关键是不支持函数,函数是一切功能运算的基础,这些问题在我的程序里已经建立了解决框架.可以迅速加入纠错功能,各种函数.
再说缺点,开N次方没找到好的算法,我的算法是通过利用手算公式取初值,再迭代,     三角函数是用的泰勒公式,速度也是超慢,阶乘针对硬乘优化,10000!3秒左右,准备调用别人的大数库,
乘法是硬乘和FFT乘法分工,具体如何分工没有好的方法,目前采用的是1100位*1100位以内的数用硬乘,和任一单数400位之内用硬乘,
乘方采用网上找到的快速算法,不知道名字,就是A的8次方转换成A的平方的平方,然后再平方,速度非常快.很满意.
说了这么多主要是希望坛主和其他大大们能给出好的算法和建议,
另外我再谈谈我对除法的一些看法,不能利用计算硬件的算法不是好算法,所以设计算法时要充分考虑硬件的效能,简单说一下我知道的算法,第一,用减法实现,从高位向低位减,我最早用的是这个方法并实现,    第二个方法是用的估商(不是试商,试商非常慢),估出一个不大于正确商的数,做乘法,再按减法算法继续,一直觉得这个算法不错,在位数少时,速度还满意,
第三个方法是一个论文上说的方法,用的是十进制,,提前算出除数和1-9相乘的得数,其它步骤和方法二差不多,低位时速度比第二个方法稍慢,高位时甩开了我,另外这个论文中说到了迭代法,说很慢,
第四个方法是迭代法,我实现后就没管它了,和上面说的一致,超慢!
第五个方法就是我现在用的方法,为什么我一直没实现第三种方法,就是因为它采用的是十进制,不能充分利用CPU的乘法功能,需要说明的一点是现在乘法指令和加法一样快,你不相信?,自已写个代码去试试,所以不要再想一些用加减代乘法,或除法的算法,整数除法指令也只比乘法慢一点,所以为了提高性能,32位机上要用万进制乘,用万进制的难点就是商更加难估了,解决要点就是通过找一个参数使除数接近10000,当然通用除法的调试是相当麻烦的,所以上网上你只能看到人讲原理,却看不到实现,其实我们论坛看到一个人用C写估商法,因为他的程序运行困难,我没有调试,不过他在采用万进制的情况下,第一位除数采用五位数估商,应该有问题,我曾经有过类似的经历,不过没有成功,也可能他是对的,必竟我没调试过他的程序.
希望坛友们给予指导,建议!


补充内容 (2016-3-24 14:36):
下面有我回别人的贴子，在里面有对估商法较详细的看法，有不同意见，看法发出来一起计论！

liangbch

楼主，加油

Ickiverar

嗯。我的大数库实现除法，主要就用浮点数除法（牛顿迭代法）和估商两种算法。在位数比较小的时候（几千位之内），估商比较快，再往上就是牛顿迭代比较快了，你觉得慢肯定是乘法没实现好（实现好后是O（n ln n ）的），或者控制精度没做到位。
说估商。我现在用的是2^64进制，计算机硬件支持128bit整数除以64bit整数，只要商不大于64位即可。于是算P/Q的时候，我把Q的最高64位(记作 H )拿出来，加上1，除2^127。得到 R=floor(2^127/(H+1))，满足 2^63<=R<2^64。用这个R去乘P的高位，就能得到商的高位。具体的位对齐需要仔细算。实现后我做了实验，这样估商，一次能估出64个二进制位（19.3个十进制位），估的值总是小于实际值，且差不超过5.

乘方那个算法名称是“快速幂”。

三角函数/幂函数/对数函数，用泰勒公式是没问题的，有加速方法。binary spliting是一种通用方法，还有各种公式也能用来加速，比如sin(a+b)=sinacosb+sinbcosa。这方面中文的资料比较少。

给个建议，能用C/C++就尽量用……如果考虑到效率的话。而且指针是非常强大的工具。换语言很耗时间，所以只是建议。

liangbch

我正在设计一个全新的高精度对数算法，期望其性能达到一个新的高度，甚至超过目前这个星球上最快的软件或者库。请静候佳音，到时我会公布我的结果的。

落叶

前天用手机回贴，乱码
今天才回，很感谢三楼给予的详细解答，我的FFT乘法也是前几天才完成，一直想自已写FFT转换，一直没搞懂细节，没办法，正好网上有C源码，修改了一下做成DLL，前几天才搞好，
因为原先的测试偏见，没想起从新测试一下迭代除法。
仔细看了一下你的估商实现，你把除数直接提出来的方法，相当于左移操作，算完结果后你应该还位（如改小数点位置，改指数，右移等等），我还是用万进制举个例子比较直观：
例：123456换成万进制后的形式是0012  3456，乘100后变成1234 5600，而你估商时是直接用1234，相当于除数乘以100，做完结果后得数要乘100，当然你也可以用其它方法，
对于你用除数加一估商，我在方法二时，也是这样用的，不过我用的是十进制，被除法的前两位除以除数的前一位加一，但我当时的想法是加在除数的最高位，而你是加在最低位，
虽然十进制中，我所取数的最高位就是最低位。
Knuth(高德纳)对估商做出了一个数学推导，公式一大堆，一句话，就是左移除数，乘2乘2...，例：0011 0101 左移成1101 0100,当时硬件差，他以左移替代乘法，估商最大误差2，
平均误差1，这种算法最佳实现是一次单精对多精乘法，一次多精减法，一次多精加法修正（我是直接加余数），书上说的竟然是用积减除数修正，哎，不说了，代沟啊！
  其实不管是你的方法，还是老K的，还有我的估商，总的思想就是整理除数，使除数尽量接近你所采用的进制的最大值，2^32进制时，就是65535，万进制就是10000，是不是不好理解，
举例：万进制时：下面的例子中因被除数第一位小于除数第一位，合并被除数第一位和第二位，整除除数第一位估商，如下：
0999 5678 1000/1000 9999 ,估商9995，正确商为9985,把被除数，除数乘9，8996 1102 9000/9008 9991,估商9986，另外不要问我被除数变了有没有影响，几乎没影响，经测试除数整理
好后，估商误差0.3左右，不超过0.5，我也是不明白，老K的主思路和我一样，为什么误差要大些，最后发现原因就在于左移，只有乘2，乘4，乘8，除数没有调整到最佳状态，为什么这么说呢？
如果一个除数需要调整，乘以3最合适，按老K的算法你就只能乘2了，或乘4了。我觉得老K的左移操作现在对高精度数没多大用处，高精度数本身是有数据头的，向这种左移右移的操作，
在数据头里把指数改一下就可以了，如我的万进制：1234*0100，只需要把123的指数位加2就行了
再说我的算法：
0999 5678 1000/0010 0099，除数乘100（不是真乘，只是调整数的位长度，和小数点位置），jw=2（这个用于后面调整小数位），变成0999 5678 1000/1000 9900，用（1000 9900）*1-9
反复试，找出一个最大的数，使它俩的积不超过100000000，这里是9，把被除数，除数乘9，好吧，进入主程序。
我的所有代码除了总控程序，其它都准备公开源码，没有很高的效率，就是比较完善，调试通过，支持正负，小数，科学计数。
但是现在没注释好，最关键是不好运行，因为没总控程序，只有等把总控做成DLL再说吧，当然如果你们不在意，现在就可贴出来。

关于你的程序我有一点小的看法：你用的进制是满进制，也就是2^64*2^64后,寄存器已经满了，后面的进位你不能一起处理，单独处理很耗时，你可能觉得满进制乘2,4,8,16不耗时，其实高精度数本身
是有数据头的，向这种左移右移的操作，在数据头里把指数改一下就可以了，如我的万进制：1234*0100，只需要把1234的指数位加2就行了。另外你的除数加1 ，你试过没加和加过后的区别吗？

关于三楼坛友对我的建议，我很感谢，其时我最早看过c,c++,mfc,汇编，还有很多很杂的计算机方面的书，就是没怎么看VB，当然因为自已的学习方法错误，一直是动眼不动手，所以都只是混了个脸熟，
知道个大概意思，已经N多年没有看计算机方面的书了，这次是下决心通过写个程序整合以往的知识，并增补新知识，因为英语0，又因VB的中文参考较好，心里其实已是准备两个语言都看看，我只是觉得语言
就是个工具，大同小异，大部分都相通。MFC虽不同，真心觉得不好用，但是你要想看类的应用，简直太棒了，请看它吧？说个笑话：你是喜欢开5档的车，还是喜欢开一百档的车，编程不是编程语言，只是
它们联系很紧而已。

落叶

再来一贴，我刚才忘了问，为什么国内几乎看不到像样的表达式计算机，高精度更是没有，难道那些工程师，科学家，计算数据是算一项记一项，或M存等，真心觉得麻烦，
有熟知这方面的坛友回一下。

wayne

落叶 发表于 2016-3-23 13:55
再来一贴，我刚才忘了问，为什么国内几乎看不到像样的表达式计算机，高精度更是没有，难道那些工程师，科学 ...

我觉得到处都是车轮子,您怎会找不到软件呢, 能否具体说一下你的需求?

落叶

我说的是这种类型的表达式计算:
-(-1)(-2)-(-3)(-4)(-5) (-6)(-7)
(-8)(-9)cuberoot(27)*3.39
+(27/3.39^2-3.39)/3  =1.8452424531466833737959119743128E6

落叶

不好意思，没说清：这个是一百精度的，目前我设计的是一万精度，各种函数可以方便加上：下面是我的计算器运行效果，
(-(-1)(-2)-(-3) (-4)(-5)(-6)(-7)(-8)(-9)cuberoot(27)*3.39+(27/ 3.39^2-3.39)/3+sin(5))^999   =6.1321082202579006232825733921632263575771187476670971901178450674213340042091847280844090084002659067E6259
如果你需要处理比较复杂的算式，你怎么办，一个一个算吗？如果是高精度数呢？
比较先进的计算器是这样的吗？，x=****,y=***,然后x+Y=***.还有其它类似编程的功能，虽然工程师，科学家都很聪明，但是他们都会使用里面类似于程序员编程的功能吗？
这个程序可以加我们日常省略的运算符，去空格，去回车，你可能说这不难呀，上述做的意思，就是我已经建立了一个处理语法错误的程序框架，实现函数
的意义是：大多的运算可以通过函数表达，如果为了直观可以做一个格式转换函数，把一些公式等转换成内在的函数，必竟计算机并不认识公式，但可以
方便的识别函数。

具体你问我的需求，我没需求，只是我开始学写程序时（去年11月份左右，之前只看过书，几乎没动过手），准备写一个计算机程序练练手时，当时突发奇想做一个别人没做过的，
因为我用网上的计算器感到很不习惯，觉得要是能一起算多好啊，分开计算太麻烦，计算出一个得数还要想办法存，当时查到一个贴子，日本的一个手持计算器不能识别括号，
我想也许是手持计算器不好编程吧，然后上网找电脑版计算器，不过也是很失望，最后在excel的功能里面找到了表达式混合运算，并支持函数，但不支持高精度，
然后我就开始构思这个计算器，并加上了多精度功能。

liangbch

"(-(-1)(-2)-(-3) (-4)(-5)(-6)(-7)(-8)(-9)cuberoot(27)*3.39+(27/  3.39^2-3.39)/3+sin(5))^999   =6.1321082202579006232825733921632263575771187476670971901178450674213340042091847280844090084002659067E6259"

你上面的这个表达式语法很奇怪，不太符合一般的语法规则。
一般说来，计算机中的表达式和通常我们数学中的表达式还是有些区别的。在计算机表达式中，乘号是不能省略的。将上面的表达式补上乘号，并将cuberoot(27)，替换为sqrtn(27,3)。在PARI/GP可以很容易算出来。
其结果为-3.288935247083505843230773751 E6560