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