;;; Copyright (c) 2015 Taylor R. Campbell
;;; All rights reserved.
;;;
;;; Redistribution and use in source and binary forms, with or without
;;; modification, are permitted provided that the following conditions
;;; are met:
;;; 1. Redistributions of source code must retain the above copyright
;;;    notice, this list of conditions and the following disclaimer.
;;; 2. Redistributions in binary form must reproduce the above copyright
;;;    notice, this list of conditions and the following disclaimer in the
;;;    documentation and/or other materials provided with the distribution.
;;;
;;; THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
;;; ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
;;; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
;;; ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
;;; FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
;;; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
;;; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
;;; HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
;;; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
;;; OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
;;; SUCH DAMAGE.

;;; Fragments related to modular inverses.  See
;;; <http://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two>
;;; for explanation.

(define (euclid a b)
  (define (loop s s* t t* r r*)
    ;; a s + b t - r = 0
    (if (zero? r*)
        (values s t r)
        (let ((q (quotient r r*))
              (r** (remainder r r*)))
          (loop s* (- s (* q s*))
                t* (- t (* q t*))
                r* r**))))
  (if (< a b)
      (loop 0 1 1 0 b a)
      (loop 1 0 0 1 a b)))

;; Bezout's identity:
(let ((a 3) (b (expt 2 32)))
  (receive (s t r) (euclid a b)
    (- (+ (* a s) (* b t)) r)))
;Value: 0

;; Modular inverse:
(receive (s t r) (euclid 3 (expt 2 32))
  (modulo s (expt 2 32)))
;Value: 2863311531

(modulo (* 3 2863311531) (expt 2 32))
;Value: 1

(define (modexp a e n)
  (if (zero? e)
      1
      (let loop ((r 1) (b (modulo a n)) (f e))
        ;; r b^f = a^e (mod n), f > 0
        (let ((r (if (even? f) r (modulo (* r b) n)))
              (f (quotient f 2)))
          (if (zero? f)
              r
              (loop r (modulo (* b b) n) f))))))

(modexp 3 (- (expt 2 31) 1) (expt 2 32))
;Value: 2863311531

(define (modinv2^2^i a i)
  (let loop ((r 0) (b a))
    ;; a b = 1 (mod 2^(2^r))
    (if (< r i)
        (loop (+ r 1) (* b (- 2 (* a b))))
        b)))

(modulo (modinv2^2^i 3 5) (expt 2 32))
;Value: 2863311531