rnd.sa
Generated by gen_html_sa_files from ICSI. Contact gomes@icsi.berkeley.edu for details
---------------------------> Sather 1.1 source file <--------------------------
-- Copyright (C) International Computer Science Institute, 1994. COPYRIGHT --
-- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
-- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in --
-- the file "Doc/License" of the Sather distribution. The license is also --
-- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA. --
--------> Please email comments to sather-bugs@icsi.berkeley.edu. <----------
-- 1995/Aug7 Blind conversion of 0.2 library by gomes.
-- to_i changes to .floor.int (old behavior of .to_i was = floor)
-- Involves DOUBLE->FLTD, adding of .0d after literals,
-- change from log(v) exp(v) etc to v.log v.exp etc.
-- New tests are all *UNCHECKED!* (to be fair, so were the tests in 0.2!)
-- changes 1994/11/03 by Erik Schnetter
-- rnd.sa: Random numbers.
-- Originally from the 0.2 version by S. Omohundro. Converted to 1.0 by
-- gomes The standard version RND::uniform and RND::int have been tested.
abstract class $RANDOM_GEN
abstract class $RANDOM_GEN is
-- Supertype from which all "raw" random number generators subtype.
-- Currently only supports the minimimal standard generator.
-- RND uses these subtypes to build more complex random numbers
create: $RANDOM_GEN;
init(seed: INT);
-- Initialize the generator using nseed. Any `INT' value should be legal.
get:FLTD;
-- The next random value. Should be in `[0.,1.)'.
end;
class RND{GEN < $RANDOM_GEN}
class RND{GEN < $RANDOM_GEN} is
-- Random number source, which uses a "raw" random number generator
-- stream of type "GEN"
-- If this is used through a regular variable, which has been created,
-- it will act as an individual generator.
-- Otherwise - if the variable is void, or routines are called on
-- the class directly - a default generator will be used
-- act as an independant random number generator.
-- Tries to satisfy a number of needs. In the most common situation,
-- a single generator will supply all random numbers for the system.
-- Calls of the form: `RND::uniform' will use this shared generator
-- (as will `r.uniform' where `r:RANDOM' is `void'). To use multiple
-- generators (perhaps of different kinds) one creates objects
-- of type `RND'. These objects have an attribute `own_gen' of type
-- `$RANDOM_GEN' which holds the actual generator. Several types of
-- generator are provided and new ones may be added by having them
-- subtype from the specification class `$RANDOM_GEN'.
-- Classes are provided to combine
-- generators and to randomly permute the output of a generator.
private shared default_gen: GEN; -- Shared generator
private attr own_gen: GEN; -- Individualized generator
gen: GEN is
-- Interface to the generator. If self is void, use the shared
-- generator, otherwise use the created generator
if void(self) then
if void(default_gen) then default_gen := make_generator; end;
return default_gen;
else return own_gen end;
end;
create: SAME is
-- Creates a new stream of random values.
res::=new;
res.own_gen := make_generator;
return(res)
end;
private make_generator: GEN is
-- Since $RANDOM_GEN::create returns a $RANDOM_GEN, we have to
-- typecase to make sure that it actually is of type GEN (abstract
-- classes cannot return SAME)
g ::= #GEN;
typecase g
when GEN then return g
else raise("Generator's create does not return SAME"); end;
end;
seed(seed: INT) is
-- Initialize the generator with "seed"
-- Can use with assignment notation
gen.init(seed);
end; -- init
int(l,u:INT):INT pre l <= u is
-- A uniform random member of `{l,...,u}'.
return(l+(((u-l+1).fltd)*uniform).floor.int)
end;
uniform: FLTD is
-- A uniformly distributed double in `[0.,1.)'.
return gen.get
end; -- uniform
uniform_range(l,u: FLTD): FLTD pre l <= u is
-- A uniformly distributed double in the range `[l,u)'.
return l + (u-l)*uniform
end;
bit(p:FLTD):BOOL pre p >= 0.0d and p <= 1.0d is
-- True with probability `p'.
return (uniform<p)
end; -- bit
-- Standard normal generator
private attr norm_cache:FLTD; -- Cache of normal deviate
-- generated earlier.
private shared snorm_cache:FLTD; -- Shared version of above.
private attr norm_cache_valid:BOOL; -- `True' if a value is cached.
private shared snorm_cache_valid:BOOL; -- Shared version of above.
standard_normal:FLTD is
-- A normally distributed `FLTD' with mean `0.' and unit variance.
-- Uses Box-Muller method which generates two normal deviates and
-- requires two calls to uniform. One of these is cached for the next
-- call. Based on Knuth, Vol. 2, p. 117.
res: FLTD;
if void(self) then -- use shareds
if snorm_cache_valid then -- use stored value
snorm_cache_valid:=false;
res:=snorm_cache;
else -- must recompute
s:FLTD:=10.0d; v1,v2:FLTD;
loop
-- get two samples inside unit circle
until!(s<1.0d and s/=0.0d);
v1:=uniform_range(-1.0d,1.0d);
v2:=uniform_range(-1.0d,1.0d);
s:=v1*v1+v2*v2;
end; -- loop, executed 1.27 times on average
rt:FLTD:= (-2.0d*(s).log/s).sqrt;
snorm_cache:=v1*rt;
snorm_cache_valid:=true;
res:=v2*rt;
end; -- if
else -- use locals
if norm_cache_valid then -- use stored value
res:=norm_cache; norm_cache_valid:=false;
else -- must recompute
s:FLTD:=10.0d; v1,v2:FLTD;
loop -- get two samples inside unit circle
until!(s<1.0d and s/=0.0d);
v1:=uniform_range(-1.0d,1.0d); v2:=uniform_range(-1.0d,1.0d);
s:=v1*v1+v2*v2;
end; -- loop, executed 1.27 times on average
rt:FLTD:= (-2.0d*(s).log/s).sqrt;
norm_cache:=v1*rt; norm_cache_valid:=true;
res:=v2*rt;
end; -- if
end; -- if
return res;
end; -- standard_normal
exponential(m:FLTD):FLTD pre m>=0.0d is
-- Samples from an exponential distribution: `p(x)=e^(-x/m)/m'.
-- Mean is `m'. Eg: Time between emissions when one per `m' seconds on
-- average. Based on Knuth, Vol. 2, p. 128.
u:FLTD; loop until!(u/=0.0d); u:=uniform end;
return -m*u.log;
end; -- exponential
gamma(a:FLTD):FLTD pre a > 0.0d is
-- Samples from the Gamma distribution: `p(x)=(x^(a-1)e^(-x))/G(a)'.
-- Mean is `a>0'. Uses Cheng's algorithm from Devroye, "Non-Uniform
-- Random Variate Generation", Springer-Verlag, (1986) p. 413.
-- This is much less restrictive than the Numerical Recipes algorithm.
res: FLTD;
b:FLTD:=a-(4.0d).log;
c:FLTD:=a+(2.0d*a-1.0d).sqrt;
loop
u:FLTD:=uniform; v:FLTD:=uniform;
if v=0.0d or v=1.0d then v:=0.5d end;
y:FLTD:=a*(v/(1.0d-v)).log;
res:=a*(v).exp;
z:FLTD:=u*v*v;
r:FLTD:=b+c*y-res;
if r>=(9.0d/2.0d)*z-(1.0d+(9.0d/2.0d).log) then break! end;
if r>=(z).log then break! end;
end; -- loop
return res;
end; -- gamma
beta(a,b:INT):FLTD pre a >0 and b > 0 is
-- Samples from the Beta distribution
-- `p(x)=(G(a+b)/(G(a)G(b)))x^(a-1)(1-x)^(b-1)', `0<x<1', `a,b>0'.
-- Mean is `a/(a+b)'. Knuth, Vol. 2, p. 129.
x1:FLTD:=gamma(a.fltd);
if x1=0.0d then return 0.0d;
else return x1/(x1+gamma(b.fltd)); end;
end; -- beta
chi_square(v:INT):FLTD is
-- Samples from the Chi-Square distribution:
-- `p(x)=e^(-x/2)x^((v/2)-1)/(2^(v/2)G(v/2))'. Mean is `v'. `x>=0'.
-- Based on Knuth, Vol. 2, p. 130.
return 2.0d*gamma(v.fltd/2.0d);
end; -- chi_square
f_dist(v1,v2:INT):FLTD is
-- Samples from the F-distribution with `v1' and `v2' degrees of
-- freedom (see any statistics text). Based on Knuth, Vol. 2, p. 130.
return chi_square(v1)*v2.fltd/(chi_square(v2)*v1.fltd);
end; -- f_dist
t_dist(v:INT):FLTD is
-- Samples from the t-distribution with `v' degrees of freedom (see
-- any statistics text). Based on Knuth, Vol. 2, p. 130.
return standard_normal/((chi_square(v)/v.fltd).sqrt);
end; -- t_dist
geometric(p:FLTD):INT pre p >= 0.0d and p <= 1.0d is
-- Samples from the geometric distribution. Is the number of trials
-- needed until an event of probability `p' first occurs.
-- `P(n)=(1-p)^(n-1)p'. Mean is `1/p'. Based on Knuth, Vol. 2, p. 131.
if p=1.0d then return 1
else return (uniform.log/(1.0d-p).log).ceiling.int; end;
end; -- geometric
binomial(n:INT, p:FLTD):INT pre p >= 0.0d and p <= 1.0d is
-- Samples from the binomial distribution:
-- `P(t)=(n t) p^t (1-p)^(n-t)'. Mean is `np'. Is number of
-- occurrences of a probability `p' event in `n' trials.
-- Based on Knuth, Vol. 2, p. 131. and Numerical Recipes in C, p. 223.
res: INT;
dr:FLTD; -- FLTD version of res.
-- can flip p->1-p if we also do res->n-res
lp:FLTD; if p<0.5d then lp:=p else lp:=1.0d-p end;
am:FLTD:=n.fltd*lp; -- mean of res
if n<25 then -- Use direct method if n small enough.
i:INT; loop until!(i=n);
if uniform<lp then dr:=dr+1.0d end;
i:=i+1;
end; -- loop
elsif am<1.0d then -- such small prob, approximate by Poisson
dr:=poisson(am).fltd; -- poisson distribution with the same mean
if dr>n.fltd then dr:=n.fltd end; -- chop the tail
else -- use rejection method
y:FLTD; t:FLTD;
sq:FLTD:=(2.0d*am*(1.0d-lp));
loop
loop -- rejection with Lorenzian comparison
y:=cauchy; dr:=sq*y+am;
if dr>=0.0d and dr<n.fltd+1.0d then break! end;
end; -- loop
dr:=dr.floor;
if uniform<=1.2d*sq*(1.0d+y*y)
*((n.fltd+1.0d).log_gamma
-(dr+1.0d).log_gamma
-(n.fltd-dr+1.0d).log_gamma
+dr*(lp).log+(n.fltd-dr)*(1.0d-lp).log)
.exp
then break! end;
end; -- loop
end; -- if
if lp/=p then res:=n-dr.floor.int else res:=dr.floor.int end;
-- undo symmetry
return res;
end; -- binomial
poisson(x:FLTD):INT is
-- Samples from the Poisson distribution: `p(j)=x^j e^(-x)/j!'.
-- Mean is `x'. Is number of events in time interval `x'.
-- Based on Numerical Recipes in C, p. 221.
res: INT;
if (x<12.0d) then -- use direct method
res:=-1; t:FLTD:=1.0d; g:FLTD:=(-x).exp;
loop
res:=res+1; t:=t*uniform; if t<=g then return res end;
end; -- loop
else -- use rejection method
rd:FLTD; -- FLTD version of res
loop
y:FLTD; g:FLTD:=x*(x.log)-(x+1.0d).log_gamma;
sq:FLTD:=(2.0d*x).sqrt;
loop -- use Lorentzian comparison
y:=cauchy; rd:=sq*y+x; if rd>=0.0d then break! end;
end; -- loop
rd:=rd.floor;
if uniform<=0.9d*(1.0d+y*y)*(rd*x.log-(rd+1.0d).log_gamma-g).exp
then break! end;
end; -- loop
res:=rd.floor.int;
end; -- if
return res;
end; -- poisson
cauchy:FLTD is
-- Samples from the Cauchy distribution: `p(x)=1/pi(1+x^2)'.
-- (Also called the Lorentzian.) No mean. Based on NR in C. p. 220.
return (FLTD::pi*uniform).tan
end; -- cauchy
cantor(s:FLTD):FLTD is
-- A uniform random sample from a Cantor set with scaling `s'.
res: FLTD;
r:FLTD:=uniform; m:FLTD:=1.0d-s;
i:INT:=0; loop until!(i>33);
if r>=0.5d then r:=r-0.5d; res:=res+m end;
r:=r*2.0d; m := m*s; i:=i+1
end; -- loop
return res;
end; -- cantor
standard_cantor:FLTD is
-- A uniform random sample from the canonical Cantor set in `[0,1)'.
return cantor(0.3333333333d)
end; -- standard_cantor
end; -- class RND
class RND
class RND is
include RND{MS_RANDOM_GEN};
end;
class MS_RANDOM_GEN < $RANDOM_GEN
class MS_RANDOM_GEN < $RANDOM_GEN is
-- The "minimal standard" generator described in "Random Number
-- Generators: Good Ones are Hard to Find" by Stephen Park and
-- Keith Miller, Communications of the ACM, October 1988, Volume
-- 31, Number 10, p. 1192. Linear congruential, produces a value
-- in `[0.,1.)' including `0.' but not `1.' Any seed value in
-- the range `[1,2147483646]' is equally good.
-- Constants used in generator:
-- BEWARE!! Problems with order of initialization and
-- and FLTD literals
-- const ms_a:FLTD:=16807.0d; -- `7^5'
-- const ms_m:FLTD:=2147483647.0d; -- `(2^31)-1', prime
-- const ms_md:FLTD:=(ms_m-1.0d); -- to avoid continual recomputation
const ms_a:FLTD:=16807.0d; -- `7^5'
shared ms_m:FLTD; -- `(2^31)-1', prime
shared ms_md:FLTD; -- to avoid continual recomputation
attr seed:INT; -- Current state of generator.
create:SAME is
-- A minimal standard generator with `seed=1'.
res ::= new; res.init(1); return(res)
end;
init(nseed:INT) is
-- Initialize the generator.
seed:=1+(nseed-1).mod(2147483645);
ms_m := (1.lshift(31)).fltd-1.0d;
ms_md := ms_m - 1.0d;
end; -- keep in legal range
get:FLTD is
-- Pseudo-random value in `[0.,1.)' generated by minimal
-- std generator.
tmp:FLTD:=ms_a*(seed.fltd);
-- Original: tmp:FLTD:=ms_a*(seed);
seed:=(tmp-(ms_m*((tmp/ms_m).floor))).floor.int;
-- Original: seed:=(tmp-(ms_m*((tmp/ms_m).to_i))).to_i
return(((seed-1).fltd)/ms_md);
-- Original: return(((seed-1))/ms_md);
end; -- get
end; -- class MS_RANDOM_GEN