Initial population.

1 files changed, 348 insertions, 0 deletions

diff --git a/src/algebra/random.spad.pamphlet b/src/algebra/random.spad.pamphlet
new file mode 100644
index 00000000..9699789c
--- /dev/null
+++ b/src/algebra/random.spad.pamphlet
@@ -0,0 +1,348 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra random.spad}
+\author{Stephen M. Watt, Mike Dewar}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package RANDSRC RandomNumberSource}
+<<package RANDSRC RandomNumberSource>>=
+)abbrev package RANDSRC RandomNumberSource
+++ Author:S.M.Watt
+++ Date Created: April 87
+++ Date Last Updated:Jan 92, May 1995 (MCD)
+++ Basic Operations:
+++ Related Domains:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description:Random number generators
+--% RandomNumberSource
+++  All random numbers used in the system should originate from
+++  the same generator.  This package is intended to be the source.
+--
+--  Possible improvements:
+--  1) Start where the user left off
+--  2) Be able to switch between methods in the random number source.
+RandomNumberSource(): with
+    -- If r := randnum() then  0 <= r < size().
+        randnum:  () -> Integer
+           ++ randnum() is a random number between 0 and size().
+    --   If r := randnum() then  0 <= r < size().
+        size:     () -> Integer
+          ++ size() is the base of the random number generator
+ 
+        -- If r := randnum n and n <= size()  then  0 <= r < n.
+        randnum:  Integer -> Integer
+           ++ randnum(n) is a random number between 0 and n.
+        reseed:   Integer -> Void
+           ++ reseed(n) restarts the random number generator at n.
+        seed : () -> Integer
+           ++ seed() returns the current seed value.
+ 
+    == add
+        -- This random number generator passes the spectral test
+        -- with flying colours. [Knuth vol2, 2nd ed, p105]
+        ranbase: Integer := 2**31-1
+        x0:   Integer := 1231231231
+        x1:   Integer := 3243232987
+ 
+        randnum() ==
+            t := (271828183 * x1 - 314159269 * x0) rem ranbase
+            if t < 0 then t := t + ranbase
+            x0:= x1
+            x1:= t
+ 
+        size() == ranbase
+        reseed n ==
+            x0 := n rem ranbase
+            -- x1 := (n quo ranbase) rem ranbase
+            x1 := n quo ranbase
+
+        seed() == x1*ranbase + x0
+ 
+        -- Compute an integer in 0..n-1.
+        randnum n ==
+            (n * randnum()) quo ranbase
+
+@
+\section{package RDIST RandomDistributions}
+<<package RDIST RandomDistributions>>=
+)abbrev package RDIST RandomDistributions
+++ Description:
+++ This package exports random distributions
+RandomDistributions(S: SetCategory): with
+        uniform:  Set S -> (() -> S)
+		++ uniform(s) \undocumented
+        weighted: List Record(value: S, weight: Integer) -> (()->S)
+		++ weighted(l) \undocumented
+        rdHack1:  (Vector S,Vector Integer,Integer)->(()->S)
+		++ rdHack1(v,u,n) \undocumented
+    == add
+        import RandomNumberSource()
+
+        weighted lvw ==
+            -- Collapse duplicates, adding weights.
+            t: Table(S, Integer) := table()
+            for r in lvw repeat
+                u := search(r.value,t)
+                w := (u case "failed" => 0; u::Integer)
+                t r.value := w + r.weight
+
+            -- Construct vectors of values and cumulative weights.
+            kl := keys t
+            n  := (#kl)::NonNegativeInteger
+            n = 0 => error "Cannot select from empty set"
+            kv: Vector(S)       := new(n, kl.0)
+            wv: Vector(Integer) := new(n, 0)
+
+            totwt: Integer := 0
+            for k in kl for i in 1..n repeat
+                kv.i := k
+                totwt:= totwt + t k
+                wv.i := totwt
+
+            -- Function to generate an integer and lookup.
+            rdHack1(kv, wv, totwt)
+
+        rdHack1(kv, wv, totwt) ==
+            w := randnum totwt
+            -- do binary search in wv
+            kv.1
+
+        uniform fset ==
+            l := members fset
+            n := #l
+            l.(randnum(n)+1)
+
+@
+\section{package INTBIT IntegerBits}
+<<package INTBIT IntegerBits>>=
+)abbrev package INTBIT IntegerBits
+----> Bug! Cannot precompute params and return a function which
+----> simpy computes the last call.  e.g. ridHack1, below.
+
+--% IntegerBits
+--  Functions related to the binary representation of integers.
+--  These functions directly access the bits in the big integer
+--  representation and so are much facter than using a quotient loop.
+--  SMW Sept 86.
+
+ 
+++ Description:
+++ This  package provides functions to lookup bits in integers 
+IntegerBits: with
+    --  bitLength(n)  == # of bits to represent abs(n)
+    --  bitCoef (n,i) == coef of 2**i in abs(n)
+    --  bitTruth(n,i) == true if coef of 2**i in abs(n) is 1
+ 
+        bitLength: Integer -> Integer
+		++ bitLength(n) returns the number of bits to represent abs(n)
+        bitCoef:   (Integer, Integer) -> Integer
+		++ bitCoef(n,m) returns the coefficient of 2**m  in abs(n)
+        bitTruth:  (Integer, Integer) -> Boolean
+		++ bitTruth(n,m) returns true if coefficient of 2**m in abs(n) is 1
+ 
+    == add
+        bitLength n   == INTEGER_-LENGTH(n)$Lisp
+        bitCoef (n,i) == if INTEGER_-BIT(n,i)$Lisp then 1 else 0
+        bitTruth(n,i) == INTEGER_-BIT(n,i)$Lisp
+
+@
+\section{package RIDIST RandomIntegerDistributions}
+<<package RIDIST RandomIntegerDistributions>>=
+)abbrev package RIDIST RandomIntegerDistributions
+++ Description:
+++ This package exports integer distributions
+RandomIntegerDistributions(): with
+        uniform:   Segment Integer           -> (() -> Integer)
+		++ uniform(s) \undocumented
+        binomial:  (Integer, RationalNumber) -> (() -> Integer)
+		++ binomial(n,f) \undocumented
+        poisson:   RationalNumber          -> (() -> Integer)
+		++ poisson(f) \undocumented
+        geometric: RationalNumber          -> (() -> Integer)
+		++ geometric(f) \undocumented
+
+        ridHack1:  (Integer,Integer,Integer,Integer) -> Integer
+		++ ridHack1(i,j,k,l) \undocumented
+    == add
+        import RandomNumberSource()
+        import IntegerBits()
+
+        -- Compute uniform(a..b) as
+        --
+        --     l + U0 + w*U1 + w**2*U2 +...+ w**(n-1)*U-1 + w**n*M
+        --
+        -- where
+        --     l = min(a,b)
+        --     m = abs(b-a) + 1
+        --     w**n < m < w**(n+1)
+        --     U0,...,Un-1  are uniform on  0..w-1
+        --     M            is  uniform on  0..(m quo w**n)-1
+
+        uniform aTob ==
+            a := lo aTob;  b := hi aTob
+            l := min(a,b); m := abs(a-b) + 1
+
+            w := 2**(bitLength size() quo 2)::NonNegativeInteger
+
+            n  := 0
+            mq := m  -- m quo w**n
+            while (mqnext := mq quo w) > 0 repeat
+                n  := n + 1
+                mq := mqnext
+            ridHack1(mq, n, w, l)
+
+        ridHack1(mq, n, w, l) ==
+            r := randnum mq
+            for i in 1..n repeat r := r*w + randnum w
+            r + l
+
+@
+\section{package RFDIST RandomFloatDistributions}
+<<package RFDIST RandomFloatDistributions>>=
+)abbrev package RFDIST RandomFloatDistributions
+++ Description:
+++ This package exports random floating-point distributions
+RationalNumber==> Fraction Integer
+RandomFloatDistributions(): Cat == Body where
+    NNI ==> NonNegativeInteger
+
+    Cat ==> with
+        uniform01:   ()  -> Float
+		++ uniform01() \undocumented
+        normal01:    ()  -> Float
+		++ normal01() \undocumented
+        exponential1:()  -> Float
+		++ exponential1() \undocumented
+        chiSquare1:  NNI -> Float
+		++ chiSquare1(n) \undocumented
+
+        uniform:     (Float, Float) -> (() -> Float)
+		++ uniform(f,g) \undocumented
+        normal:      (Float, Float) -> (() -> Float)
+		++ normal(f,g) \undocumented
+        exponential: (Float)        -> (() -> Float)
+		++ exponential(f) \undocumented
+        chiSquare:   (NNI)          -> (() -> Float)
+		++ chiSquare(n) \undocumented
+        Beta:        (NNI, NNI)     -> (() -> Float)
+		++ Beta(n,m) \undocumented
+        F:           (NNI, NNI)     -> (() -> Float)
+		++ F(n,m) \undocumented
+        t:           (NNI)          -> (() -> Float)
+		++ t(n) \undocumented
+
+
+    Body ==> add
+        import RandomNumberSource()
+--      FloatPackage0()
+
+        -- random()  generates numbers in 0..rnmax
+        rnmax := (size()$RandomNumberSource() - 1)::Float
+
+        uniform01() ==
+            randnum()::Float/rnmax
+        uniform(a,b) ==
+            a + uniform01()*(b-a)
+
+        exponential1() ==
+            u: Float := 0
+            -- This test should really be  u < m where m is
+            -- the minumum acceptible argument to log.
+            while u = 0 repeat u := uniform01()
+            - log u
+        exponential(mean) ==
+            mean*exponential1()
+
+        -- This method is correct but slow.
+        normal01() ==
+            s := 2::Float
+            while s >= 1 repeat
+                v1 := 2 * uniform01() - 1
+                v2 := 2 * uniform01() - 1
+                s  := v1**2 + v2**2
+            v1 * sqrt(-2 * log s/s)
+        normal(mean, stdev) ==
+            mean + stdev*normal01()
+
+        chiSquare1 dgfree ==
+            x: Float := 0
+            for i in 1..dgfree quo 2 repeat
+                x := x + 2*exponential1()
+            if odd? dgfree then
+                x := x + normal01()**2
+            x
+        chiSquare dgfree ==
+            chiSquare1 dgfree
+
+        Beta(dgfree1, dgfree2) ==
+            y1 := chiSquare1 dgfree1
+            y2 := chiSquare1 dgfree2
+            y1/(y1 + y2)
+
+        F(dgfree1, dgfree2) ==
+            y1 := chiSquare1 dgfree1
+            y2 := chiSquare1 dgfree2
+            (dgfree2 * y1)/(dgfree1 * y2)
+
+        t dgfree ==
+            n := normal01()
+            d := chiSquare1(dgfree) / (dgfree::Float)
+            n / sqrt d
+
+@
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+ 
+<<package RANDSRC RandomNumberSource>>
+<<package RDIST RandomDistributions>>
+<<package INTBIT IntegerBits>>
+<<package RIDIST RandomIntegerDistributions>>
+<<package RFDIST RandomFloatDistributions>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}