From 8ffbc2df0fde83082610149d24e594c1cd879f4a Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Mon, 20 Feb 2017 00:05:25 -0800 Subject: Import Upstream version 2a6 --- randinex.scm | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 randinex.scm (limited to 'randinex.scm') diff --git a/randinex.scm b/randinex.scm new file mode 100644 index 0000000..1c2b702 --- /dev/null +++ b/randinex.scm @@ -0,0 +1,99 @@ +;;;"randinex.scm" Pseudo-Random inexact real numbers for scheme. +;;; Copyright (C) 1991, 1993 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. + +;This file is loaded by random.scm if inexact numbers are supported by +;the implementation. + +;;; Fixed sphere and normal functions from: Harald Hanche-Olsen + +(define random:float-radix + (+ 1 (exact->inexact random:MASK))) + +;;; This determines how many chunks will be neccessary to completely +;;; fill up an inexact real. +(define (random:size-float l x) + (cond ((= 1.0 (+ 1 x)) l) + ((= 4 l) l) + (else (random:size-float (+ l 1) (/ x random:float-radix))))) +(define random:chunks/float (random:size-float 0 1.0)) + +(define (random:uniform-chunk n state) + (if (= 1 n) + (/ (exact->inexact (random:chunk state)) + random:float-radix) + (/ (+ (random:uniform-chunk (- n 1) state) + (exact->inexact (random:chunk state))) + random:float-radix))) + +;;; Generate an inexact real between 0 and 1. +(define (random:uniform state) + (random:uniform-chunk random:chunks/float state)) + +;;; If x and y are independent standard normal variables, then with +;;; x=r*cos(t), y=r*sin(t), we find that t is uniformly distributed +;;; over [0,2*pi] and the cumulative distribution of r is +;;; 1-exp(-r^2/2). This latter means that u=exp(-r^2/2) is uniformly +;;; distributed on [0,1], so r=sqrt(-2 log u) can be used to generate r. + +(define (random:normal-vector! vect . args) + (let ((state (if (null? args) *random-state* (car args))) + (sum2 0)) + (let ((do! (lambda (k x) + (vector-set! vect k x) + (set! sum2 (+ sum2 (* x x)))))) + (do ((n (- (vector-length vect) 1) (- n 2))) + ((negative? n) sum2) + (let ((t (* 6.28318530717958 (random:uniform state))) + (r (sqrt (* -2 (log (random:uniform state)))))) + (do! n (* r (cos t))) + (if (positive? n) (do! (- n 1) (* r (sin t))))))))) + +(define random:normal + (let ((vect (make-vector 1))) + (lambda args + (apply random:normal-vector! vect args) + (vector-ref vect 0)))) + +;;; For the uniform distibution on the hollow sphere, pick a normal +;;; family and scale. + +(define (random:hollow-sphere! vect . args) + (let ((ms (sqrt (apply random:normal-vector! vect args)))) + (do ((n (- (vector-length vect) 1) (- n 1))) + ((negative? n)) + (vector-set! vect n (/ (vector-ref vect n) ms))))) + +;;; For the uniform distribution on the solid sphere, note that in +;;; this distribution the length r of the vector has cumulative +;;; distribution r^n; i.e., u=r^n is uniform [0,1], so r kan be +;;; generated as r=u^(1/n). + +(define (random:solid-sphere! vect . args) + (apply random:hollow-sphere! vect args) + (let ((r (expt (random:uniform (if (null? args) *random-state* (car args))) + (/ (vector-length vect))))) + (do ((n (- (vector-length vect) 1) (- n 1))) + ((negative? n)) + (vector-set! vect n (* r (vector-ref vect n)))))) + +(define (random:exp . args) + (let ((state (if (null? args) *random-state* (car args)))) + (- (log (random:uniform state))))) + +(require 'random) -- cgit v1.2.3