Modular (Binary) Exponentiation example

[cbruce]

.
Modular Exponentiation is an operation that is used quite frequently
in areas such as cryptography. The Diffie-Hellman method is a prime
(pun intended) example.

I actually use the below function with the parameters shown for some
unique purposes... but it has limited use in most real world
programming; due to the fact that it can only handle the relatively
small values shown. Remember, something like the Diffie-Hellman
method is normally using 300 digit or larger numbers... that means
that DH calls for a bignum library.

So just consider this an educational example of how binary
exponentiation works.

Note:
If anybody wants to help me out by explaining why the intermediate
result inside the function is only updated when working with the ones
bits of the exponent - feel free! I wrote this a few years ago, and
every time I use it I wonder about that.

I've seen lots of examples on the web, but I've never seen code that
even comes close to commenting correctly about why the individual
steps are doing what they do... which is ridiculous considering that
there are so few steps involved and that it is a function that is in
great use. Hopefully, this will help someone else.

Code: Select all

CONST DBGPRINT  = false

' -------------------------------
'    Modular Exponentiation
' -------------------------------
Function Modular_Exponentiation(base_num as ulong, exponent as ulongint, modulus as ulong) as ulong
    ' Modular Exponentiation by Repeated Squaring - (Binary Exponentiation)
    '
    ' Modular Exponentiation - Wikipedia
    ' https://en.wikipedia.org/wiki/Modular_exponentiation
    '
    ' We are able to directly calculate things, like a 32-bit number 
    ' raised to the power of a 64-bit number and then take the 
    ' modulo of that result by a 32-bit number, by utilizing a 
    ' numerical method called "repeated squaring" to break the 
    ' exponentiation down into a series of smaller exponentiations 
    ' that are manageable without the need for a BIGNUM library.
    '    
    dim b1 as ulongint = base_num
    dim e1 as ulongint = exponent
    dim result as ulongint = 1
'                                                ' This temp var is just so I can print which 
'                                                ' bit of the exponent I'm working on below.
'                                                dim i_bit as Integer     = 0
'                                                IF DBGPRINT THEN
'                                                    print "INIT base             = "; b1
'                                                    print "INIT exponent (priv)  = "; e1
'                                                    print "INIT modulus          = "; modulus
'                                                    print "INIT result           = "; result
'                                                END IF
    '
    b1 = b1 MOD modulus
'                                                IF DBGPRINT THEN
'                                                    print "b1 MOD modulus        = "; b1
'                                                    print
'                                                END IF
    '
    ' Binary Exponentiation... 
    ' ... We're working with the exponent as a series of bits. 
    ' We will be doing a SHR on the exponent in every pass of the 
    ' loop in order to get at the exponent's next bit.
    ' When all of the ones bits are gone (leaving exponent = 0), 
    ' then we are done with the exponentiation.
    do while exponent > 0
'                                                IF DBGPRINT THEN
'                                                    print "DO exponent           = "; bin(exponent,32)
'                                                END IF
        ' Binary Exponentiation... 
        ' ... We're working with the exponent as a series of bits. 
        ' If bit zero of the remaining exponent bits is a one bit...
        if (exponent and 1) then
            ' I'm not a maths guy... I don't know why we are only 
            ' updating the result when working with the ones bits
            ' of the exponent.
            ' The equation just says to do this now.
            result = (result * b1) MOD modulus
'                                                IF DBGPRINT THEN
'                                                    print "result*b1 MOD modulus = "; result
'                                                END IF
        end if
'                                                IF DBGPRINT THEN
'                                                    print "b1                    = "; b1
'                                                END IF
        '
'                                                IF DBGPRINT THEN
'                                                    ' We want to see both the square and the MOD results.
'                                                    print "(b1 ^ 2)              = "; (b1 * b1)
'                                                END IF
        b1 = (b1 * b1) MOD modulus
'                                                IF DBGPRINT THEN
'                                                    ' We want to see both the square and the MOD results.
'                                                    print "b1 MOD modulus        = "; b1
'                                                END IF
        '
        ' Binary Exponentiation... 
        ' ... We're working with the exponent as a series of bits. 
        ' Shift the exponent right one bit to get at the next bit.
        exponent = exponent shr 1
'                                                IF DBGPRINT THEN
'                                                    ' Track which bit of the exponent is next...
'                                                    i_bit += 1
'                                                    print i_bit; tab(3); " - exponent SHR     = "; bin(exponent,32)
'                                                    print
'                                                END IF
    loop
'                                                IF DBGPRINT THEN
'                                                    print "RESULT ................ "; result; "  "; ucase(hex(result))
'                                                    print
'                                                    print
'                                                END IF
    '
    return result
end function

' ------------------------------------------------------------

dim base_num as ulong    = 5
dim exponent as ulongint = 2982377858398747
' This is the largest 32-bit prime. It will result in the largest 
' 32-bit result set.
dim modulus as ulong     = 4294967291
'
dim result as ulong      = 0
'
result = Modular_Exponentiation(base_num, exponent, modulus)
print "result = "; result
end

[dodicat]

Dunno cbruce, it just has to be an odd number.
I suppose the ones bits of the exponent is a fancy way of saying an odd number?
Took me a minute to get past that bit, and spent a little while with google, like yourself, couldn't find an easy answer.

Here is a string way.
Agrees with your own.
Also solves the rosetta code question https://rosettacode.org/wiki/Modular_exponentiation#C without gmp or whatever.

Code: Select all

 

Function _divide(n1 As String,n2 As String,decimal_places As Integer=10,dpflag As String="s") As String
    Dim As String number=n1,divisor=n2
    dpflag=Lcase(dpflag)
    'For MOD
    Dim As Integer modstop
    If dpflag="mod" Then 
        If Len(n1)<Len(n2) Then Return n1
        If Len(n1)=Len(n2) Then
            If n1<n2 Then Return n1
        End If
        modstop=Len(n1)-Len(n2)+1
    End If
    If dpflag<>"mod" Then
        If dpflag<>"s"  Then dpflag="raw" 
    End If
    Dim runcount As Integer
    '_______  LOOK UP TABLES ______________
    Dim Qmod(0 To 19) As Ubyte
    Dim bool(0 To 19) As Ubyte
    For z As Integer=0 To 19
        Qmod(z)=(z Mod 10+48)
        bool(z)=(-(10>z))
    Next z
    Dim answer As String   'THE ANSWER STRING  
    
    '_______ SET THE DECIMAL WHERE IT SHOULD BE AT _______
    Dim As String part1,part2
    #macro set(decimal)
    #macro insert(s,char,position)
    If position > 0 And position <=Len(s) Then
        part1=Mid(s,1,position-1)
        part2=Mid(s,position)
        s=part1+char+part2
    End If
    #endmacro
    insert(answer,".",decpos)
    answer=thepoint+zeros+answer
    If dpflag="raw" Then
        answer=Mid(answer,1,decimal_places)
    End If
    #endmacro
    '______________________________________________
    '__________ SPLIT A STRING ABOUT A CHARACTRR __________
    Dim As String var1,var2
    Dim pst As Integer
    #macro split(stri,char,var1,var2)
    pst=Instr(stri,char)
    var1="":var2=""
    If pst<>0 Then
        var1=Rtrim(Mid(stri,1,pst),".")
        var2=Ltrim(Mid(stri,pst),".")
    Else
        var1=stri
    End If
    #endmacro
    
    #macro Removepoint(s)
    split(s,".",var1,var2)
    #endmacro
    '__________ GET THE SIGN AND CLEAR THE -ve __________________
    Dim sign As String
    If Left(number,1)="-" Xor Left (divisor,1)="-" Then sign="-"
    If Left(number,1)="-" Then  number=Ltrim(number,"-")
    If Left (divisor,1)="-" Then divisor=Ltrim(divisor,"-")
    
    'DETERMINE THE DECIMAL POSITION BEFORE THE DIVISION
    Dim As Integer lennint,lenddec,lend,lenn,difflen
    split(number,".",var1,var2)
    lennint=Len(var1)
    split(divisor,".",var1,var2)
    lenddec=Len(var2)
    
    If Instr(number,".") Then 
        Removepoint(number)
        number=var1+var2
    End If
    If Instr(divisor,".") Then 
        Removepoint(divisor)
        divisor=var1+var2
    End If
    Dim As Integer numzeros
    numzeros=Len(number)
    number=Ltrim(number,"0"):divisor=Ltrim (divisor,"0")
    numzeros=numzeros-Len(number)
    lend=Len(divisor):lenn=Len(number)
    If lend>lenn Then difflen=lend-lenn
    Dim decpos As Integer=lenddec+lennint-lend+2-numzeros 'THE POSITION INDICATOR
    Dim _sgn As Byte=-Sgn(decpos)
    If _sgn=0 Then _sgn=1
    Dim As String thepoint=String(_sgn,".") 'DECIMAL AT START (IF)
    Dim As String zeros=String(-decpos+1,"0")'ZEROS AT START (IF) e.g. .0009
    If dpflag<>"mod" Then
        If Len(zeros) =0 Then dpflag="s"
    End If
    Dim As Integer runlength
    If Len(zeros) Then 
        runlength=decimal_places
        answer=String(Len(zeros)+runlength+10,"0")
        If dpflag="raw" Then 
            runlength=1
            answer=String(Len(zeros)+runlength+10,"0")
            If decimal_places>Len(zeros) Then
                runlength=runlength+(decimal_places-Len(zeros))
                answer=String(Len(zeros)+runlength+10,"0")
            End If
        End If
        
    Else
        decimal_places=decimal_places+decpos
        runlength=decimal_places
        answer=String(Len(zeros)+runlength+10,"0")
    End If
    '___________DECIMAL POSITION DETERMINED  _____________
    
    'SET UP THE VARIABLES AND START UP CONDITIONS
    number=number+String(difflen+decimal_places,"0")
    Dim count As Integer
    Dim temp As String
    Dim copytemp As String
    Dim topstring As String
    Dim copytopstring As String
    Dim As Integer lenf,lens
    Dim As Ubyte takeaway,subtractcarry
    Dim As Integer n3,diff
    If Ltrim(divisor,"0")="" Then Return "Error :division by zero"   
    lens=Len(divisor)
    topstring=Left(number,lend)
    copytopstring=topstring
    Do
        count=0
        Do
            count=count+1
            copytemp=temp
            
            Do
                '___________________ QUICK SUBTRACTION loop _________________              
                
                lenf=Len(topstring)
                If  lens<lenf=0 Then 'not
                    If Lens>lenf Then
                        temp= "done"
                        Exit Do
                    End If
                    If divisor>topstring Then 
                        temp= "done"
                        Exit Do
                    End If
                End If
                
                diff=lenf-lens
                temp=topstring
                subtractcarry=0
                
                For n3=lenf-1 To diff Step -1
                    takeaway= topstring[n3]-divisor[n3-diff]+10-subtractcarry
                    temp[n3]=Qmod(takeaway)
                    subtractcarry=bool(takeaway)
                Next n3 
                If subtractcarry=0 Then Exit Do
                If n3=-1 Then Exit Do
                For n3=n3 To 0 Step -1 
                    takeaway= topstring[n3]-38-subtractcarry
                    temp[n3]=Qmod(takeaway)
                    subtractcarry=bool(takeaway)
                    If subtractcarry=0 Then Exit Do
                Next n3
                Exit Do
                
            Loop 'single run
            temp=Ltrim(temp,"0")
            If temp="" Then temp= "0"
            topstring=temp
        Loop Until temp="done"
        ' INDIVIDUAL CHARACTERS CARVED OFF ________________       
        runcount=runcount+1
        If count=1 Then
            topstring=copytopstring+Mid(number,lend+runcount,1)
        Else
            topstring=copytemp+Mid(number,lend+runcount,1)
        End If
        copytopstring=topstring
        topstring=Ltrim(topstring,"0")
        If dpflag="mod" Then
            If runcount=modstop Then 
                If topstring="" Then Return "0"
                Return Mid(topstring,1,Len(topstring)-1)
            End If
        End If
        answer[runcount-1]=count+47
        If topstring="" And runcount>Len(n1)+1 Then
            Exit Do
        End If
    Loop Until runcount=runlength+1
    
    ' END OF RUN TO REQUIRED DECIMAL PLACES
    set(decimal) 'PUT IN THE DECIMAL POINT
    'THERE IS ALWAYS A DECIMAL POINT SOMEWHERE IN THE ANSWER
    'NOW GET RID OF IT IF IT IS REDUNDANT
    answer=Rtrim(answer,"0")
    answer=Rtrim(answer,".")
    answer=Ltrim(answer,"0")
    If answer="" Then Return "0"
    Return sign+answer
End Function

Dim Shared As Integer _Mod(0 To 99),_Div(0 To 99)
For z As Integer=0 To 99:_Mod(z)=(z Mod 10+48):_Div(z)=z\10:Next
    
    Function Qmult(a As String,b As String) As String
        Var flag=0,la = Len(a),lb = Len(b)
        If Len(b)>Len(a) Then flag=1:Swap a,b:Swap la,lb
        Dim As Ubyte n,carry,ai
        Var c =String(la+lb,"0")
        For i As Integer =la-1 To 0 Step -1
            carry=0:ai=a[i]-48
            For j As Integer =lb-1 To 0 Step -1
                n = ai * (b[j]-48) + (c[i+j+1]-48) + carry
                carry =_Div(n):c[i+j+1]=_Mod(n)
            Next j
            c[i]+=carry
        Next i
        If flag Then Swap a,b
        Return Ltrim(c,"0")
    End Function
    '=======================================================================
    #define mod_(a,b) _divide((a),(b),,"mod")
    #define div_(a,b) iif(len((a))<len((b)),"0",_divide((a),(b),-2))
    
    Function Modular_Exponentiation(base_num As String, exponent As String, modulus As String) As String
        Dim b1 As String = base_num
        Dim e1 As String = exponent
        Dim As String result="1"
        b1 = mod_(b1,modulus)
        Do While e1 <> "0"
            Var L=Len(e1)-1
            If e1[L] And 1 Then
                result=mod_(Qmult(result,b1),modulus)
            End If
            b1=mod_(qmult(b1,b1),modulus)
            e1=div_(e1,"2")
        Loop
        Return result
    End Function
    
    Dim base_num As String    = "5"
    Dim exponent As String = "2982377858398747"
    Dim modulus As String     = "4294967291"
    Var ans=Modular_Exponentiation(base_num,exponent,modulus)
    Print
    Print ans
    
    base_num="2988348162058574136915891421498819466320163312926952423791023078876139"
    exponent="2351399303373464486466122544523690094744975233415544072992656881240319"
    modulus="10000000000000000000000000000000000000000"
    ans=Modular_Exponentiation(base_num,exponent,modulus)
    Print ans
    Print "Done"
    Sleep
     
 

[cbruce]

dodicat, that string version is a really cool one! Thanks!

[integer]

This is not a full (reasonable) explanation, but it is a start.

positional notation:
decimal: 70506 --> 7*10^4 + 0*10^3 + 5*10^2 + 0*10^1 + 6*10^0
7*10000 + 0*1000 + 8*100 + 0*10 + 6*1 == 70506

binary: 10101 --> 1*2^4 + 0*2^3 + 1*2^2 + 0*2^1 + 1*2^0
16 + 0 + 4 + 0 + 1 == 21

This is the SUM of the values.


Exponentiation is the PRODUCT of the positional values

21 decimal == 10101 binary

in decimal: 3^(21)
in binary: 3^(10101)

1 * 3^16 = 43046721
0 * 3^8 = 0 not used
1 * 3^4 = 81
0 * 3^2 = 0 not used
1 * 3^1 = 3

43046721 * 81 * 3 == 3^21 = 10460353203

[cbruce]

integer, I don't "grok" it, but at least I can see it and can start to wrap my head around it. Thanks!