(16C) Modular Exponentiation

[wynen]

Computing \(base^{exponent} \bmod modulus\) for positive integer base, exponent and modulus is a common task in cryptography. The HP-16C is ideally suited to implement a binary exponentiation algorithm.

Code:


function modular_pow(base, exponent, modulus)
    if modulus = 1 then return 0
    Assert :: (modulus - 1) * (modulus - 1) does not overflow word size
    result := 1
    base := base mod modulus
    while exponent > 0
        if (exponent mod 2 = 1):
           result := (result * base) mod modulus
        exponent := exponent >> 1
        base := (base * base) mod modulus
    return result

from wikipedia: https://en.wikipedia.org/wiki/Modular_exponentiation

The HP-16C has a maximum word size of 64bit, so the maximum possible modulus should be \(2^{32}-1\) (about \(10^9\)). For larger values, the two (modular) multiplications could be implemented in a similar binary multiplication algorithm, but fortunately the HP-16C allows double precision multiplication and division. Using this feature, the maximum values for the modulus are \(2^{63}-1\).

Code:


001 43,22,A   g LBL A       x: modulus, y: exponent, z: base
002 44,0        STO 0       R0 = modulus
003 33          Rdown
004 34          x<->y       x: base, y: exponent, z: ?, t: modulus
005 43,33     g Rup
006 42,9      f RMD
007 44,1        STO 1       R1 = base mod modulus
008 43,35     g CLx
009 1           [1]         x: result = 1, y: exponent
010 44,2        STO 2       R2 = result
011 33          Rdown
012 43,22,9   g LBL 9       x: exponent
013 43,40     g x=0 
014 22,1        GTO 1       if exponent == 0 return result
015 42,b      f SR          exponent = exponent >> 1
016 43,6,4    g F? 4        Carry set?
017 43,2      g x<0         yes, but x<0 is always false (because of shift right)
018 22,2        GTO 2       Carry not set
019 45,2        RCL 2       Carry set
020 45,1        RCL 1       x: modulus, y: result, z: exponent
021 43,20     g DBLx        x: xy..., y: ...xy, z: exponent
022 45,0        RCL 0
023 43,9      g DBLR        x: xy mod modulus, z: exponent
024 44,2        STO 2       result = (result * base) mod modulus
025 33          Rdown       x: exponent
026 43,22,2   g LBL 2
027 45,1        RCL 1
028 36          ENTER       x: base, y: base, z: exponent
029 43,20     g DBLx        x: xy..., y: ...xy, z: exponent
030 45,0        RCL 0
031 43,9      g DBLR        x: xy mod modulus, z: exponent
032 44,1        STO 1       base = (base * base) mod modulus
033 33          Rdown
034 22,9        GTO 9       x: exponent
035 43,22,1   g LBL 1
036 45,2        RCL 2
037 43,21       RTN

Quote:usage:
[f][UNSGN] 64 [f][WSIZE] base [ENTER] exponent [ENTER] modulus [GSB][A]

Register:
R0: modulus
R1: base
R2: result

Best regards
Hartmut