summaryrefslogtreecommitdiffstats
path: root/modular.scm
diff options
context:
space:
mode:
Diffstat (limited to 'modular.scm')
-rw-r--r--modular.scm178
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)