summaryrefslogtreecommitdiffstats
path: root/modular.scm
diff options
context:
space:
mode:
Diffstat (limited to 'modular.scm')
-rw-r--r--modular.scm180
1 files changed, 143 insertions, 37 deletions
diff --git a/modular.scm b/modular.scm
index a653739..e836100 100644
--- a/modular.scm
+++ b/modular.scm
@@ -1,5 +1,5 @@
;;;; "modular.scm", modular fixnum arithmetic for Scheme
-;;; Copyright (C) 1991, 1993, 1995 Aubrey Jaffer
+;;; 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
@@ -8,7 +8,7 @@
;1. Any copy made of this software must include this copyright notice
;in full.
;
-;2. I have made no warrantee or representation that the operation of
+;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.
;
@@ -17,45 +17,128 @@
;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)))
-(define (modular:normalize m k)
- (cond ((positive? m) (modulo k m))
- ((zero? m) k)
- ((<= m k (- m)) k)
- ((or (provided? 'bignum)
- (<= 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))))
+;;@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!
-(require 'logical)
+;;@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
-(define (modular: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))))
+;;@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
@@ -68,12 +151,17 @@
(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))
@@ -91,6 +179,8 @@
(+ 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))
@@ -102,7 +192,15 @@
;;; modular:r = 2**((nb-2)/2) where nb = number of bits in a word.
(define modular:r
- (ash 1 (quotient (integer-length most-positive-fixnum) 2)))
+ (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)
@@ -115,9 +213,8 @@
(cond
((zero? m) (* a b))
((negative? m)
- "This doesn't work for the full range of modulus M;"
- "Someone please create or convert the following"
- "algorighm to work with symmetric representation"
+ ;; 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
@@ -146,13 +243,22 @@
(modulo (+ (if (positive? p) (- p m) p)
(* a0 (modulo b q))) m)))))))))
-(define (modular:expt m a b)
- (cond ((= a 1) 1)
- ((= a (- m 1)) (if (odd? b) a 1))
- ((zero? a) 0)
- ((zero? m) (integer-expt a b))
- (else
- (logical:ipow-by-squaring a b 1
- (lambda (c d) (modular:* m c d))))))
-
-(define extended-euclid modular:extended-euclid)
+;;@args modulus n2 n3
+;;Returns (@var{n2} ^ @var{n3}) mod @var{modulus}.
+(define modular:expt
+ (let ((integer-expt (and (provided? 'inexact) expt)))
+ (lambda (m n xpn)
+ (cond ((= n 1) 1)
+ ((= n (- m 1)) (if (odd? xpn) n 1))
+ ((and (zero? m) integer-expt) (integer-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))))))))))