6

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

  1. 解释基于几个技巧:
    • 输入数字始终具有相同的长度;
    • 抽象 CPU 有带符号的寄存器,可以同时保存数字和符号;
    • 抽象 CPU 具有半无限的寄存器,即它永远不会溢出。
  2. 只提供了基本的 GCD 算法,没有关注逆辅因子。

至于第一个问题,我最初对找不到任何真实世界的实现感到惊讶(不要将我指向 GNU MP 库——它不是学习的来源),但最终通过反编译微软的实现获得了灵感来自.Net 4.0,这显然是基于论文“<a href="http://www.csie.nuk.edu.tw/~cychen/gcd/A%20double-digit%20Lehmer-Euclid% 20algorithm%20for%20finding%20the%20GCD%20of%20long%20integers.pdf" rel="nofollow">一种用于查找长整数 GCD 的两位数 Lehmer-Euclid 算法,作者 Jebelean。生成的函数很大,看起来很吓人,但效果很好。

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

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

因此,我的问题归结为:

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

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

4

1 回答 1

5

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

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

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


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

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

我的图书馆中数字表示的概述。数字是 32 位无符号整数,因此可以在需要时使用 64 位无符号和有符号算术。数字按从最低到最高有效 (LSB) 的顺序存储在普通数组 ( ValueDigits) 中,显式存储实际大小 ( ValueLength),即函数尝试预测结果大小,但不优化计算后的内存消耗。对象是类型(struct在 .Net 中),但它们引用数字数组;因此,对象是不变的,即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# 翻译器,因为它们在这样的主题中非常糟糕;只靠自己的脑袋。

于 2013-08-23T17:46:58.063 回答