Import Upstream version 3a4upstream/3a4

1 files changed, 77 insertions, 101 deletions

diff --git a/modular.scm b/modular.scm
index 052bf92..e77ced4 100644
--- a/modular.scm
+++ b/modular.scm
@@ -1,5 +1,5 @@
 ;;;; "modular.scm", modular fixnum arithmetic for Scheme
-;;; Copyright (C) 1991, 1993, 1995, 2001, 2002 Aubrey Jaffer
+;;; Copyright (C) 1991, 1993, 1995, 2001, 2002, 2006 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
@@ -17,40 +17,9 @@
 ;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.
@@ -59,30 +28,33 @@
   (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))))
+;;For odd positive integer @1, returns an object suitable for passing
+;;as the first argument to @code{modular:} procedures, directing them
+;;to return a symmetric modular number, ie. an @var{n} such that
+;;@example
+;;(<= (quotient @var{m} -2) @var{n} (quotient @var{m} 2)
+;;@end example
+(define (symmetric:modulus m)
+  (cond ((or (not (number? m)) (not (positive? m)) (even? m))
+	 (slib:error 'symmetric:modulus m))
+	(else (quotient (+ -1 m) -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))
+(define (modular:characteristic m)
+  (cond ((negative? m) (- 1 (+ m m)))
 	((zero? m) #f)
 	(else m)))
 
 ;;@args modulus n
-;;Returns the integer @code{(modulo @var{n} (modulus->integer
+;;Returns the integer @code{(modulo @var{n} (modular:characteristic
 ;;@var{modulus}))} in the representation specified by @var{modulus}.
 (define modular:normalize
   (if (provided? 'bignum)
@@ -115,17 +87,21 @@
 ;;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}.
+;;Integers mod @var{modulus}.  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
+;;Otherwise, if @var{modulus} is a value returned by
+;;@code{(symmetric:modulus @var{radix})}, then the arguments and
+;;result are treated as members of the integers modulo @var{radix},
+;;but with @dfn{symmetric} representation; i.e.
+;;@example
+;;(<= (quotient @var{radix} 2) @var{n} (quotient (- -1 @var{radix}) 2)
+;;@end example
 
 ;;@noindent
 ;;If all the arguments are fixnums the computation will use only fixnums.
@@ -134,43 +110,44 @@
 ;;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)))
+  (eqv? 1 (gcd (or (modular:characteristic 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)
+  (define (barf) (slib:error 'modular:invert "can't invert" m a))
   (cond ((eqv? 1 (abs a)) a)		; unit
 	(else
-	 (let ((pm (modulus->integer m)))
+	 (let ((pm (modular:characteristic 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)))))))
+		   (barf))))
+	    (else (barf)))))))
 
 ;;@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))))
+  (cond ((zero? a) 0)
+	((negative? m) (- a))
+	(else (- 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))
+  (cond ((positive? m) (modulo (+ (- a m) b) m))
 	((zero? m) (+ a b))
+	;; m is negative
 	((negative? a)
 	 (if (negative? b)
 	     (let ((s (+ (- a m) b)))
 	       (if (negative? s)
-		   (- s -1 m)
+		   (- s (+ -1 m))
 		   (+ s m)))
 	     (+ a b)))
 	((negative? b) (+ a b))
@@ -182,9 +159,7 @@
 ;;@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)))))
+  (modular:+ m a (modular:negate m b)))
 
 ;;; See: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 ;;; with Splitting Facilities." ACM Transactions on Mathematical
@@ -208,52 +183,53 @@
 	      ((positive? m) (modulo (* a b) m))
 	      (else (modular:normalize m (* a b)))))
       (lambda (m a b)
-	(let ((a0 a)
-	      (p 0))
+	(define a0 a)
+	(define p 0)
+	
+	(cond
+	 ((zero? m) (* a b))
+	 ((negative? m)
+	  ;; Need algorighm to work with symmetric representation.
+	  (modular:normalize m (* a b)))
+	 (else
 	  (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)))
+	   ((< a modular:r))
+	   ((< b modular:r) (set! a b) (set! b a0) (set! a0 a))
 	   (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)))))))))
+	    (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))
+(define (modular:expt m base xpn)
+  (cond ((zero? m) (expt base xpn))
+	((= base 1) 1)
+	((if (negative? m) (= -1 base) (= (- m 1) base))
+	 (if (odd? xpn) base 1))
 	((negative? xpn)
-	 (modular:expt m (modular:invert m n) (- xpn)))
-	((zero? n) 0)
+	 (modular:expt m (modular:invert m base) (- xpn)))
+	((zero? base) 0)
 	(else
-	 (do ((x n (modular:* m x x))
+	 (do ((x base (modular:* m x x))
 	      (j xpn (quotient j 2))
 	      (acc 1 (if (even? j) acc (modular:* m x acc))))
 	     ((<= j 1)