\documentclass{article}
\usepackage{open-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 negative? t 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
        import %ilength: Integer -> Integer        from Foreign Builtin
        import %ibit: (Integer,Integer) -> Boolean from Foreign Builtin

        bitLength n   == %ilength n
        bitCoef (n,i) == if %ibit(n,i) then 1 else 0
        bitTruth(n,i) == %ibit(n,i)

@
\section{package RIDIST RandomIntegerDistributions}
<<package RIDIST RandomIntegerDistributions>>=
)abbrev package RIDIST RandomIntegerDistributions
++ Description:
++ This package exports integer distributions
RationalNumber==> Fraction Integer
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 positive?(mqnext := mq quo w) 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
            v1: 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}