1;;;; "modular.scm", modular fixnum arithmetic for Scheme 2;;; Copyright (C) 1991, 1993, 1995, 2001, 2002, 2006 Aubrey Jaffer 3; 4;Permission to copy this software, to modify it, to redistribute it, 5;to distribute modified versions, and to use it for any purpose is 6;granted, subject to the following restrictions and understandings. 7; 8;1. Any copy made of this software must include this copyright notice 9;in full. 10; 11;2. I have made no warranty or representation that the operation of 12;this software will be error-free, and I am under no obligation to 13;provide any services, by way of maintenance, update, or otherwise. 14; 15;3. In conjunction with products arising from the use of this 16;material, there shall be no use of my name in any advertising, 17;promotional, or sales literature without prior written consent in 18;each case. 19 20;;@code{(require 'modular)} 21;;@ftindex modular 22 23;;@args n1 n2 24;;Returns a list of 3 integers @code{(d x y)} such that d = gcd(@var{n1}, 25;;@var{n2}) = @var{n1} * x + @var{n2} * y. 26(define (extended-euclid x y) 27 (define q 0) 28 (do ((r0 x r1) (r1 y (remainder r0 r1)) 29 (u0 1 u1) (u1 0 (- u0 (* q u1))) 30 (v0 0 v1) (v1 1 (- v0 (* q v1)))) 31 ((zero? r1) (list r0 u0 v0)) 32 (set! q (quotient r0 r1)))) 33 34(define modular:extended-euclid extended-euclid) 35 36;;@body 37;;For odd positive integer @1, returns an object suitable for passing 38;;as the first argument to @code{modular:} procedures, directing them 39;;to return a symmetric modular number, ie. an @var{n} such that 40;;@example 41;;(<= (quotient @var{m} -2) @var{n} (quotient @var{m} 2) 42;;@end example 43(define (symmetric:modulus m) 44 (cond ((or (not (number? m)) (not (positive? m)) (even? m)) 45 (slib:error 'symmetric:modulus m)) 46 (else (quotient (+ -1 m) -2)))) 47 48;;@args modulus 49;;Returns the non-negative integer characteristic of the ring formed when 50;;@var{modulus} is used with @code{modular:} procedures. 51(define (modular:characteristic m) 52 (cond ((negative? m) (- 1 (+ m m))) 53 ((zero? m) #f) 54 (else m))) 55 56;;@args modulus n 57;;Returns the integer @code{(modulo @var{n} (modular:characteristic 58;;@var{modulus}))} in the representation specified by @var{modulus}. 59(define modular:normalize 60 (if (provided? 'bignum) 61 (lambda (m k) 62 (cond ((positive? m) (modulo k m)) 63 ((zero? m) k) 64 ((<= m k (- m)) k) 65 (else 66 (let* ((pm (+ 1 (* -2 m))) 67 (s (modulo k pm))) 68 (if (<= s (- m)) s (- s pm)))))) 69 (lambda (m k) 70 (cond ((positive? m) (modulo k m)) 71 ((zero? m) k) 72 ((<= m k (- m)) k) 73 ((<= m (quotient (+ -1 most-positive-fixnum) 2)) 74 (let* ((pm (+ 1 (* -2 m))) 75 (s (modulo k pm))) 76 (if (<= s (- m)) s (- s pm)))) 77 ((positive? k) (+ (+ (+ k -1) m) m)) 78 (else (- (- (+ k 1) m) m)))))) 79 80;;;; NOTE: The rest of these functions assume normalized arguments! 81 82;;@noindent 83;;The rest of these functions assume normalized arguments; That is, the 84;;arguments are constrained by the following table: 85;; 86;;@noindent 87;;For all of these functions, if the first argument (@var{modulus}) is: 88;;@table @code 89;;@item positive? 90;;Integers mod @var{modulus}. The result is between 0 and 91;;@var{modulus}. 92;; 93;;@item zero? 94;;The arguments are treated as integers. An integer is returned. 95;;@end table 96;; 97;;@noindent 98;;Otherwise, if @var{modulus} is a value returned by 99;;@code{(symmetric:modulus @var{radix})}, then the arguments and 100;;result are treated as members of the integers modulo @var{radix}, 101;;but with @dfn{symmetric} representation; i.e. 102;;@example 103;;(<= (quotient @var{radix} 2) @var{n} (quotient (- -1 @var{radix}) 2) 104;;@end example 105 106;;@noindent 107;;If all the arguments are fixnums the computation will use only fixnums. 108 109;;@args modulus k 110;;Returns @code{#t} if there exists an integer n such that @var{k} * n 111;;@equiv{} 1 mod @var{modulus}, and @code{#f} otherwise. 112(define (modular:invertable? m a) 113 (eqv? 1 (gcd (or (modular:characteristic m) 0) a))) 114 115;;@args modulus n2 116;;Returns an integer n such that 1 = (n * @var{n2}) mod @var{modulus}. If 117;;@var{n2} has no inverse mod @var{modulus} an error is signaled. 118(define (modular:invert m a) 119 (define (barf) (slib:error 'modular:invert "can't invert" m a)) 120 (cond ((eqv? 1 (abs a)) a) ; unit 121 (else 122 (let ((pm (modular:characteristic m))) 123 (cond 124 (pm 125 (let ((d (modular:extended-euclid (modular:normalize pm a) pm))) 126 (if (= 1 (car d)) 127 (modular:normalize m (cadr d)) 128 (barf)))) 129 (else (barf))))))) 130 131;;@args modulus n2 132;;Returns (@minus{}@var{n2}) mod @var{modulus}. 133(define (modular:negate m a) 134 (cond ((zero? a) 0) 135 ((negative? m) (- a)) 136 (else (- m a)))) 137 138;;; Being careful about overflow here 139 140;;@args modulus n2 n3 141;;Returns (@var{n2} + @var{n3}) mod @var{modulus}. 142(define (modular:+ m a b) 143 (cond ((positive? m) (modulo (+ (- a m) b) m)) 144 ((zero? m) (+ a b)) 145 ;; m is negative 146 ((negative? a) 147 (if (negative? b) 148 (let ((s (+ (- a m) b))) 149 (if (negative? s) 150 (- s (+ -1 m)) 151 (+ s m))) 152 (+ a b))) 153 ((negative? b) (+ a b)) 154 (else (let ((s (+ (+ a m) b))) 155 (if (positive? s) 156 (+ s -1 m) 157 (- s m)))))) 158 159;;@args modulus n2 n3 160;;Returns (@var{n2} @minus{} @var{n3}) mod @var{modulus}. 161(define (modular:- m a b) 162 (modular:+ m a (modular:negate m b))) 163 164;;; See: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package 165;;; with Splitting Facilities." ACM Transactions on Mathematical 166;;; Software, 17:98-111 (1991) 167 168;;; modular:r = 2**((nb-2)/2) where nb = number of bits in a word. 169(define modular:r 170 (do ((mpf most-positive-fixnum (quotient mpf 4)) 171 (r 1 (* 2 r))) 172 ((<= mpf 0) (quotient r 2)))) 173 174;;@args modulus n2 n3 175;;Returns (@var{n2} * @var{n3}) mod @var{modulus}. 176;; 177;;The Scheme code for @code{modular:*} with negative @var{modulus} is 178;;not completed for fixnum-only implementations. 179(define modular:* 180 (if (provided? 'bignum) 181 (lambda (m a b) 182 (cond ((zero? m) (* a b)) 183 ((positive? m) (modulo (* a b) m)) 184 (else (modular:normalize m (* a b))))) 185 (lambda (m a b) 186 (define a0 a) 187 (define p 0) 188 (cond 189 ((zero? m) (* a b)) 190 ((negative? m) 191 ;; Need algorighm to work with symmetric representation. 192 (modular:normalize m (* (modular:normalize m a) 193 (modular:normalize m b)))) 194 (else 195 (set! a (modulo a m)) 196 (set! b (modulo b m)) 197 (set! a0 a) 198 (cond ((< a modular:r)) 199 ((< b modular:r) (set! a b) (set! b a0) (set! a0 a)) 200 (else 201 (set! a0 (modulo a modular:r)) 202 (let ((a1 (quotient a modular:r)) 203 (qh (quotient m modular:r)) 204 (rh (modulo m modular:r))) 205 (cond ((>= a1 modular:r) 206 (set! a1 (- a1 modular:r)) 207 (set! p (modulo (- (* modular:r (modulo b qh)) 208 (* (quotient b qh) rh)) m)))) 209 (cond ((not (zero? a1)) 210 (let ((q (quotient m a1))) 211 (set! p (- p (* (quotient b q) (modulo m a1)))) 212 (set! p (modulo (+ (if (positive? p) (- p m) p) 213 (* a1 (modulo b q))) m))))) 214 (set! p (modulo (- (* modular:r (modulo p qh)) 215 (* (quotient p qh) rh)) m))))) 216 (if (zero? a0) 217 p 218 (let ((q (quotient m a0))) 219 (set! p (- p (* (quotient b q) (modulo m a0)))) 220 (modulo (+ (if (positive? p) (- p m) p) 221 (* a0 (modulo b q))) m)))))))) 222 223;;@args modulus n2 n3 224;;Returns (@var{n2} ^ @var{n3}) mod @var{modulus}. 225(define (modular:expt m base xpn) 226 (cond ((zero? m) (expt base xpn)) 227 ((= base 1) 1) 228 ((if (negative? m) (= -1 base) (= (- m 1) base)) 229 (if (odd? xpn) base 1)) 230 ((negative? xpn) 231 (modular:expt m (modular:invert m base) (- xpn))) 232 ((zero? base) 0) 233 (else 234 (do ((x base (modular:* m x x)) 235 (j xpn (quotient j 2)) 236 (acc 1 (if (even? j) acc (modular:* m x acc)))) 237 ((<= j 1) 238 (case j 239 ((0) acc) 240 ((1) (modular:* m x acc)))))))) 241