Re: [Axiom-mail] Random numbers in axiom

[root]

Chris,

The function you want is called "reseed". 
You can get the current value using "seed".

(1) -> seed()

    (1) 694953969372946172
                                                 Type: PositiveInteger

(2) -> reseed(6)
                                                 Type: Void

(3) seed()

    (3) 6
                                                 Type: PositiveInteger

I've attached the code for the domain "random.spad" which contains
some of the random functions in Axiom. Axiom can also generate many
other type of random objects such as permuations and polynomials. 
To see the functions named "random" in Axiom type:

(4) -> )d op random

There are 8 exposed functions called random :
   [1]  -> D from D if D has FINITE
   [2] D -> D from D if D has INS
   [3]  -> D from D if D has INS
   [4] Integer -> Integer from Integer
   [5] NonNegativeInteger -> NonNegativeInteger from NonNegativeInteger
            
   [6] PermutationGroup D2 -> Permutation D2 from PermutationGroup D2
            if D2 has SETCAT
   [7] (PermutationGroup D3,Integer) -> Permutation D3
            from PermutationGroup D3 if D3 has SETCAT
   [8]  -> D from D if D has QFCAT D1 and D1 has INS and D1 has INTDOM
            

There are 3 unexposed functions called random :
   [1] PositiveInteger -> SparseUnivariatePolynomial D3
            from FiniteFieldPolynomialPackage D3 if D3 has FFIELDC
   [2] (PositiveInteger,PositiveInteger) -> SparseUnivariatePolynomial 
            D3
            from FiniteFieldPolynomialPackage D3 if D3 has FFIELDC
   [3] PositiveInteger -> Vector D3 from InnerNormalBasisFieldFunctions
            D3
            if D3 has FFIELDC


Tim Daly
address@hidden
address@hidden


=========================================================================
=== random.spad
=========================================================================


)abbrev package RANDSRC RandomNumberSource
++ Description:Random number generators
++  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

)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)

)abbrev package INTBIT IntegerBits
++ 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

)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

)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()

        -- 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