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