summaryrefslogtreecommitdiffstats
path: root/modular.scm
diff options
context:
space:
mode:
Diffstat (limited to 'modular.scm')
-rw-r--r--modular.scm158
1 files changed, 158 insertions, 0 deletions
diff --git a/modular.scm b/modular.scm
new file mode 100644
index 0000000..357ce77
--- /dev/null
+++ b/modular.scm
@@ -0,0 +1,158 @@
+;;;; "modular.scm", modular fixnum arithmetic for Scheme
+;;; Copyright (C) 1991, 1993, 1995 Aubrey Jaffer.
+;
+;Permission to copy this software, to redistribute it, and to use it
+;for any purpose is granted, subject to the following restrictions and
+;understandings.
+;
+;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
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3. In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+(define (symmetric:modulus n)
+ (cond ((or (not (number? n)) (not (positive? n)) (even? n))
+ (slib:error 'symmetric:modulus n))
+ (else (quotient (+ -1 n) -2))))
+
+(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))))
+
+;;;; NOTE: The rest of these functions assume normalized arguments!
+
+(require 'logical)
+
+(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))))
+
+(define (modular:invertable? m a)
+ (eqv? 1 (gcd (or (modulus->integer m) 0) a)))
+
+(define (modular:invert m a)
+ (cond ((eqv? 1 (abs a)) a) ; unit
+ (else
+ (let ((pm (modulus->integer 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)))))))
+
+(define (modular:negate m a)
+ (if (zero? a) 0
+ (if (negative? m) (- a)
+ (- m a))))
+
+;;; Being careful about overflow here
+(define (modular:+ m a b)
+ (cond ((positive? m)
+ (modulo (+ (- a m) b) m))
+ ((zero? m) (+ a b))
+ ((negative? a)
+ (if (negative? b)
+ (let ((s (+ (- a m) b)))
+ (if (negative? s)
+ (- s -1 m)
+ (+ s m)))
+ (+ a b)))
+ ((negative? b) (+ a b))
+ (else (let ((s (+ (+ a m) b)))
+ (if (positive? s)
+ (+ s -1 m)
+ (- s m))))))
+
+(define (modular:- m a b)
+ (cond ((positive? m) (modulo (- a b) m))
+ ((zero? m) (- a b))
+ (else (modular:+ m a (- b)))))
+
+;;; See: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
+;;; with Splitting Facilities." ACM Transactions on Mathematical
+;;; Software, 17:98-111 (1991)
+
+;;; 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)))
+(define modular:*
+ (if (provided? 'bignum)
+ (lambda (m a b)
+ (cond ((zero? m) (* a b))
+ ((positive? m) (modulo (* a b) m))
+ (else (modular:normalize m (* a b)))))
+ (lambda (m a b)
+ (let ((a0 a)
+ (p 0))
+ (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"
+ (modular:normalize m (* a b)))
+ (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)))))))))
+
+(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)