; Computing (a^b) mod n, assuming a, b, n are integers, and that
; b >= 0.
;
;
; First, the legalese:
;
; Author: David Neto
; Copyright (c) 1996, David Neto.
; 
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License
; as published by the Free Software Foundation; either version 2
; of the License, or (at your option) any later version.
; 
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
; 
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
; 
; The GNU General Public License is available at
; http://www.gnu.ai.mit.edu/copyleft/gnu.html 



; This is code for expmod in Ideal Turing, i.e. Turing with support
; for arthmetic over arbitrarily large numbers.  
; Note that all values remain bounded from above by n^2. 
; So addition or subtraction takes O(log n) time.  
; (This is linear in the size of the representation of the numbers.)
; Each multiply takes O((log n)^2) time.
; Assuming a<= n and b<= n, the total running time is
; is O((log n)^3) because at most (log_2 n)+1 iterations are performed.
;
; Invariant:
;  (a^(x*2^i) * accum) mod n = (a^b mod n)
;  and  y = (a^(2^i) mod n)
;
;   function expmod(a :integer,b:integer,n:integer) : integer
;       x := b mod n
;       y := a mod n
;       accum := 1
;       while x > 0
;           if (x mod 2) = 1 then
;               accum := (accum * y) mod n
;           end if
;           y := (y*y) mod n
;           x := x div 2
;       end while
;       result accum
;   end expmod



; Here is the code in Scheme.  Scheme has built-in support for arithmetic
; on large integers.  :)
;
; (expmod a b n) computes (a^b) mod n.
; The expmod function uses a helper function expmod-h.  It passes the right
; values to satisfy the variant initially.
;
(define (expmod a b n) 
    (expmod-h (remainder a n)   ; initial value of a.
              (remainder b n)   ; initial value of b.
              (remainder b n)   ; initial value of x.
              (remainder a n)   ; initial value of y.
              1                 ; initial value of accum
              0                 ; initial value of i.
              n ))
                

(define (expmod-h a b x y accum i n)
; invariant:
;  (a^(x*2^i) * accum) mod n = (a^b mod n)
;  and  y = a^(2^i)
; The code goes from iteration number i=0 and increases i by 1, satisfying
; the invariant.
; Note that in a more careful implementation, we would remove arguments a, b,
; and i because they aren't used in this code.  They are only used to state
; the invariant.
; Uncomment the following lines to see how a changes and how x quickly
; decays to 0.  
;(display a)(display " ")
;(display x)(newline)
    (if (= x 0)
        accum
        (expmod-h   ; Tail-recursive call.
            a 
            b 
            (/ (if (odd x) (- x 1) x) 2) 
                ; shift x right one bit, losing bottom bit
            (remainder (* y y) n)   ; square y
                ; The next argument is the accumulator.  
                ; It is multiplied by a^(2^i) if there is a 1 bit in position
                ; i of the binary expansion of b.
                ; Otherwise, accumulator stays the same.
            (if (odd x) (remainder (* accum y) n) accum)
            (+ i 1)
            n )))


; We've used the following function that answers the question whether 
; a given number is odd or not.
(define (odd x) (= 1 (remainder x 2)))