home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / lisp / guile-ii.src / guile-ii / guile-src / gls / randinex.scm < prev    next >
Encoding:
Text File  |  1995-08-07  |  3.5 KB  |  100 lines
  1. ;;;"randinex.scm" Pseudo-Random inexact real numbers 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. ;This file is loaded by random.scm if inexact numbers are supported by
  21. ;the implementation.
  22.  
  23. ;;; Fixed sphere and normal functions from: Harald Hanche-Olsen
  24.  
  25. (define random:float-radix
  26.   (+ 1 (exact->inexact random:MASK)))
  27.  
  28. ;;; This determines how many chunks will be neccessary to completely
  29. ;;; fill up an inexact real.
  30. (define (random:size-float l x)
  31.   (cond ((= 1.0 (+ 1 x)) l)
  32.     ((= 4 l) l)
  33.     (else (random:size-float (+ l 1) (/ x random:float-radix)))))
  34. (define random:chunks/float (random:size-float 0 1.0))
  35.  
  36. (define (random:uniform-chunk n state)
  37.   (if (= 1 n)
  38.       (/ (exact->inexact (random:chunk state))
  39.      random:float-radix)
  40.       (/ (+ (random:uniform-chunk (- n 1) state)
  41.         (exact->inexact (random:chunk state)))
  42.      random:float-radix)))
  43.  
  44. ;;; Generate an inexact real between 0 and 1.
  45. (define (random:uniform state)
  46.   (random:uniform-chunk random:chunks/float state))
  47.  
  48. ;;; If x and y are independent standard normal variables, then with
  49. ;;; x=r*cos(t), y=r*sin(t), we find that t is uniformly distributed
  50. ;;; over [0,2*pi] and the cumulative distribution of r is
  51. ;;; 1-exp(-r^2/2).  This latter means that u=exp(-r^2/2) is uniformly
  52. ;;; distributed on [0,1], so r=sqrt(-2 log u) can be used to generate r.
  53.  
  54. (define (random:normal-vector! vect . args)
  55.   (let ((state (if (null? args) *random-state* (car args)))
  56.     (sum2 0))
  57.     (let ((do! (lambda (k x)
  58.          (vector-set! vect k x)
  59.          (set! sum2 (+ sum2 (* x x))))))
  60.       (do ((n (- (vector-length vect) 1) (- n 2)))
  61.       ((negative? n) sum2)
  62.     (let ((t (* 6.28318530717958 (random:uniform state)))
  63.           (r (sqrt (* -2 (log (random:uniform state))))))
  64.       (do! n (* r (cos t)))
  65.       (if (positive? n) (do! (- n 1) (* r (sin t)))))))))
  66.  
  67. (define random:normal
  68.   (let ((vect (make-vector 1)))
  69.     (lambda args 
  70.       (apply random:normal-vector! vect args)
  71.       (vector-ref vect 0))))
  72.  
  73. ;;; For the uniform distibution on the hollow sphere, pick a normal
  74. ;;; family and scale.
  75.  
  76. (define (random:hollow-sphere! vect . args)
  77.   (let ((ms (sqrt (apply random:normal-vector! vect args))))
  78.     (do ((n (- (vector-length vect) 1) (- n 1)))
  79.     ((negative? n))
  80.       (vector-set! vect n (/ (vector-ref vect n) ms)))))
  81.  
  82. ;;; For the uniform distribution on the solid sphere, note that in
  83. ;;; this distribution the length r of the vector has cumulative
  84. ;;; distribution r^n; i.e., u=r^n is uniform [0,1], so r kan be
  85. ;;; generated as r=u^(1/n).
  86.  
  87. (define (random:solid-sphere! vect . args)
  88.   (apply random:hollow-sphere! vect args)
  89.   (let ((r (expt (random:uniform (if (null? args) *random-state* (car args)))
  90.          (/ (vector-length vect)))))
  91.     (do ((n (- (vector-length vect) 1) (- n 1)))
  92.     ((negative? n))
  93.       (vector-set! vect n (* r (vector-ref vect n))))))
  94.  
  95. (define (random:exp . args)
  96.   (let ((state (if (null? args) *random-state* (car args))))
  97.     (- (log (random:uniform state)))))
  98.  
  99. (require 'random)
  100.