algorithm - Lehmer 的扩展 GCD 算法实现

Question

在执行我自己的 BigInteger 实现时，我遇到了扩展 GCD 算法，该算法是查找模乘逆的基础。由于众所周知的欧几里德方法执行速度太慢，而混合和二进制算法仅快 5-10 倍，因此选择了 Lehmer 对经典算法的修改。但困难在于，在描述 Lehmer 时，我找到的所有书籍(Knuth、Handbook of Applied Cryptography、Internet 等)都有相同的缺点:

1. 解释基于几个技巧:
   • 输入的数字总是相同的长度；
   • 抽象CPU有signed寄存器，可以同时保存数字和符号；
   • 抽象 CPU 具有半无限寄存器，即。 e.它永远不会溢出。
2. 仅提供基本的GCD算法，不关注逆余因子。

至于第一个问题，我最初对无法找到任何真实世界的实现感到惊讶(不要让我指向 GNU MP 库——它不是学习的来源)，但最终通过反编译获得了灵感Microsoft 从 .Net 4.0 开始的实现，这显然是基于论文“A double-digit Lehmer-Euclid algorithm for finding the GCD of long integers”中的思想。 ” 杰贝林。生成的函数很大，看起来很吓人，但效果很好。

但是微软的库只提供了基本的功能，没有计算辅助因子。好吧，准确地说，一些辅助因子是在速记步骤中计算的，并且在第一步中这些辅助因子只是初始的，但是在执行速记步骤之后然后他们不再匹配了。我目前的解决方案是将“真实”辅助因子与“替代”辅助因子并行更新(第一步除外)，但这会使性能降到零以下:函数现在完成速度仅比二进制快 25-50%基本模式下的方法。因此，问题在于，虽然输入数字仅在速记步骤中完全更新，但辅助因子也在每个速记步骤的迭代中更新，从而几乎破坏了 Lehmer 方法的任何好处。

为了稍微加快速度，我实现了一个“融合乘加”函数，因为“融合乘-乘-减”确实有助于更新输入数字，但这次的影响可以忽略不计。另一个改进是基于这样一个事实，即通常只需要一个辅助因子，所以另一个可以根本不计算。这应该会使开销减半(甚至更多，因为第二个数字通常比第一个数字小得多)，但实际上开销仅减少了预期的 25% 到 50%。

因此，我的问题归结为:

1. 是否有关于 Lehmer 算法的全面解释，与现实世界硬件上的实际实现相关(使用有限大小的无符号字)？
2. 同上，但关于扩展 GCD 计算。

因此，尽管我对基本算法的性能感到满意，但扩展操作模式却恰恰相反，这在我的案例中是主要的。

Answer 1

最后，我咨询了一位数学家，他很快找出了正确的公式——与我自己尝试的公式非常相似，但还是略有不同。这允许在输入数字完全更新的同时，仅在手写步骤上更新辅助因子。

然而，令我大吃一惊的是，仅此一项措施对性能的影响很小。只有当我将其重新实现为“融合 (A×X + B×Y)”时，速度提升才变得明显。在计算两个辅助因子时，与纯 Lehmer GCD 算法相比，它现在对于 5 位数字的运行率为 56%，对于 32K 位数字的运行率为 34%；对于单个辅助因子，比率分别为 70% 和 52%。在以前的实现中，两个辅助因子的比例仅为 33% 到 7%，单个辅助因子的比例为 47% 到 14%，所以我的不满是显而易见的。

至于像andy256推荐的那样写一篇论文，让其他实现者不至于遇到同样的麻烦，这可不是件容易的事。在向数学家解释我当前的实现时，我已经写了一篇“小”论文，它很快就超过了三张 A4 纸——同时只包含基本思想，没有详细解释、溢出检查、分支、展开等。现在我部分理解为什么 Knuth 和其他人使用卑鄙的把戏来缩短故事。目前，我不知道如何在不失去彻底性的情况下达到同样的简单程度；也许，它需要结合评论的几个大流程图。

---

更新。看起来我将无法在不久的将来完成和发布该库(仍然没有运气理解牛顿-拉夫森除法和蒙哥马利减少)，所以我将简单地发布结果函数给那些有兴趣的人。

该代码不包括明显的函数，如 ComputeGCD_Euclid 和通用例程，如 ComputeDivisionLonghand。该代码也没有任何注释(除了少数异常(exception))——如果您想理解它并将其集成到您自己的库中，您应该已经大致熟悉 Lehmer 的想法和上面提到的两位数速记技术。

我的图书馆中数字表示的概述。数字是 32 位无符号整数，因此在需要时可以使用 64 位无符号和有符号算术。数字从最低位到最高位 (LSB) 存储在普通数组 (ValueDigits) 中，实际大小显式存储 (ValueLength)，即。 e.函数尝试预测结果大小，但不优化计算后的内存消耗。对象是value 类型(.Net 中的struct)，但它们引用 数字数组；因此，对象是不变的，i。 e. a = a + 1 创建一个新对象而不是改变现有对象。

Public Shared Function ComputeGCD(ByVal uLeft As BigUInteger, ByVal uRight As BigUInteger,
        ByRef uLeftInverse As BigUInteger, ByRef uRightInverse As BigUInteger, ByVal fComputeLeftInverse As Boolean, ByVal fComputeRightInverse As Boolean) As BigUInteger

    Dim fSwap As Boolean = False
    Select Case uLeft.CompareTo(uRight)
        Case 0
            uLeftInverse = Instance.Zero : uRightInverse = Instance.One : Return uRight
        Case Is < 0
            fSwap = fComputeLeftInverse : fComputeLeftInverse = fComputeRightInverse : fComputeRightInverse = fSwap
            fSwap = True : Swap(uLeft, uRight)
    End Select

    Dim uResult As BigUInteger
    If (uLeft.ValueLength = 1) AndAlso (uRight.ValueLength = 1) Then
        Dim wLeftInverse As UInt32, wRightInverse As UInt32
        uResult = ComputeGCD_Euclid(uLeft.DigitLowest, uRight.DigitLowest, wLeftInverse, wRightInverse)
        uLeftInverse = wLeftInverse : uRightInverse = wRightInverse
    ElseIf uLeft.ValueLength <= 2 Then
        uResult = ComputeGCD_Euclid(uLeft, uRight, uLeftInverse, uRightInverse)
    Else
        uResult = ComputeGCD_Lehmer(uLeft, uRight, uLeftInverse, uRightInverse, fComputeLeftInverse, fComputeRightInverse)
    End If

    If fSwap Then Swap(uLeftInverse, uRightInverse)

    Return uResult
End Function

Private Shared Function ComputeGCD_Lehmer(ByVal uLeft As BigUInteger, ByVal uRight As BigUInteger,
        ByRef uLeftInverse As BigUInteger, ByRef uRightInverse As BigUInteger, ByVal fComputeLeftInverse As Boolean, ByVal fComputeRightInverse As Boolean) As BigUInteger


    Dim uLeftCur As BigUInteger = uLeft, uRightCur As BigUInteger = uRight
    Dim uLeftInvPrev As BigUInteger = Instance.One, uRightInvPrev As BigUInteger = Instance.Zero,
        uLeftInvCur As BigUInteger = uRightInvPrev, uRightInvCur As BigUInteger = uLeftInvPrev,
        fInvInit As Boolean = False, fIterationIsEven As Boolean = True

    Dim dwLeftCur, dwRightCur As UInt64
    Dim wLeftInvPrev, wRightInvPrev, wLeftInvCur, wRightInvCur As UInt32
    Dim dwNumeratorMore, dwNumeratorLess, dwDenominatorMore, dwDenominatorLess, dwQuotientMore, dwQuotientLess As UInt64,
        wQuotient As UInt32
    Const nSubtractionThresholdBits As Byte = (5 - 1)

    Dim ndxDigitMax As Integer, fRightIsShorter As Boolean

    Dim fResultFound As Boolean = False
    Dim uRemainder As BigUInteger = uRightCur, uQuotient As BigUInteger
    Dim uTemp As BigUInteger = Nothing, dwTemp, dwTemp2 As UInt64

    Do While uLeftCur.ValueLength > 2

        ndxDigitMax = uLeftCur.ValueLength - 1 : fRightIsShorter = (uRightCur.ValueLength < uLeftCur.ValueLength)

        Dim fShorthandStep As Boolean = True, fShorthandIterationIsEven As Boolean
        If fRightIsShorter AndAlso (uLeftCur.ValueLength - uRightCur.ValueLength > 1) Then fShorthandStep = False

        If fShorthandStep Then

            dwLeftCur = uLeftCur.ValueDigits(ndxDigitMax - 1) Or (CULng(uLeftCur.ValueDigits(ndxDigitMax)) << DigitSize.Bits)
            dwRightCur = uRightCur.ValueDigits(ndxDigitMax - 1) Or If(fRightIsShorter, DigitValue.Zero, CULng(uRightCur.ValueDigits(ndxDigitMax)) << DigitSize.Bits)
            If ndxDigitMax >= 2 Then
                Dim nNormHead As Byte = GetNormalizationHead(uLeftCur.ValueDigits(ndxDigitMax))
                If nNormHead <> ByteValue.Zero Then
                    dwLeftCur = (dwLeftCur << nNormHead) Or (uLeftCur.ValueDigits(ndxDigitMax - 2) >> (DigitSize.Bits - nNormHead))
                    dwRightCur = (dwRightCur << nNormHead) Or (uRightCur.ValueDigits(ndxDigitMax - 2) >> (DigitSize.Bits - nNormHead))
                End If
            End If

            If CUInt(dwRightCur >> DigitSize.Bits) = DigitValue.Zero Then fShorthandStep = False

        End If

        If fShorthandStep Then

            ' First iteration, where overflow may occur in general formulae.

            If dwLeftCur = dwRightCur Then
                fShorthandStep = False
            Else
                If dwLeftCur = DoubleValue.Full Then dwLeftCur >>= 1 : dwRightCur >>= 1
                dwDenominatorMore = dwRightCur : dwDenominatorLess = dwRightCur + DigitValue.One
                dwNumeratorMore = dwLeftCur + DigitValue.One : dwNumeratorLess = dwLeftCur

                If (dwNumeratorMore >> nSubtractionThresholdBits) <= dwDenominatorMore Then
                    wQuotient = DigitValue.Zero
                    Do
                        wQuotient += DigitValue.One : dwNumeratorMore -= dwDenominatorMore
                    Loop While dwNumeratorMore >= dwDenominatorMore
                    dwQuotientMore = wQuotient
                Else
                    dwQuotientMore = dwNumeratorMore \ dwDenominatorMore
                    If dwQuotientMore >= DigitValue.BitHi Then fShorthandStep = False
                    wQuotient = CUInt(dwQuotientMore)
                End If

                If fShorthandStep Then
                    If (dwNumeratorLess >> nSubtractionThresholdBits) <= dwDenominatorLess Then
                        wQuotient = DigitValue.Zero
                        Do
                            wQuotient += DigitValue.One : dwNumeratorLess -= dwDenominatorLess
                        Loop While dwNumeratorLess >= dwDenominatorLess
                        dwQuotientLess = wQuotient
                    Else
                        dwQuotientLess = dwNumeratorLess \ dwDenominatorLess
                    End If
                    If dwQuotientMore <> dwQuotientLess Then fShorthandStep = False
                End If

            End If

        End If

        If fShorthandStep Then

            ' Prepare for the second iteration.
            wLeftInvPrev = DigitValue.Zero : wLeftInvCur = DigitValue.One
            wRightInvPrev = DigitValue.One : wRightInvCur = wQuotient
            dwTemp = dwLeftCur - wQuotient * dwRightCur : dwLeftCur = dwRightCur : dwRightCur = dwTemp
            fShorthandIterationIsEven = True

            fIterationIsEven = Not fIterationIsEven

            ' Other iterations, no overflow possible(?).
            Do

                If fShorthandIterationIsEven Then
                    If dwRightCur = wRightInvCur Then Exit Do
                    dwDenominatorMore = dwRightCur - wRightInvCur : dwDenominatorLess = dwRightCur + wLeftInvCur
                    dwNumeratorMore = dwLeftCur + wRightInvPrev : dwNumeratorLess = dwLeftCur - wLeftInvPrev
                Else
                    If dwRightCur = wLeftInvCur Then Exit Do
                    dwDenominatorMore = dwRightCur - wLeftInvCur : dwDenominatorLess = dwRightCur + wRightInvCur
                    dwNumeratorMore = dwLeftCur + wLeftInvPrev : dwNumeratorLess = dwLeftCur - wRightInvPrev
                End If

                If (dwNumeratorMore >> nSubtractionThresholdBits) <= dwDenominatorMore Then
                    wQuotient = DigitValue.Zero
                    Do
                        wQuotient += DigitValue.One : dwNumeratorMore -= dwDenominatorMore
                    Loop While dwNumeratorMore >= dwDenominatorMore
                    dwQuotientMore = wQuotient
                Else
                    dwQuotientMore = dwNumeratorMore \ dwDenominatorMore
                    If dwQuotientMore >= DigitValue.BitHi Then Exit Do
                    wQuotient = CUInt(dwQuotientMore)
                End If

                If (dwNumeratorLess >> nSubtractionThresholdBits) <= dwDenominatorLess Then
                    wQuotient = DigitValue.Zero
                    Do
                        wQuotient += DigitValue.One : dwNumeratorLess -= dwDenominatorLess
                    Loop While dwNumeratorLess >= dwDenominatorLess
                    dwQuotientLess = wQuotient
                Else
                    dwQuotientLess = dwNumeratorLess \ dwDenominatorLess
                End If
                If dwQuotientMore <> dwQuotientLess Then Exit Do

                dwTemp = wLeftInvPrev + wQuotient * wLeftInvCur : dwTemp2 = wRightInvPrev + wQuotient * wRightInvCur
                If (dwTemp >= DigitValue.BitHi) OrElse (dwTemp2 >= DigitValue.BitHi) Then Exit Do
                wLeftInvPrev = wLeftInvCur : wLeftInvCur = CUInt(dwTemp)
                wRightInvPrev = wRightInvCur : wRightInvCur = CUInt(dwTemp2)
                dwTemp = dwLeftCur - wQuotient * dwRightCur : dwLeftCur = dwRightCur : dwRightCur = dwTemp
                fShorthandIterationIsEven = Not fShorthandIterationIsEven

                fIterationIsEven = Not fIterationIsEven

            Loop

        End If

        If (Not fShorthandStep) OrElse (wRightInvPrev = DigitValue.Zero) Then
            ' Longhand step.

            uQuotient = ComputeDivisionLonghand(uLeftCur, uRightCur, uTemp) : If uTemp.IsZero Then fResultFound = True : Exit Do
            uRemainder = uTemp

            fIterationIsEven = Not fIterationIsEven
            If fComputeLeftInverse Then
                uTemp = uLeftInvPrev + uQuotient * uLeftInvCur : uLeftInvPrev = uLeftInvCur : uLeftInvCur = uTemp
            End If
            If fComputeRightInverse Then
                uTemp = uRightInvPrev + uQuotient * uRightInvCur : uRightInvPrev = uRightInvCur : uRightInvCur = uTemp
            End If
            fInvInit = True

            uLeftCur = uRightCur : uRightCur = uRemainder

        Else
            ' Shorthand step finalization.

            If Not fInvInit Then
                If fComputeLeftInverse Then uLeftInvPrev = wLeftInvPrev : uLeftInvCur = wLeftInvCur
                If fComputeRightInverse Then uRightInvPrev = wRightInvPrev : uRightInvCur = wRightInvCur
                fInvInit = True
            Else
                If fComputeLeftInverse Then ComputeFusedMulMulAdd(uLeftInvPrev, uLeftInvCur, wLeftInvPrev, wLeftInvCur, wRightInvPrev, wRightInvCur)
                If fComputeRightInverse Then ComputeFusedMulMulAdd(uRightInvPrev, uRightInvCur, wLeftInvPrev, wLeftInvCur, wRightInvPrev, wRightInvCur)
            End If

            ComputeFusedMulMulSub(uLeftCur, uRightCur, wLeftInvPrev, wLeftInvCur, wRightInvPrev, wRightInvCur, fShorthandIterationIsEven)

        End If

    Loop

    ' Final rounds: numbers are quite short now.
    If Not fResultFound Then

        ndxDigitMax = uLeftCur.ValueLength - 1 : fRightIsShorter = (uRightCur.ValueLength < uLeftCur.ValueLength)
        If ndxDigitMax = 0 Then
            dwLeftCur = uLeftCur.ValueDigits(0)
            dwRightCur = uRightCur.ValueDigits(0)
        Else
            dwLeftCur = uLeftCur.ValueDigits(0) Or (CULng(uLeftCur.ValueDigits(1)) << DigitSize.Bits)
            dwRightCur = uRightCur.ValueDigits(0) Or If(fRightIsShorter, DigitValue.Zero, CULng(uRightCur.ValueDigits(1)) << DigitSize.Bits)
        End If

        Do While dwLeftCur >= DigitValue.BitHi

            Dim dwRemainder As UInt64 = dwLeftCur

            If (dwRemainder >> nSubtractionThresholdBits) <= dwRightCur Then
                wQuotient = DigitValue.Zero
                Do
                    wQuotient += DigitValue.One : dwRemainder -= dwRightCur
                Loop While dwRemainder >= dwRightCur
                dwQuotientMore = wQuotient
            Else
                dwQuotientMore = dwLeftCur \ dwRightCur
                dwRemainder = dwLeftCur - dwQuotientMore * dwRightCur
            End If

            If dwRemainder = DigitValue.Zero Then fResultFound = True : Exit Do


            fIterationIsEven = Not fIterationIsEven
            If dwQuotientMore < DigitValue.BitHi Then
                wQuotient = CUInt(dwQuotientMore)
                If fComputeLeftInverse Then ComputeFusedMulAdd(uLeftInvPrev, uLeftInvCur, wQuotient)
                If fComputeRightInverse Then ComputeFusedMulAdd(uRightInvPrev, uRightInvCur, wQuotient)
            Else
                If fComputeLeftInverse Then
                    uTemp = uLeftInvPrev + dwQuotientMore * uLeftInvCur : uLeftInvPrev = uLeftInvCur : uLeftInvCur = uTemp
                End If
                If fComputeRightInverse Then
                    uTemp = uRightInvPrev + dwQuotientMore * uRightInvCur : uRightInvPrev = uRightInvCur : uRightInvCur = uTemp
                End If
            End If

            dwLeftCur = dwRightCur : dwRightCur = dwRemainder

        Loop

        If fResultFound Then

            uRightCur = dwRightCur

        Else

            ' Final rounds: both numbers have only one digit now, and this digit has MS-bit unset.
            Dim wLeftCur As UInt32 = CUInt(dwLeftCur), wRightCur As UInt32 = CUInt(dwRightCur)

            Do

                Dim wRemainder As UInt32 = wLeftCur

                If (wRemainder >> nSubtractionThresholdBits) <= wRightCur Then
                    wQuotient = DigitValue.Zero
                    Do
                        wQuotient += DigitValue.One : wRemainder -= wRightCur
                    Loop While wRemainder >= wRightCur
                Else
                    wQuotient = wLeftCur \ wRightCur
                    wRemainder = wLeftCur - wQuotient * wRightCur
                End If

                If wRemainder = DigitValue.Zero Then fResultFound = True : Exit Do

                fIterationIsEven = Not fIterationIsEven
                If fComputeLeftInverse Then ComputeFusedMulAdd(uLeftInvPrev, uLeftInvCur, wQuotient)
                If fComputeRightInverse Then ComputeFusedMulAdd(uRightInvPrev, uRightInvCur, wQuotient)

                wLeftCur = wRightCur : wRightCur = wRemainder

            Loop

            uRightCur = wRightCur

        End If


    End If

    If fComputeLeftInverse Then
        uLeftInverse = If(fIterationIsEven, uRight - uLeftInvCur, uLeftInvCur)
    End If
    If fComputeRightInverse Then
        uRightInverse = If(fIterationIsEven, uRightInvCur, uLeft - uRightInvCur)
    End If

    Return uRightCur
End Function

''' <remarks>All word-sized parameters must have their most-significant bit unset.</remarks>
Private Shared Sub ComputeFusedMulMulAdd(
        ByRef uLeftInvPrev As BigUInteger, ByRef uLeftInvCur As BigUInteger,
        ByVal wLeftInvPrev As UInt32, ByVal wLeftInvCur As UInt32, ByVal wRightInvPrev As UInt32, ByVal wRightInvCur As UInt32)

    Dim ndxDigitMaxPrev As Integer = uLeftInvPrev.ValueLength - 1, ndxDigitMaxCur As Integer = uLeftInvCur.ValueLength - 1,
        ndxDigitMaxNew As Integer = ndxDigitMaxCur + 1

    Dim awLeftInvPrev() As UInt32 = uLeftInvPrev.ValueDigits, awLeftInvCur() As UInt32 = uLeftInvCur.ValueDigits
    Dim awLeftInvPrevNew(0 To ndxDigitMaxNew) As UInt32, awLeftInvCurNew(0 To ndxDigitMaxNew) As UInt32
    Dim dwResult As UInt64, wCarryLeftPrev As UInt32 = DigitValue.Zero, wCarryLeftCur As UInt32 = DigitValue.Zero
    Dim wDigitLeftInvPrev, wDigitLeftInvCur As UInt32

    For ndxDigit As Integer = 0 To ndxDigitMaxPrev
        wDigitLeftInvPrev = awLeftInvPrev(ndxDigit) : wDigitLeftInvCur = awLeftInvCur(ndxDigit)

        dwResult = wCarryLeftPrev + wLeftInvPrev * CULng(wDigitLeftInvPrev) + wRightInvPrev * CULng(wDigitLeftInvCur)
        awLeftInvPrevNew(ndxDigit) = CUInt(dwResult And DigitValue.Full) : wCarryLeftPrev = CUInt(dwResult >> DigitSize.Bits)

        dwResult = wCarryLeftCur + wLeftInvCur * CULng(wDigitLeftInvPrev) + wRightInvCur * CULng(wDigitLeftInvCur)
        awLeftInvCurNew(ndxDigit) = CUInt(dwResult And DigitValue.Full) : wCarryLeftCur = CUInt(dwResult >> DigitSize.Bits)

    Next

    If ndxDigitMaxCur > ndxDigitMaxPrev Then

        For ndxDigit As Integer = ndxDigitMaxPrev + 1 To ndxDigitMaxCur
            wDigitLeftInvCur = awLeftInvCur(ndxDigit)

            dwResult = wCarryLeftPrev + wRightInvPrev * CULng(wDigitLeftInvCur)
            awLeftInvPrevNew(ndxDigit) = CUInt(dwResult And DigitValue.Full) : wCarryLeftPrev = CUInt(dwResult >> DigitSize.Bits)

            dwResult = wCarryLeftCur + wRightInvCur * CULng(wDigitLeftInvCur)
            awLeftInvCurNew(ndxDigit) = CUInt(dwResult And DigitValue.Full) : wCarryLeftCur = CUInt(dwResult >> DigitSize.Bits)

        Next

    End If

    If wCarryLeftPrev <> DigitValue.Zero Then awLeftInvPrevNew(ndxDigitMaxNew) = wCarryLeftPrev
    If wCarryLeftCur <> DigitValue.Zero Then awLeftInvCurNew(ndxDigitMaxNew) = wCarryLeftCur

    uLeftInvPrev = New BigUInteger(awLeftInvPrevNew) : uLeftInvCur = New BigUInteger(awLeftInvCurNew)

End Sub

''' <remarks>All word-sized parameters must have their most-significant bit unset.</remarks>
Private Shared Sub ComputeFusedMulMulSub(
        ByRef uLeftCur As BigUInteger, ByRef uRightCur As BigUInteger,
        ByVal wLeftInvPrev As UInt32, ByVal wLeftInvCur As UInt32, ByVal wRightInvPrev As UInt32, ByVal wRightInvCur As UInt32,
        ByVal fShorthandIterationIsEven As Boolean)

    Dim ndxDigitMax As Integer = uLeftCur.ValueLength - 1,
        fRightIsShorter As Boolean = (uRightCur.ValueLength < uLeftCur.ValueLength),
        ndxDigitStop As Integer = If(fRightIsShorter, ndxDigitMax - 1, ndxDigitMax)

    Dim awLeftCur() As UInt32 = uLeftCur.ValueDigits, awRightCur() As UInt32 = uRightCur.ValueDigits
    Dim awLeftNew(0 To ndxDigitMax) As UInt32, awRightNew(0 To ndxDigitStop) As UInt32
    Dim iTemp As Int64, wCarryLeft As Int32 = 0I, wCarryRight As Int32 = 0I
    Dim wDigitLeftCur, wDigitRightCur As UInt32

    If fShorthandIterationIsEven Then

        For ndxDigit As Integer = 0 To ndxDigitStop
            wDigitLeftCur = awLeftCur(ndxDigit) : wDigitRightCur = awRightCur(ndxDigit)
            iTemp = wCarryLeft + CLng(wDigitRightCur) * wRightInvPrev - CLng(wDigitLeftCur) * wLeftInvPrev
            awLeftNew(ndxDigit) = CUInt(iTemp And DigitValue.Full) : wCarryLeft = CInt(iTemp >> DigitSize.Bits)
            iTemp = wCarryRight + CLng(wDigitLeftCur) * wLeftInvCur - CLng(wDigitRightCur) * wRightInvCur
            awRightNew(ndxDigit) = CUInt(iTemp And DigitValue.Full) : wCarryRight = CInt(iTemp >> DigitSize.Bits)
        Next
        If fRightIsShorter Then
            wDigitLeftCur = awLeftCur(ndxDigitMax)
            iTemp = wCarryLeft - CLng(wDigitLeftCur) * wLeftInvPrev
            awLeftNew(ndxDigitMax) = CUInt(iTemp And DigitValue.Full)
        End If

    Else

        For ndxDigit As Integer = 0 To ndxDigitStop
            wDigitLeftCur = awLeftCur(ndxDigit) : wDigitRightCur = awRightCur(ndxDigit)
            iTemp = wCarryLeft + CLng(wDigitLeftCur) * wLeftInvPrev - CLng(wDigitRightCur) * wRightInvPrev
            awLeftNew(ndxDigit) = CUInt(iTemp And DigitValue.Full) : wCarryLeft = CInt(iTemp >> DigitSize.Bits)
            iTemp = wCarryRight + CLng(wDigitRightCur) * wRightInvCur - CLng(wDigitLeftCur) * wLeftInvCur
            awRightNew(ndxDigit) = CUInt(iTemp And DigitValue.Full) : wCarryRight = CInt(iTemp >> DigitSize.Bits)
        Next
        If fRightIsShorter Then
            wDigitLeftCur = awLeftCur(ndxDigitMax)
            iTemp = wCarryLeft + CLng(wDigitLeftCur) * wLeftInvPrev
            awLeftNew(ndxDigitMax) = CUInt(iTemp And DigitValue.Full)
        End If

    End If

    uLeftCur = New BigUInteger(awLeftNew) : uRightCur = New BigUInteger(awRightNew)

End Sub

''' <remarks>All word-sized parameters must have their most-significant bit unset.</remarks>
Private Shared Sub ComputeFusedMulAdd(ByRef uLeftInvPrev As BigUInteger, ByRef uLeftInvCur As BigUInteger, ByVal wQuotient As UInt32)

    Dim ndxDigitPrevMax As Integer = uLeftInvPrev.ValueLength - 1, ndxDigitCurMax As Integer = uLeftInvCur.ValueLength - 1,
        ndxDigitNewMax As Integer = ndxDigitCurMax + 1
    Dim awLeftInvPrev() As UInt32 = uLeftInvPrev.ValueDigits, awLeftInvCur() As UInt32 = uLeftInvCur.ValueDigits,
        awLeftInvNew(0 To ndxDigitNewMax) As UInt32
    Dim dwResult As UInt64 = DigitValue.Zero, wCarry As UInt32 = DigitValue.Zero

    For ndxDigit As Integer = 0 To ndxDigitPrevMax
        dwResult = CULng(wCarry) + awLeftInvPrev(ndxDigit) + CULng(wQuotient) * awLeftInvCur(ndxDigit)
        awLeftInvNew(ndxDigit) = CUInt(dwResult And DigitValue.Full) : wCarry = CUInt(dwResult >> DigitSize.Bits)
    Next

    For ndxDigit As Integer = ndxDigitPrevMax + 1 To ndxDigitCurMax
        dwResult = CULng(wCarry) + CULng(wQuotient) * awLeftInvCur(ndxDigit)
        awLeftInvNew(ndxDigit) = CUInt(dwResult And DigitValue.Full) : wCarry = CUInt(dwResult >> DigitSize.Bits)
    Next

    If wCarry <> DigitValue.Zero Then awLeftInvNew(ndxDigitNewMax) = wCarry

    uLeftInvPrev = uLeftInvCur : uLeftInvCur = New BigUInteger(awLeftInvNew)

End Sub

如果您想直接使用此代码，您可能需要 Visual Basic 2012 编译器来构建某些结构——我没有检查以前的版本；我也不知道最低 .Net 版本(至少 3.5 就足够了)；众所周知，已编译的应用程序可以在 Mono 上运行，尽管性能较差。我唯一可以肯定的是，人们不应该尝试使用自动 VB-to-C# 转换器，因为它们在此类主题中非常糟糕；只能靠自己的脑袋。