Lehmer的扩展欧几里得算法实现

6
在实现自己的BigInteger时,我被扩展的GCD算法卡住了,这对于找到模乘逆元是基本的。由于众所周知的欧几里得方法执行速度太慢,而混合和二进制算法只快了5-10倍,因此选择了Lehmer对经典算法的修改。但问题在于,当涉及到描述Lehmer时,我发现所有的书籍(Knuth、应用密码学手册、互联网等)都有同样的缺点:
1.解释基于几个技巧: -输入数字总是相同长度的; -抽象CPU具有带符号寄存器,可以容纳数字和符号; -抽象CPU具有半无限的寄存器,即它永远不会溢出。
2.只提供基本的GCD算法,没有关注逆因子。
关于第一个问题,我最初很惊讶地发现找不到任何实际的实现(不要指向GNU MP库-它不是学习的来源),但最终通过反编译Microsoft的.Net 4.0实现得到了灵感,这显然是基于Jebelean的论文“一种用于查找长整数GCD的双位Lehmer-Euclid算法”的思想。结果函数很大,看起来吓人,但效果非常好。
但是微软的库只提供基本功能,不计算任何余因子。准确地说,在简写步骤中会计算一些余因子,并且在第一步时这些余因子仅是初始余因子,但在执行长手步骤之后它们就不再匹配了。我目前的解决方案是在“替代”余因子与“真实”余因子(除了第一步)并行更新,但这使得性能下降到了零以下:该函数现在仅比基本模式的二进制方法快25-50%。因此,问题在于,尽管输入数字仅在长手步骤中完全更新,但余因子也在每个简写步骤的迭代中更新,从而破坏了Lehmer方法的任何好处。
为了加快速度,我实现了一个“融合乘加”函数,因为“融合乘-乘-减”确实有助于更新输入数字,但这次影响微不足道。另一个改进基于通常只需要一个余子式的事实,所以另一个余子式可以完全不计算。这应该将开销减少一半(甚至更多,因为第二个数字通常比第一个数字小得多),但在实践中,开销仅减少了预期的25%到50%。
因此,我的问题归结为以下内容:
1. 是否有关于Lehmer算法的全面解释,与在现实世界硬件上的实际实现相关(使用有限大小的无符号字)? 2. 同上,但涉及扩展GCD计算。
尽管我对基本算法的性能感到满意,但相反的情况适用于我主要使用的扩展操作模式。

只是问一下:这是你正在写的论文,你想让我们提供文献搜索吗? - andy256
1
不,这是我正在编写的一个辅助库(BigUInteger / BigInteger 实现),作为需要验证输入文档数字签名的应用程序的一部分。它已经可以工作并具有可接受的性能,但当然速度越快越好。我甚至编写了一个彻底的基准测试套件,以比较每个运算符和函数的性能与参考 Microsoft 的 .Net 实现。 - Anton Samsonov
好的。完成后,它应该是一篇论文 :-) 我对这种算法研究的担忧在于测试和认证。 - andy256
我并不聪明,不能发明新算法——这就是为什么我要求帮助。至于实现已知方法,当然会出错,但单元测试可以帮助解决这个问题,混合使用预计算的情况和纯随机数据(如果可能的话,在运行时与另一个实现进行验证)。在我的情况下不需要认证,所以这对我来说不是问题:如果应用程序按照预期工作,则假定它是无错误的,直到有证据表明存在问题。 - Anton Samsonov
你说你“通过反编译微软的 .Net 4.0 实现获得了灵感,...结果函数很大,...但非常有效”。请记住,微软的代码既不是公共领域,也不是开源的,它是他们的知识产权。查看它的反编译以了解其功能可能是可以的(我不确定),但未经他们许可提取和使用它绝对是不可以的。 - RBarryYoung
“取灵感”并不是以任何方式“复制粘贴”。在最初的阶段,我研究了Mono和.Net,但因为它们都没有实现高级技术,如速记和预计算,所以我选择了自己的方式。Lehmer's是唯一的例外(仅限于Microsoft),但它缺乏余子计算。此外,我没有看到过确切的Microsoft代码,而是一个糟糕的反编译源代码(请自行尝试Mono.Cecil)。更不用说.Net中的内部数字表示与我的迥然不同。简而言之,实现完全不同,没有什么可担心的。 - Anton Samsonov
1个回答

5
最后,我咨询了一位数学家,他很快就找到了正确的公式——非常类似于我自己尝试的那些,但仍然略有不同。这使得只需在手写步骤上更新余因子,同时输入数字完全更新。
然而,令我大为惊讶的是,这个措施单独对性能影响很小。只有当我将其重新实现为“融合(A×X + B×Y)”时,速度提升才变得明显。当计算两个余因子时,它现在以纯Lehmer GCD算法为基础,对于5位数字和32K位数字,它运行速度分别为56%和34%; 对于单个余因子,速率分别为70%和52%。使用之前的实现,对于两个余因子,它仅为33%至7%,对于单个余因子,它为47%至14%,所以我的不满显而易见。
至于像andy256建议的撰写论文,以便其他实现者不会遇到同样的问题,这并不容易。当我向数学家解释我的当前实现时,我已经写了一篇“小”论文,它相当快地超过了三张A4大小的纸——仅包含基本思想,没有详细的解释、溢出检查、分支、展开等。现在我部分理解了为什么Knuth和其他人使用肮脏的技巧来保持故事简短。目前,我不知道如何在不失彻底性的情况下实现相同水平的简洁;也许需要几个大流程图结合评论。
更新。看起来我无法在近期完成并发布库(仍然没有运气理解Newton-Raphson除法和Montgomery约简),所以我将简单地发布结果函数供感兴趣的人使用。
代码不包括明显的函数,如ComputeGCD_Euclid和通用例程,如ComputeDivisionLonghand。代码也缺少任何注释(有少数例外)——如果您想要理解它并将其集成到自己的库中,您应该已经熟悉Lehmer的思想和上面提到的双位数手写技术。
我的库中数字表示的概述。数字是32位无符号整数,因此可以在需要时使用64位无符号和有符号算术。数字按从最低有效位(LSB)到最高有效位的顺序存储在一个普通数组(ValueDigits)中,实际大小明确存储(ValueLength),即函数尝试预测结果大小,但不会在计算后优化内存消耗。对象是类型(.Net中的struct),但它们引用数字数组;因此,对象是不变的,即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到C#翻译器,因为它们在这方面非常糟糕; 只能依靠您自己的头脑。

非常感谢您发布这篇文章。尽管它很复杂且难以理解,但仍然是互联网上可以找到的极少数关于如何实现这种东西的工作示例之一。 - RBarryYoung
1
干得好!我也需要一个更快的反转算法,但是我将经典的Lehmer算法(对我来说更直观)与建议的两位数调整相结合(这是一个很好的提升)。如果您感兴趣,可以在这里找到它。 - Axel Heer
有人将这个示例转换为PowerShell了吗?还是我应该另外提一个问题? - Carsten

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接