高精度除法之估商法

落叶

网友 （只是呼吸），你的理论是对的，当时我模拟你的程序时，估出来的商是三位数（被除数8位，除数5位只能估出三位商），我个人认为万进制时估出四位数才能效率更高，所以我就没向下模拟了，昨天我又推算了一下，你是对的，我曾经看过和你的思路类似的思路，用被除数的前三个进制位（这里假设万进制时：一个进制位为四位）除以除数的前两位估商，估商精度要高于Knuth(高纳德)的算法（最大误差2，平均1），不过32位机存不下三个进制位，要单独处理，我放弃了，看了你的回贴，我又重新思考这个问题，得出如下结论：估商时并不一定要用三个进制位，估四位商时，参入估商的除数位大于四位时，可以大大增加正确率，但是你的程序,如果把被除数参加估商的位数设为9位时（而不是8），可以得出四位商，下面这个是我的估商源码：
方法一：ysjg(i) = yShu(2) \ b(2)        ‘ysjg(i)存放估出的商， yShu(2)存放被除数的前8位， b(2) 存放优化过的除数前四位，估商误差50%左右，应该小于50%，只是粗略估计，没有数学论证（关键是我不会）
这个是我借鉴上述思路后改进的代码：
方法二：ysjg(i) = (yShu(2) * 10 + yShu(3) \ 1000) \ (b(2) * 10 + b(3) \ 1000)，现在估商误差5.2%，每计算出2500（相当于一万位）个万进制得数，估商错误修正函数运行130次左右（余数加一修正）
下面是我的一组测试数据：
被除数一万位，精度一万位：里面有2500次估商或求商操作：一般来说，第一个进制位越小,第二个进制位越大，对估商影响越大，反之估商越准确。
除数                                 方法一修正次数       方法二修正次数                                   除数                                 方法一修正次数       方法二修正次数                     除数                  方法二修正次数
100099999999               1425                         162                                                     111111111111                1256                        111                                       987654321111      52
200099999999               1563                         148                                                     222222222222                1267                        125                                       100000000000      0
300099999999               1346                         128                                                     333333333333                1284                        126                                       200010001000      0
400099999999               1562                         149                                                     444444444444                1206                        129                                       444422222222      70
500099999999               2520                         270                                                     555555555555                1219                        113                                       666622221111      40
600099999999               2064                         180                                                     666666666666                1263                        119                                       700000000000      0
700099999999               1811                         172                                                     777777777777                1303                        131                                       777744444444      73
800099999999               1557                         142                                                     888888888888                1257                        118                                       845612345678      32
900099999999               1435                         150                                                     999999999999                1251                        113                                       555522222222      50
现在估商误差已经很小，甚至可以忽略不计.,感谢你的回贴,思路!
看了你的代码,感觉你的cpu硬件,汇编学的好,编程的控制竟然和CPU控制寄存器类似.
看了宝宝的代码,感觉数学超强,汇编,硬件方面超强. ,感谢宝宝的奉献！

补充内容 (2016-4-22 10:01):
估商的误差主要是除数第二个进制位的进位造成的，当第二个进制位全部为零时，估商准确率接近百分之百，例：22220000，因为是进位造成的误差，所以我们用整除方式估商时，错误的商只会大于正确商

落叶

下面是我的估商除法源代码:
Public Function myEXCEPT(str1 As StrToZx, str2 As StrToZx) As StrToZx        '除法子程序
        Dim a() As Long
        Dim b() As Long
        Dim i As Long
        Dim ysjg() As Long                                            '存放商
        Dim ysjg4to1() As Long                                        '把万进制商转换为十进制商时，存放在这儿。
        Dim Bjs() As Long                                             '存放减数，后面的核心运算有大数减法
        Dim yShu() As Long                                            '存放运行中的余数
        Dim jg As StrToZx                                             '存放最后结果
        Dim string1 As StrToZx
        Dim string2 As StrToZx

        string1 = str1
        string2 = str2
        Dim cszd As Integer
        cszd = 1
        string1 = Zjzzh(string1, cszd)                                                 ' 转换成10进制数
        string2 = Zjzzh(string2, cszd)
        string1 = jdkz(string1, jd + 4)                                                精度控制,减少不必要的运算
        string2 = jdkz(string2, jd + 4)
        string2jw1 = string2.XsdWz - 1                                                  '解决除数0.000...01的问题  如不调整0.001也可以运行
        string2.XsdWz = 1
        
        string1 = qslqwl(string1)                                                    '除去数首和数尾部的零
        string2 = qslqwl(string2)
        For i = 1 To string2.strlen                                                  '除数为零提示返回
                If string2.Zx(i) <> 0 Then
                    Exit For
                ElseIf i = string2.strlen Then
                    MsgBox "除法子程序中除数为零", vbExclamation, "错误提示"
                    Exit Function
                End If
        Next

        cszd = 4
        string1 = Zjzzh(string1, cszd)                            '转换为万进制        
        string1 = qslqwl(string1)                                  '除去数首和数尾部的零
         
        Dim tempQszl As Long                             '获取除数最佳倍乘数(qszl ),下面一大段是关键步骤,运行时间可忽略不计
        Dim tempSz(-1 To 10) As Byte                     '清零
        For i = -1 To 10
            tempSz(i) = 0
        Next
        tempQszl1 = string2.strlen                       '后面几步是判断应该获取的除数长度，大于8，取8位，小于8，取实际长度
        If tempQszl1 > 8 Then
           tempQszl1 = 8
        End If
        For i = 1 To tempQszl1
            tempSz(i) = string2.Zx(i)
        Next
        tempQszl = 0
        For i = 0 To 7
            tempQszl = tempQszl + tempSz(8 - i) * 10 ^ i
        Next
        qszl = 0
        For i = 1 To 10
            temp = tempQszl * i
            If temp >= 100000000 Then
                qszl = i - 1                             '除数整理，qszl 参数乘以除数后，使除数的第一个十进制位接近10,例:1234中的1变为8
                Exit For
            End If
        Next


        str2len = string2.strlen - string2.XsdWz                               这一段就是把例:0012 3456转换成1234 5600,运行时间忽略不计,下面一段是关键步骤
        LL = str2len Mod 4
        string2jw = 0
        If LL <> 0 Then
            string2jw = 4 - LL
            string2.strlen = string2.strlen + string2jw
            ReDim Preserve string2.Zx(-1 To string2.strlen + 2) As Long
        End If
                                                                    
         JzZh1to4 string2                                                               '十进制转换成万进制

        If qszl > 1 Then                                                                  '下面这段代码会占用一定的计算时间(单精对多精乘两次),下面一段是是关键步骤
                jw4 = 0                                                                        '除数乘（调整参数qszl）
                For i = string2.strlen To 1 Step -1
                    string2.Zx(i) = string2.Zx(i) * qszl + jw4
                    jw4 = string2.Zx(i) \ 10000
                    string2.Zx(i) = string2.Zx(i) Mod 10000
                Next              
                jw5 = 0                                                                  '被除数乘（调整参数qszl）
                For i = string1.strlen To 0 Step -1
                    string1.Zx(i) = string1.Zx(i) * qszl + jw5
                    jw5 = string1.Zx(i) \ 10000
                    string1.Zx(i) = string1.Zx(i) Mod 10000
                Next
                If string1.Zx(0) <> 0 Then                                               '如果产生了进位，把被除数整体后移一位
                    For i = string1.strlen To 1 Step -1
                        string1.Zx(i) = string1.Zx(i - 1)
                    Next
                    string1.strlen = string1.strlen + 1
                    ReDim Preserve string1.Zx(-1 To string1.strlen + 2) As Long
                End If
                string1.Zx(0) = 0               
        End If


        string2 = qslqwl(string2)
        GzBzsLen = jd \ 4 + string2.strlen + 4                 '扩展被除数长度
        a() = string1.Zx()
        b() = string2.Zx()
        tempstr2 = string2.strlen + 2
        ReDim Preserve a(-1 To GzBzsLen + 2) As Long             '存放扩展后的被除数
        ReDim Preserve b(-1 To tempstr2) As Long                 '存放除数
        ReDim ysjg(-1 To GzBzsLen + 2) As Long                   '存放商
        ReDim ysjg4to1(-1 To GzBzsLen * 4 + 2) As Long           '存放十进制商
        ReDim Bjs(-1 To tempstr2)                                '存放减数
        ReDim yShu(-1 To tempstr2)                               '存放余数
        For i = string1.strlen To 1 Step -1                 '左前方预留一位零 ,                  是关键步骤
            a(i + 1) = a(i)
        Next
        a(1) = 0
        string1.strlen = string1.strlen + 1                  '左前方预留一位零 ,                是关键步骤
        For i = string2.strlen To 1 Step -1
            b(i + 1) = b(i)
        Next
        b(1) = 0
        string2.strlen = string2.strlen + 1


        azd = string2.strlen                                 '被除数是取出和除数等长的数到余数参加运算，这时这个标志的所标位置，即是被除数待取数位置标志。也作为退出运算标志。是关键步骤
        For i = 1 To string2.strlen                          '从被除数中取出和除数等长的数，被除数和除数左方必须加一个零例：987/12 转换成 0987/012     
            yShu(i) = a(i)
        Next

        i = 1                                                                  下面是核心运算代码
        Do
line2:
            If yShu(2) > b(2) Then                                   '当余数大于除数时的处理
                ysjg(i) = (yShu(2) * 10 + yShu(3) \ 1000) \ (b(2) * 10 + b(3) \ 1000)        '估商            是关键步骤
line1:
                jw1 = 0
                For j = string2.strlen To 1 Step -1                        '单精对多精乘
                    Bjs(j) = ysjg(i) * b(j) + jw1
                    jw1 = Bjs(j) \ 10000
                    Bjs(j) = Bjs(j) Mod 10000
                Next
                Bjs(2) = Bjs(1) * 10000 + Bjs(2)
                Bjs(1) = 0
               
                jw2 = 0                                                    '多精对多精减
                For u = string2.strlen To 1 Step -1
                    yShu(u) = yShu(u) - Bjs(u) - jw2
                    If yShu(u) < 0 Then
                        yShu(u) = yShu(u) + 10000
                        jw2 = 1
                    Else
                        jw2 = 0
                    End If
                Next
               
           '商修正，余数加一，估出的商只会大于或等于正确商,当商不正确时，做减法时会从 yShu(1)借位， yShu(1)从零变为9999，修正后 yShu(1)变回零,   是关键步骤,重点
                jw3 = 0         
                Do While yShu(1) = 9999
                    For u = string2.strlen To 1 Step -1
                        yShu(u) = yShu(u) + b(u) + jw3
                        If yShu(u) >= 10000 Then
                           jw3 = 1
                           yShu(u) = yShu(u) - 10000
                        Else
                           jw3 = 0
                        End If
                    Next
                    ysjg(i) = ysjg(i) - 1                               '商修正
                Loop
                yShu(1) = 0
            End If

            If yShu(2) < b(2) Then                                 '当余数小于除数时的处理
line3:
                yShu(2) = yShu(2) * 10000 + yShu(3)
                For v = 3 To string2.strlen - 1
                    yShu(v) = yShu(v + 1)
                Next
                If azd = GzBzsLen Then                             '这里控制运算长度，从这里退出循环
                    Exit Do
                End If
                azd = azd + 1
                yShu(string2.strlen) = a(azd)
                i = i + 1
                GoTo line2
            End If

            If yShu(2) = b(2) Then                                 '当余数和除数相等时的处理
                For w = 3 To string2.strlen
                    If yShu(w) < b(w) Then
                        GoTo line3
                    End If
                    If yShu(w) > b(w) Then
                        ysjg(i) = 1
                        GoTo line1
                    End If
                    If w = string2.strlen Then
                        ysjg(i) = 1
                        GoTo line1
                    End If
                Next
            End If
        Loop
        
        bz = 0                                                    '转换成十进制，因为小数点的对位要在十进制下进行
        j = (jd \ 4) * 4 + 16
        For i = (jd \ 4 + 4) To 1 Step -1                                                           '
            ysjg4to1(j) = ysjg(i)
            Do
                ysjg4to1(j - 1) = ysjg4to1(j) \ 10
                ysjg4to1(j) = ysjg4to1(j) Mod 10
                j = j - 1
                bz = bz + 1
            Loop While bz <> 4
            bz = 0
        Next

         GzBzsLen4 = GzBzsLen * 4                                                     '下面几句代码主要是确定小数点位置，数字长度，正负,指数
         string1.strlen = string1.strlen * 4 - string1.XsdWz + string2.XsdWz
         jg.JzBz = 1
         jg.strlen = GzBzsLen4
         jg.XsdWz = GzBzsLen4 - string1.strlen + string2.strlen * 4 - string2jw - 4 - string2jw1
         jg.Zx() = ysjg4to1()
         temp1 = jg.XsdWz - jg.strlen + 1
         jg.eE = jg.eE - temp1                                      
         jg.eE = jg.eE + string1.eE - string2.eE                          '计算最终结果的指数
         jg.XsdWz = jg.strlen - 1
         If string1.ZhFhBz <> string2.ZhFhBz Then                    '根据两运算数的正负符号确定最终得数的正负符号
            jg.ZhFhBz = False
         Else
            jg.ZhFhBz = True
         End If
      
        cszd = 1
        jg = Zjzzh(jg, cszd)                                           转换成十进制数
        jg = jdkz(jg, (jd + 4))                                        '精度控制
        jg = qslqwl(jg)                                                '除去数的首零和尾零
        myEXCEPT = jg
    End Function


上面的程序并没在代码易读性上优化,某些方面可能重复,也可能误删了代码,(源代码中插有一些调用迭代除法的代码,测试代码,为了不影响阅读,作了删除操作)不明白的再问我.
程序支持整数,小数,采用科学记数法的数的运算,并通过测试,并非凭空猜测.

只是呼吸

估出来的商是三位数（被除数8位，除数5位只能估出三位商），

我们来探讨一下，请（落叶）不要有其他想法，（被除数8位，除数5位只能估出三位商），实际上，这个估出来的商应该是四位商，一个例子：8341 8341/12345=6765 （小数部分省略），这个例子已经是四位商，这不用多说。
我们再看第二个例子  8341 8341/86321=966（小数部分省略）。这个商看上去是三位，其实不然，这个966仍然是四位商，它的真实值是0966，当这个0966是第一位商时（最高位商），它前面的0就会多余，会被省略，成了966。当这个0966不是第一位商时，它一定是0966，这个0不能省了。这也许会和你的想法不一致，这种不一致会干扰你写代码的质量。


另外：我的汇编是学习自宝宝的，我基本是原样照搬。
再有：看了宝宝和你的贴，我还得重新思考一下除法。我说的那个估商的证明，在这几天中，我会写出来，这样大家用起来信心会大增。

落叶

只是呼吸 发表于 2016-4-21 00:15
我们来探讨一下，请（落叶）不要有其他想法，（被除数8位，除数5位只能估出三位商），实际上，这个估出来 ...

真没想到会出现四位的商，不好意思！ 其实我的代码也是会出现四位商，和其它位的商，不足四位的也是前面加零。
你的思路是对的，我已采用，和我原先的方法叠加，获得了很好的估商效果。