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