[Arjen Markus] 19 november 2002. In reaction to a discussion in the Tclers' chat room:

A function that is used quite regularly in cryptographic applications is:

    P(n,r,m) = n^r mod m

(n^r denotes n to the power r, all three parameters are positive integers)

This function does not exist as such within the [expr] command, but it is
simple to create a procedure that evaluates this function efficiently.

The parameters for this function are usually large enough to make a direct
calculation impossible. But the trick is to recognise that n^r can be written as n^(r/2)*n^(r/2) when r is even. If r is odd, then write it as:

    n^r = n * (n^(r-1))

and then we have the first case back.

Furthermore, the modulo function is very well behaved:

   (h * g) mod m = ( (h mod m) * (g mod m) ) mod m

Given these observations, here is a recursive proc that evaluates the
above function in at most 2*log(r) steps. This can probably be improved,
but I am a bit lazy.

----
  proc powm {n r m} {
     #
     # Take care of the trivial cases first (they also stop the recursion)
     #
     if { $r == 0 } {
        return 1
     }
     if { $r == 1 } {
        expr {$n % $m}
     }

     if { $r%2 == 0 } {
        set nn [powm $n [expr {$r/2}] $m]
        expr { ($nn*$nn) % $m}
     } else {
        set nn [powm $n [expr {$r-1}] $m]
        expr { ($n*$nn) % $m}
     }
  }

  #
  # Testing the function
  #
  # The following cases are easy to verify manually:
  # powm(1,10,2) = 1
  # powm(2,10,3) = 1024 mod 3 = 1
  # powm(3,4,5) = 81 mod 5 = 1

  puts "powm(1,10,2) = [powm 1 10 2 ]"
  puts "powm(2,10,3) = [powm 2 10 3 ]"
  puts "powm(3, 4,5) = [powm 3  4 5 ]"

  #
  # Larger numbers (results respectively, 1, 71, 381, 12685)
  #
  puts "powm(31,24,15)          = [powm 31 24 15 ]"
  puts "powm(131,124,115)       = [powm 131 124 115 ]"
  puts "powm(2131,3124,4115)    = [powm 2131 3124 4115 ]"
  puts "powm(52131,53124,54115) = [powm 52131 53124 54115 ]"

----
Note that in the last test cases the intermediate numbers would be
much too big to deal with directly:

    52131^53124 = O(10^577007)

or, to put it in another way, a figure with half a million digits

----
[[
[Category Cryptography] |
[Category Mathematics]
]]