[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Axiom-mail] Random numbers in axiom
From: |
root |
Subject: |
Re: [Axiom-mail] Random numbers in axiom |
Date: |
Tue, 27 Apr 2004 07:40:24 -0400 |
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