home *** CD-ROM | disk | FTP | other *** search
/ Dream 44 / Amiga_Dream_44.iso / RiscPc / programmation / scm4e2.arc / !Scm / slib / modular < prev    next >
Text File  |  1994-05-25  |  3KB  |  91 lines
  1. ;;;; "modular.scm", modular fixnum arithmetic for Scheme
  2. ;;; Copyright (C) 1991, 1993 Aubrey Jaffer.
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;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 warrantee 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. (require 'logical)
  21.  
  22. (define (modular:extended-euclid x y)
  23.   (define q 0)
  24.   (do ((r0 x r1) (r1 y (remainder r0 r1))
  25.        (u0 1 u1) (u1 0 (- u0 (* q u1)))
  26.        (v0 0 v1) (v1 1 (- v0 (* q v1))))
  27.       ;; (assert (= r0 (+ (* u0 x) (* v0 y))))
  28.       ;; (assert (= r1 (+ (* u1 x) (* v1 y))))
  29.       ((zero? r1) (list r0 u0 v0))
  30.     (set! q (quotient r0 r1))))
  31.  
  32. (define (modular:invert m a)
  33.   (let ((d (modular:extended-euclid a m)))
  34.     (if (= 1 (car d))
  35.     (modulo (cadr d) m)
  36.     (slib:error "modular:invert can't invert" m a))))
  37.  
  38. (define (modular:negate m a) (if (zero? a) 0 (- m a)))
  39.  
  40. ;;; Being careful about overflow here
  41. (define (modular:+ m a b) (modulo (+ (- a m) b) m))
  42.  
  43. (define (modular:- m a b) (modulo (- a b) m))
  44.  
  45. ;;; See: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  46. ;;; with Splitting Facilities." ACM Transactions on Mathematical
  47. ;;; Software, 17:98-111 (1991)
  48.  
  49. ;;; modular:r = 2**((nb-2)/2) where nb = number of bits in a word.
  50. (define modular:r
  51.   (ash 1 (quotient (integer-length most-positive-fixnum) 2)))
  52. (define modular:*
  53.   (if (provided? 'bignum) (lambda (m a b) (modulo (* a b) m))
  54.       (lambda (m a b)
  55.     (let ((a0 a)
  56.           (p 0))
  57.       (cond ((< a modular:r))
  58.         ((< b modular:r) (set! a b) (set! b a0) (set! a0 a))
  59.         (else
  60.          (set! a0 (modulo a modular:r))
  61.          (let ((a1 (quotient a modular:r))
  62.                (qh (quotient m modular:r))
  63.                (rh (modulo m modular:r)))
  64.            (cond ((>= a1 modular:r)
  65.               (set! a1 (- a1 modular:r))
  66.               (set! p (modulo (- (* modular:r (modulo b qh))
  67.                          (* (quotient b qh) rh)) m))))
  68.            (cond ((not (zero? a1))
  69.               (let ((q (quotient m a1)))
  70.                 (set! p (- p (* (quotient b q) (modulo m a1))))
  71.                 (set! p (modulo (+ (if (positive? p) (- p m) p)
  72.                            (* a1 (modulo b q))) m)))))
  73.            (set! p (modulo (- (* modular:r (modulo p qh))
  74.                       (* (quotient p qh) rh)) m)))))
  75.       (if (zero? a0)
  76.           p
  77.           (let ((q (quotient m a0)))
  78.         (set! p (- p (* (quotient b q) (modulo m a0))))
  79.         (modulo (+ (if (positive? p) (- p m) p)
  80.                (* a0 (modulo b q))) m)))))))
  81.  
  82. (define (modular:expt m a b)
  83.   (cond ((= a 1) 1)
  84.     ((= a (- m 1)) (if (odd? b) a 1))
  85.     ((zero? a) 0)
  86.     (else
  87.      (logical:ipow-by-squaring a b 1
  88.                    (lambda (c d) (modular:* m c d))))))
  89.  
  90. (define extended-euclid modular:extended-euclid)
  91.