1)abbrev package RANDSRC RandomNumberSource 2++ Author: S.M.Watt 3++ Date Created: April 87 4++ Basic Operations: 5++ Related Domains: 6++ Also See: 7++ AMS Classifications: 8++ Keywords: 9++ Examples: 10++ References: 11++ Description: Random number generators 12--% RandomNumberSource 13++ All random numbers used in the system should originate from 14++ the same generator. This package is intended to be the source. 15-- 16-- Possible improvements: 17-- 1) Start where the user left off 18-- 2) Be able to switch between methods in the random number source. 19-- Currently unused. 20RandomNumberSource() : with 21 -- If r := randnum() then 0 <= r < size(). 22 randnum : () -> Integer 23 ++ randnum() is a random number between 0 and size(). 24 -- If r := randnum() then 0 <= r < size(). 25 size : () -> Integer 26 ++ size() is the base of the random number generator 27 28 -- If r := randnum n and n <= size() then 0 <= r < n. 29 randnum : Integer -> Integer 30 ++ randnum(n) is a random number between 0 and n-1. 31 reseed : Integer -> Void 32 ++ reseed(n) restarts the random number generator at n. 33 seed : () -> Integer 34 ++ seed() returns the current seed value. 35 36 == add 37 -- This random number generator passes the spectral test 38 -- with flying colours. [Knuth vol2, 2nd ed, p105] 39 ranbase : Integer := 2^31-1 40 x0 : Integer := 1231231231 41 x1 : Integer := 3243232987 42 43 randnum() == 44 t := (271828183 * x1 - 314159269 * x0) rem ranbase 45 if t < 0 then t := t + ranbase 46 x0:= x1 47 x1:= t 48 49 size() == ranbase 50 reseed n == 51 x0 := n rem ranbase 52 -- x1 := (n quo ranbase) rem ranbase 53 x1 := n quo ranbase 54 55 seed() == x1*ranbase + x0 56 57 -- Compute an integer in 0..n-1. 58 randnum n == 59 (n * randnum()) quo ranbase 60 61)abbrev package RDIST RandomDistributions 62++ Description: 63++ This package exports random distributions 64-- Currently unused. 65RandomDistributions(S : SetCategory) : with 66 uniform : Set S -> (() -> S) 67 ++ uniform(s) \undocumented 68 weighted : List Record(value : S, weight : Integer) -> (()->S) 69 ++ weighted(l) \undocumented 70 rdHack1 : (Vector S, Vector Integer, Integer)->(()->S) 71 ++ rdHack1(v, u, n) \undocumented 72 == add 73 import from RandomNumberSource() 74 75 weighted lvw == 76 -- Collapse duplicates, adding weights. 77 t : Table(S, Integer) := table() 78 for r in lvw repeat 79 u := search(r.value, t) 80 w := (u case "failed" => 0; u::Integer) 81 t r.value := w + r.weight 82 83 -- Construct vectors of values and cumulative weights. 84 kl := keys t 85 n := (#kl)::NonNegativeInteger 86 n = 0 => error "Cannot select from empty set" 87 kv : Vector(S) := new(n, kl.0) 88 wv : Vector(Integer) := new(n, 0) 89 90 totwt : Integer := 0 91 for k in kl for i in 1..n repeat 92 kv.i := k 93 totwt:= totwt + t k 94 wv.i := totwt 95 96 -- Function to generate an integer and lookup. 97 rdHack1(kv, wv, totwt) 98 99 rdHack1(kv, wv, totwt) == 100 w := randnum totwt 101 -- do binary search in wv 102 kv.1 103 104 uniform fset == 105 l := members fset 106 n := #l 107 l.(randnum(n)+1) 108 109)abbrev package INTBIT IntegerBits 110----> Bug! Cannot precompute params and return a function which 111----> simply computes the last call. e.g. ridHack1, below. 112 113--% IntegerBits 114-- Functions related to the binary representation of integers. 115-- These functions directly access the bits in the big integer 116-- representation and so are much facter than using a quotient loop. 117-- SMW Sept 86. 118 119 120++ Description: 121++ This package provides functions to lookup bits in integers. 122-- Currently unused. 123IntegerBits : with 124 bitCoef : (Integer, NonNegativeInteger) -> Integer 125 ++ bitCoef(n, m) returns the coefficient of 2^m in two 126 ++ complement representation of n. 127 bitTruth : (Integer, NonNegativeInteger) -> Boolean 128 ++ bitTruth(n, m) returns true if coefficient of 2^m in 129 ++ two complement representation of n is 1. 130 131 == add 132 133 bitTruth(n, i) == LOGBITP(i, n)$Lisp 134 bitCoef(n, i) == if bitTruth(n, i) then 1 else 0 135 136)abbrev package RIDIST RandomIntegerDistributions 137++ Description: 138++ This package exports integer distributions 139RationalNumber==> Fraction Integer 140-- Currently unused. 141RandomIntegerDistributions() : with 142 uniform : Segment Integer -> (() -> Integer) 143 ++ uniform(s) \undocumented 144 binomial : (Integer, RationalNumber) -> (() -> Integer) 145 ++ binomial(n, f) \undocumented 146 poisson : RationalNumber -> (() -> Integer) 147 ++ poisson(f) \undocumented 148 geometric : RationalNumber -> (() -> Integer) 149 ++ geometric(f) \undocumented 150 151 ridHack1 : (Integer, Integer, Integer, Integer) -> Integer 152 ++ ridHack1(i, j, k, l) \undocumented 153 == add 154 import from RandomNumberSource() 155 156 -- Compute uniform(a..b) as 157 -- 158 -- l + U0 + w*U1 + w^2*U2 +...+ w^(n-1)*U-1 + w^n*M 159 -- 160 -- where 161 -- l = min(a, b) 162 -- m = abs(b-a) + 1 163 -- w^n < m < w^(n+1) 164 -- U0, ..., Un-1 are uniform on 0..w-1 165 -- M is uniform on 0..(m quo w^n)-1 166 167 uniform aTob == 168 a := low(aTob); b := high(aTob) 169 l := min(a, b); m := abs(a-b) + 1 170 171 w := 2^(length size() quo 2)::NonNegativeInteger 172 173 n := 0 174 mq := m -- m quo w^n 175 while (mqnext := mq quo w) > 0 repeat 176 n := n + 1 177 mq := mqnext 178 ridHack1(mq, n, w, l) 179 180 ridHack1(mq, n, w, l) == 181 r := randnum mq 182 for i in 1..n repeat r := r*w + randnum w 183 r + l 184 185)abbrev package RFDIST RandomFloatDistributions 186++ Description: 187++ This package exports random floating-point distributions 188-- Currently unused. 189RandomFloatDistributions() : Cat == Body where 190 NNI ==> NonNegativeInteger 191 192 Cat ==> with 193 uniform01 : () -> Float 194 ++ uniform01() \undocumented 195 normal01 : () -> Float 196 ++ normal01() \undocumented 197 exponential1 : () -> Float 198 ++ exponential1() \undocumented 199 chiSquare1 : NNI -> Float 200 ++ chiSquare1(n) \undocumented 201 202 uniform : (Float, Float) -> (() -> Float) 203 ++ uniform(f, g) \undocumented 204 normal : (Float, Float) -> (() -> Float) 205 ++ normal(f, g) \undocumented 206 exponential : (Float) -> (() -> Float) 207 ++ exponential(f) \undocumented 208 chiSquare : (NNI) -> (() -> Float) 209 ++ chiSquare(n) \undocumented 210 Beta : (NNI, NNI) -> (() -> Float) 211 ++ Beta(n, m) \undocumented 212 F : (NNI, NNI) -> (() -> Float) 213 ++ F(n, m) \undocumented 214 t : (NNI) -> (() -> Float) 215 ++ t(n) \undocumented 216 217 218 Body ==> add 219 import from RandomNumberSource() 220-- FloatPackage0() 221 222 -- random() generates numbers in 0..rnmax 223 rnmax := (size()$RandomNumberSource() - 1)::Float 224 225 uniform01() == 226 randnum()::Float/rnmax 227 uniform(a, b) == 228 a + uniform01()*(b-a) 229 230 exponential1() == 231 u : Float := 0 232 -- This test should really be u < m where m is 233 -- the minumum acceptible argument to log. 234 while u = 0 repeat u := uniform01() 235 - log u 236 exponential(mean) == 237 mean*exponential1() 238 239 -- This method is correct but slow. 240 normal01() == 241 s := 2::Float 242 while s >= 1 repeat 243 v1 := 2 * uniform01() - 1 244 v2 := 2 * uniform01() - 1 245 s := v1^2 + v2^2 246 v1 * sqrt(-2 * log s/s) 247 normal(mean, stdev) == 248 mean + stdev*normal01() 249 250 chiSquare1 dgfree == 251 x : Float := 0 252 for i in 1..dgfree quo 2 repeat 253 x := x + 2*exponential1() 254 if odd? dgfree then 255 x := x + normal01()^2 256 x 257 chiSquare dgfree == 258 chiSquare1 dgfree 259 260 Beta(dgfree1, dgfree2) == 261 y1 := chiSquare1 dgfree1 262 y2 := chiSquare1 dgfree2 263 y1/(y1 + y2) 264 265 F(dgfree1, dgfree2) == 266 y1 := chiSquare1 dgfree1 267 y2 := chiSquare1 dgfree2 268 (dgfree2 * y1)/(dgfree1 * y2) 269 270 t dgfree == 271 n := normal01() 272 d := chiSquare1(dgfree) / (dgfree::Float) 273 n / sqrt d 274 275--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. 276--All rights reserved. 277-- 278--Redistribution and use in source and binary forms, with or without 279--modification, are permitted provided that the following conditions are 280--met: 281-- 282-- - Redistributions of source code must retain the above copyright 283-- notice, this list of conditions and the following disclaimer. 284-- 285-- - Redistributions in binary form must reproduce the above copyright 286-- notice, this list of conditions and the following disclaimer in 287-- the documentation and/or other materials provided with the 288-- distribution. 289-- 290-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the 291-- names of its contributors may be used to endorse or promote products 292-- derived from this software without specific prior written permission. 293-- 294--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 295--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 296--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 297--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 298--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 299--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 300--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 301--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 302--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 303--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 304--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 305