diff options
Diffstat (limited to 'modular.scm')
| -rw-r--r-- | modular.scm | 178 |
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) |