;;;; "modular.scm", modular fixnum arithmetic for Scheme
;;; Copyright (C) 1991, 1993, 1995, 2001, 2002 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, and to use it for any purpose is
;granted, subject to the following restrictions and understandings.
;
;1.  Any copy made of this software must include this copyright notice
;in full.
;
;2.  I have made no warranty or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3.  In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.

(require 'multiarg/and-)

;;@code{(require 'modular)}
;;@ftindex modular

;;@body
;;These procedures implement the Common-Lisp functions of the same names.
;;The real number @var{x2} must be non-zero.
;;@code{mod} returns @code{(- @var{x1} (* @var{x2} (floor (/ @var{x1} @var{x2}))))}.
;;@code{rem} returns @code{(- @var{x1} (* @var{x2} (truncate (/ @var{x1} @var{x2}))))}.
;;
;;If @var{x1} and @var{x2} are integers, then @code{mod} behaves like
;;@code{modulo} and @code{rem} behaves like @code{remainder}.
;;
;;@format
;;@t{(mod -90 360)                          @result{} 270
;;(rem -90 180)                          @result{} -90
;;
;;(mod 540 360)                          @result{} 180
;;(rem 540 360)                          @result{} 180
;;
;;(mod (* 5/2 pi) (* 2 pi))              @result{} 1.5707963267948965
;;(rem (* -5/2 pi) (* 2 pi))             @result{} -1.5707963267948965
;;}
;;@end format
(define (mod x1 x2)
  (if (and (integer? x1) (exact? x1) (integer? x2) (exact? x2))
      (modulo x1 x2)
      (- x1 (* x2 (floor (/ x1 x2))))))
(define (rem x1 x2)
  (if (and (integer? x1) (exact? x1) (integer? x2) (exact? x2))
      (remainder x1 x2)
      (- x1 (* x2 (truncate (/ x1 x2))))))

;;@args n1 n2
;;Returns a list of 3 integers @code{(d x y)} such that d = gcd(@var{n1},
;;@var{n2}) = @var{n1} * x + @var{n2} * y.
(define (extended-euclid x y)
  (define q 0)
  (do ((r0 x r1) (r1 y (remainder r0 r1))
       (u0 1 u1) (u1 0 (- u0 (* q u1)))
       (v0 0 v1) (v1 1 (- v0 (* q v1))))
      ;; (assert (= r0 (+ (* u0 x) (* v0 y))))
      ;; (assert (= r1 (+ (* u1 x) (* v1 y))))
      ((zero? r1) (list r0 u0 v0))
    (set! q (quotient r0 r1))))

(define modular:extended-euclid extended-euclid)

;;@body
;;Returns @code{(quotient (+ -1 n) -2)} for positive odd integer @var{n}.
(define (symmetric:modulus n)
  (cond ((or (not (number? n)) (not (positive? n)) (even? n))
	 (slib:error 'symmetric:modulus n))
	(else (quotient (+ -1 n) -2))))

;;@args modulus
;;Returns the non-negative integer characteristic of the ring formed when
;;@var{modulus} is used with @code{modular:} procedures.
(define (modulus->integer m)
  (cond ((negative? m) (- 1 m m))
	((zero? m) #f)
	(else m)))

;;@args modulus n
;;Returns the integer @code{(modulo @var{n} (modulus->integer
;;@var{modulus}))} in the representation specified by @var{modulus}.
(define modular:normalize
  (if (provided? 'bignum)
      (lambda (m k)
	(cond ((positive? m) (modulo k m))
	      ((zero? m) k)
	      ((<= m k (- m)) k)
	      (else
	       (let* ((pm (+ 1 (* -2 m)))
		      (s (modulo k pm)))
		 (if (<= s (- m)) s (- s pm))))))
      (lambda (m k)
	(cond ((positive? m) (modulo k m))
	      ((zero? m) k)
	      ((<= m k (- m)) k)
	      ((<= m (quotient (+ -1 most-positive-fixnum) 2))
	       (let* ((pm (+ 1 (* -2 m)))
		      (s (modulo k pm)))
		 (if (<= s (- m)) s (- s pm))))
	      ((positive? k) (+ (+ (+ k -1) m) m))
	      (else  (- (- (+ k 1) m) m))))))

;;;; NOTE: The rest of these functions assume normalized arguments!

;;@noindent
;;The rest of these functions assume normalized arguments; That is, the
;;arguments are constrained by the following table:
;;
;;@noindent
;;For all of these functions, if the first argument (@var{modulus}) is:
;;@table @code
;;@item positive?
;;Work as before.  The result is between 0 and @var{modulus}.
;;
;;@item zero?
;;The arguments are treated as integers.  An integer is returned.
;;
;;@item negative?
;;The arguments and result are treated as members of the integers modulo
;;@code{(+ 1 (* -2 @var{modulus}))}, but with @dfn{symmetric}
;;representation; i.e. @code{(<= (- @var{modulus}) @var{n}
;;@var{modulus})}.
;;@end table

;;@noindent
;;If all the arguments are fixnums the computation will use only fixnums.

;;@args modulus k
;;Returns @code{#t} if there exists an integer n such that @var{k} * n
;;@equiv{} 1 mod @var{modulus}, and @code{#f} otherwise.
(define (modular:invertable? m a)
  (eqv? 1 (gcd (or (modulus->integer m) 0) a)))

;;@args modulus n2
;;Returns an integer n such that 1 = (n * @var{n2}) mod @var{modulus}.  If
;;@var{n2} has no inverse mod @var{modulus} an error is signaled.
(define (modular:invert m a)
  (cond ((eqv? 1 (abs a)) a)		; unit
	(else
	 (let ((pm (modulus->integer m)))
	   (cond
	    (pm
	     (let ((d (modular:extended-euclid (modular:normalize pm a) pm)))
	       (if (= 1 (car d))
		   (modular:normalize m (cadr d))
		   (slib:error 'modular:invert "can't invert" m a))))
	    (else (slib:error 'modular:invert "can't invert" m a)))))))

;;@args modulus n2
;;Returns (@minus{}@var{n2}) mod @var{modulus}.
(define (modular:negate m a)
  (if (zero? a) 0
      (if (negative? m) (- a)
	  (- m a))))

;;; Being careful about overflow here

;;@args modulus n2 n3
;;Returns (@var{n2} + @var{n3}) mod @var{modulus}.
(define (modular:+ m a b)
  (cond ((positive? m)
	 (modulo (+ (- a m) b) m))
	((zero? m) (+ a b))
	((negative? a)
	 (if (negative? b)
	     (let ((s (+ (- a m) b)))
	       (if (negative? s)
		   (- s -1 m)
		   (+ s m)))
	     (+ a b)))
	((negative? b) (+ a b))
	(else (let ((s (+ (+ a m) b)))
		(if (positive? s)
		    (+ s -1 m)
		    (- s m))))))

;;@args modulus n2 n3
;;Returns (@var{n2} @minus{} @var{n3}) mod @var{modulus}.
(define (modular:- m a b)
  (cond ((positive? m) (modulo (- a b) m))
	((zero? m) (- a b))
	(else (modular:+ m a (- b)))))

;;; See: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
;;; with Splitting Facilities." ACM Transactions on Mathematical
;;; Software, 17:98-111 (1991)

;;; modular:r = 2**((nb-2)/2) where nb = number of bits in a word.
(define modular:r
  (do ((mpf most-positive-fixnum (quotient mpf 4))
       (r 1 (* 2 r)))
      ((<= mpf 0) (quotient r 2))))

;;@args modulus n2 n3
;;Returns (@var{n2} * @var{n3}) mod @var{modulus}.
;;
;;The Scheme code for @code{modular:*} with negative @var{modulus} is
;;not completed for fixnum-only implementations.
(define modular:*
  (if (provided? 'bignum)
      (lambda (m a b)
	(cond ((zero? m) (* a b))
	      ((positive? m) (modulo (* a b) m))
	      (else (modular:normalize m (* a b)))))
      (lambda (m a b)
	(let ((a0 a)
	      (p 0))
	  (cond
	   ((zero? m) (* a b))
	   ((negative? m)
	    ;; This doesn't work for the full range of modulus M.
	    ;; Need algorighm to work with symmetric representation.
	    (modular:normalize m (* a b)))
	   (else
	    (cond
	     ((< a modular:r))
	     ((< b modular:r) (set! a b) (set! b a0) (set! a0 a))
	     (else
	      (set! a0 (modulo a modular:r))
	      (let ((a1 (quotient a modular:r))
		    (qh (quotient m modular:r))
		    (rh (modulo m modular:r)))
		(cond ((>= a1 modular:r)
		       (set! a1 (- a1 modular:r))
		       (set! p (modulo (- (* modular:r (modulo b qh))
					  (* (quotient b qh) rh)) m))))
		(cond ((not (zero? a1))
		       (let ((q (quotient m a1)))
			 (set! p (- p (* (quotient b q) (modulo m a1))))
			 (set! p (modulo (+ (if (positive? p) (- p m) p)
					    (* a1 (modulo b q))) m)))))
		(set! p (modulo (- (* modular:r (modulo p qh))
				   (* (quotient p qh) rh)) m)))))
	    (if (zero? a0)
		p
		(let ((q (quotient m a0)))
		  (set! p (- p (* (quotient b q) (modulo m a0))))
		  (modulo (+ (if (positive? p) (- p m) p)
			     (* a0 (modulo b q))) m)))))))))

;;@args modulus n2 n3
;;Returns (@var{n2} ^ @var{n3}) mod @var{modulus}.
(define (modular:expt m n xpn)
  (cond ((= n 1) 1)
	((= n (- m 1)) (if (odd? xpn) n 1))
	((zero? m) (expt n xpn))
	((negative? xpn)
	 (modular:expt m (modular:invert m n) (- xpn)))
	((zero? n) 0)
	(else
	 (do ((x n (modular:* m x x))
	      (j xpn (quotient j 2))
	      (acc 1 (if (even? j) acc (modular:* m x acc))))
	     ((<= j 1)
	      (case j
		((0) acc)
		((1) (modular:* m x acc))))))))