高精度除法之估商法

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

落叶

这个是针对（用户的表达式），用户有可能范点小错误，我是为他们纠正一些允许的错误，方便用户体验，为了这些我可是死了好多脑细胞，而且数学中允许(-8)(-9)这种算式不写乘号，正确的表达式更好识别，
不是我要搞的怪模怪样，只是为了表达我的程序中有纠错功能，具体得数我又算了几遍，和我的结果差不多，当然这段时间一直在改代码，又出现新的错误可能还没发现，不过你那个得数是负数应该不对
-(-3) (-4)(-5)(-6)(-7)(-8)(-9)cuberoot(27)*3.39这个数是里面最大的数，是正数，而且你算的不是高精度数，高精度表达式计算要比非高精度表达式编程复杂太多了，希望你再看看。

上面是错的，sin(5)的得数对了，符号错了，我算成正号，我再改改
弧度化简程序错了，
(-(-1)*(-2)-(-3) *(-4)*(-5)*(-6)*(-7)*(-8)*(-9)*cuberoot(27)*3.39+(27/ 3.39^2-3.39)/3+sin(5))^999   =6.1257445052316170859779980021768321731143396947920769987792051076360748536399302800753283701594905922E6259

补充内容 (2016-3-24 16:06):
补充一下，上面这个是正确答案，上一贴的答案错了，原因是这里面的三角函数是才写好的，因为速度慢，一直没怎么调试，里面那个把大弧度值转换成小弧度值的程序有问题，把应该是负数的答案算成了正数。

mathe

数学里两常数相乘通常不能省略乘号的，只有常数和变量相乘以及当字母表示的变量中相乘才会省略乘号。而计算机中经常会用多字母单词代表变量（为了可读）所以乘号不能省略

mathe

现存高精度数学运算软件很多，要善于利用搜索引擎寻找

落叶

mathe 发表于 2016-3-23 22:09
数学里两常数相乘通常不能省略乘号的，只有常数和变量相乘以及当字母表示的变量中相乘才会省略乘号。而计算 ...

有乘号更好，他也能正常运算，只是为了用户方便，我可以随便把这个功能取消的。另外我当时没找到高精度表达式运算软件，低精度表达式的老外做的很好，
其实这是我的第一个程序，主要目的是熟悉编程操作。最关键是我英文零，搜索引擎寻找的英文东东不明白。
我也知道不要重复学做一个轮子的故事，要站在巨人的肩膀上，主要还是学习。自已能写的就写写，不能写的再调用别人的。
主要是希望得到帮助和交流。其时我写表达式实现和编译程序类似，通过学习写这个软件，感觉原先一直看不懂的编译原理松动了。
三楼的坛友的一翻话提醒了我，其实看网上说除法通过乘法提高，我一直在想是那个方法，原来是迭代法，三角函数太慢了他也为我指出解决方法，我也对他的除法提出
了看法，也希望你们在算法，数学上给予指导，谢谢！

liangbch

非常理解楼上。我们中的许多人也不是一开始就知道这个世界上有那么多功能强大的数值计算和符号计算软件的，都经历了一个从无知到熟悉的过程。另外，拿高精度表达式计算器来练习编程能力是个很好的主意。我想，许多人（包括我和gxq）都这么干过。楼主加油。

只是呼吸

,不过他在采用万进制的情况下,第一位除数采用五位数估商,应该有问题,我曾经有过类似的经历,不过没有成功,也可能他是对的,必竟我没调试过他的程序.

  啊啊，网友落叶可能说的是我吧~
  这个方法是我想了很久后找到的，我自己找到了严格的数学证明，既然数学上是支持的，那就不会有问题的。要出问题只能是程序中的实现过程。
当采用万进制时，除数的第一位要用5个数字，采用10^8进制时，除数的第一位要用9个数字。采用10^9进制时，除数的第一位要用10个数字。采用10^k进制，除数的第一位要用k+1位数字。
这样估出来的商最多相差1，这不是我乱讲的，这是数学证明的结果，如果大家有兴趣，我可以把证明写出来。

别外：我写的程序，基本上是函数级独立的，你只要调用中注意接口就不会有多大问题。仅仅是我写的程序只有初步的测试，对各种情况是不是正确，就不知道了。

liangbch

说说我的方法。

1.术语：
  基，就是进制，下文用R来表示，R可等于 \( 10^4，10^8, 10^9 \) 等.
  机器字：  在32位系统，一个机器字至多可表示 \( 2^{32} \)以内的数, 在64位系统，一个机器字可表示\( 2^{64} \) 以内的数。
  半机器字，在32位系统，一个半机器字至多可表示 \( 2^{16} \)以内的数, 在64位系统，一个半机器字可表示\( 2^{32} \) 以内的数。
  比特： 一个2进制位。

2.提高精度的基本原则
  一般的，对于乘法和除法。如果希望结果精确到n比特，两个源操作数至少包含n比特数据。所以，为了估商准确，我们总是尽可能让更多的比特参与运算。在32位系统，R=10000时，尽量让除数达到半机器字。当 \( R=10^8 \) 或者 \( R=10^9 \)时，尽量让除数达到一个机器字。

3. 大数除法的估商
3.1 定义：
被除数用数组a来表示，a[0]表示最重要的单元， a[1]表示次重要的单元，在通常情况，a中的所有元素都是>=0且<R的数。
除数用数组b来表示，b[0]表示最重要的单元， b[1]表示次重要的单元，在通常情况，b中的所有元素都是>=0且<R的数。
被除数有n个单元，除数有m个单元，n>=m且m>=2, 被除数的前(高）m个单元小于除数(否则可在被除数最高单元前补0)。


3.2 估商，普通的方法是 求 \( (a[0] \times R + a[1] )/ b[0] \)
3.3 更好的估商方法，使用普通的估商法求除数，在某些情况下，可带来较大的误差，所以我们需要将被除数和除数做一些处理，来减少误差。

3.3.1. 当 \( R= 10^4 \)时，当b[0]分别小于等于6,65,655,6553时，我们可以将被除数和除数同时乘以10000,1000,100,10，且最高单元不进位。这样b[0]的比特数将变多了，估商也就更准确了。

3.3.2 当 \( R=10^8 \)时，当b[0]分别小于等于4,42,429,4294，42949，429496，4294967，42949672，99999999 时，我们可以将被除数和除数同时乘以\(10^9，10^8，10^7，10^6，10^5，10^4,1000,100，10 \). 同理，由于b[0]的比特数变多了。估商也更准确。

3.3.3. 当\( R=10^9 \)时，当b[0]分别小于等于4,42,429,4294，42949，429496，4294967，42949672，429496729 时，我们可以将被除数和除数同时乘以\( 10^9，10^8，10^7，10^6，10^5，10^4,1000,100，10 \). 同理,这样使得估商也更准确。

注意，当被除数和除数扩大时，不必增加被除数和除数的长度，可保持a[0]和b[0]大于R，在整个运算过程中，a[0]可始终 >=R

liangbch

楼上的算法仍然不够好，仍可改进。
1.将被除数和除数同时乘以\( 10^k \)是不必要的。
2.乘以\( 10^k \)这种方法，可增加被除数和除数的比特数，但是仍不够彻底。


更好的方法。

1. 考察被除数，若被除数最高的m个单元小于除数，最高位补零。

2. 考察除数的最高2个单元，b[0]和b[1], 若\( [b0] \times R + b[1] \) 小于半机器或者1个机器字, 那么做如下预处理。
  2.1. 将b[0]和b[1]合并，做 \( b[1] = b[0] \times R + b[1] \),同时将b[0]清0.
  2.2. 对被除数也做同样的处理。
这样，将 n个单元的被除数 除以 m个单元的除数 的 除法 转化为 n-1个单元的被除数 除以 m-1个单元除数 的 除法.这个变换仅需做一次。

3.考察除数的最高单元b[0]，若其比特数小于半个机器字或者1个机器字，扩大\( 2^k \)倍,将其变半个机器字或者1个机器字,存入 \( bh \)，这个步骤也仅需做一次。

4.试商，并依次得到各个部分商，在这个过程中，被除数要扩大同样的倍数，这个过程要重复许多次。


下面以\( R=10^4 \)为例，详细解释下。

  例1.被除数为(5,9999,1234,5678),除数为(6,2345,5432)。

  step0，预处理：除数的前两个单元表示的数6*R+2345=62345小于半个机器字65536, 将被除数和除数的前2个单元合并。合并后
    被除数变为(59999,1234,5678),除数为(62345,5432)
   
  step 1: 除数最高单元为62345，包含16比特数据，无需将其扩大。

  step 2.试商，用被除数前2个单元除以除数的第一个单元。
         得部分商 \( 9623=(59999 \times 10000+1234) \div 62345 \)
  
  step 3.部分商9623乘以除数(62345,5432),得到部分积(59995,1162,2136)
  
  step 4,被除数减去部分积，得余数(40072,3542)

  接着，可重复step2,step3,step4，以得到更多的部分商
  

注：从这个例子中可以看出，在计算过程中，被除数，部分积和余数的最高位是可以大于R的


  例2.被除数为(599,9912,3456),除数为(623,4555)。

  step 0，除数的前2个单元表示的数 \( 6234555 > 2^{16} \),不需要做预处理。
   
  step 1. b[0]=623,只有10比特，\( 2^9 < 623 < 2^{10} \) ,乘以\( 2^6 \) 扩大至16比特. 即 \( bh= 623 \times 64 + (4555 \times 64 \div 10000)=39901 \)

  step 2. 试商，用被除数前2个单元扩大64位，然后除以除数的第一个单元,得部分商9623
       \( 9623=(599 \times 10000+9912)\times 64 \div 39901 \)

  step 3.部分商9623乘以除数(623,4555),得到部分积(599，9512，2765)
  
  step 4,被除数减去部分积，得余数(400,0691)
  
  重复 step2,step3,step4得到更多的部分商


这个估商法的优缺点。
  1.相对于普通的估商法，这个方法得到的商更加准确，部分商偏小的概率大大降低。
  2.相对于楼上的估商法，这个方法不需要将被除数和整数整个乘以 \( 10^k \)
  2.相对于普通估商法，在估商时，需要做一个额外的运算步骤，被除数要扩大\( 2^k \)倍，不过，这个操作可用左移指令来实现，增加的成本可忽略不计。

一点儿说明，在实践中，为了避免估商偏大，\( bh \)要稍微取得大一点儿。在例2中。 计算\( bh \)时，其步骤改为 \( bh= (623 \times 64 + 4555 \times 64 \div 10000)+1=39902 \)   

liangbch

关于基的选择。
我认为，
  1.基的选择不宜太小，太小则不能充分发挥CPU计算能力。
  2.基的选择不宜太大，基太大时，估商的误差会比较大，另外，计算乘积的累加很容易溢出。

我的建议是。在32位系统，基可选择 \( 2^{30} \) 或者 \( 10 ^9 \), 这样
1。估商时，最大误差不超过1
2.  做乘积的累加和时，双精度结果可承受4次乘积的累加和而不溢出。

相应的，在64位系统，我建议选择\( 2^{60} \) 或者 \( 10^{18} \) 作为基。
当基为\( 2^{60} \) 时，双精度结果可承受 256次乘积的累加和而不溢出。
当基为\( 10^{18} \) 时，双精度结果可承受 340次乘积的累加和而不溢出。