[Back to MATH SWAG index] [Back to Main SWAG index] [Original]
{$A-,B-,D+,E+,F-,I+,L+,N+,O-,R+,S+,V+}
UNIT Rangen;
{==============================================================================
MODULE: RANDOMGENERATOR (RANGEN.PAS)
===============================================================================
VERSION: V1.5 Turbo-Pascal 6.0, 7.0
COPYRIGHT: (c) 1987-1993
Dr. W. Gross, Keplerstrasse 51
D-69120 Heidelberg, FRG
gross@aecds.exchi.uni-heidelberg.de
CREATED: 10-NOV-87 Wolfgang Gross on VAX/VMS, PASCAL-2
LAST UPDATE: 02-MAY-88 Wolfgang Gross round off error in
NORRANDOM
21-JUL-88 Wolfgang Gross round off error in
UNIRANDOM
15-MAR-89 Wolfgang Gross NorRan, RanSphere
TESTED: 10-FEB-90 Wolfgang Gross
-------------------------------------------------------------------------------
TITLE: RANDOM NUMBER GENERATES
-------------------------------------------------------------------------------
DESCRIPTION: Random number generators for uniform, normal, exponential,
lognormal, poisson and uniform sphere distribution.
If the TURBO-PASCAL random generator RANDOM is used
one must set the variable SYSTEM.RandSeed before calling
RANDOM and save it after the call using the given
seed variable. This allows independent streams of
random numbers. See comment in UNIRANDOM.
If TEST8086>2 we can use 386 instructions. UNIRANDOM
will then use its own method to generate the next
univariate pseudo random number as described below.
The method is approved by many independent tests
that have been published in the literature. It has
good Bayer coefficents which means it can be used for
generating multi-dimensional random vectors.
One might also utilize a floating point processor
or switch to IEEE reals (type single, double) etc.
Everybody may elaborate on this him/herself.
If someone needs to modify and recompile this unit
and doesn't have TASM, the complete ASM code is
provided for compilation by the internal assembler.
Some special arrangements are necessary to cope
with 386 instructions. Define HAVETASM if you want
to link MD2P31.OBJ. If this symbol is not defined
the inline assembler version is used.
This unit has evolved from a similar module worked
out under PASCAL-2 in a VAX/VMS environment (including
VAX assembler for MD2P31). The module can be found
with anon ftp on aecds.exchi.uni-heidelberg.de.
CALLED BY: <module name> <title>
EXPORTS: UNIRANDOM uniformly distributed random numbers
NORRANDOM normally distributed random numbers
NORRANDOM2 normally distributed random numbers,
polar method
NORRANDOM3 normally distributed random numbers,
other mean/std-dev values than (0,1)
EXPRANDOM exponentially distributed random numbers
LOGNORRANDOM logarithmic normal distribution
POIRANDOM Poisson distribution
RanSphere uniformly distri. random numbers on sphere
INPUT: seed random number seed
mean mean value for distribution
std standard deviation
OUTPUT: function value = random number
}
{=============================================================================}
{ Exports: }
Interface
FUNCTION UNIRANDOM ( VAR seed : longint ) : real;
FUNCTION NORRANDOM ( VAR seed : longint ) : real;
FUNCTION NORRANDOM2 ( VAR seed : longint;
VAR v2 : real ) : real;
FUNCTION NORRANDOM3 ( mean, std: real;
VAR seed : longint ) : real;
FUNCTION EXPRANDOM ( mean : real;
VAR seed : longint ) : real;
FUNCTION LOGNORRANDOM ( mean, std : real;
VAR seed : longint ) : real;
FUNCTION POIRANDOM ( lambda : real;
VAR seed : longint ) : longint;
PROCEDURE RanSphere ( VAR seed : longint;
VAR x,y,z: real );
{===========================================================================}
Implementation
{DEFINE HAVETASM} {put a $ in front of DEFINE to link MD2P31.OBJ}
{$IFDEF HAVETASM}
PROCEDURE MD2P31 (A,B : longint;
VAR Q,R : longint); far; external;
{$L MD2P31.OBJ}
{$ELSE}
PROCEDURE MD2P31 (A,B : longint;
VAR Q,R : longint); far; assembler;
{rather tricky: we fake TP into generating 386 intructions
Calculate product a*b, represent product as q*2^31+r, 0<=r<2^31}
asm
push es
push di
db $66 { 386 prefix for dword operation}
push ax { = push eax }
db $66
push bx { = push ebx }
db $66
push dx { = push edx }
db $66
mov ax,word ptr ss:[bp+18] { mov eax,a }
db $66
mov bx,word ptr ss:[bp+14] { mov ebx,b }
db $66
imul bx { imul ebx ; a*b in edx:eax (= 8 bytes !) }
db $66
mov bx,ax { mov ebx,eax }
db $66,$0f,$a4,$da,$01 { SHLD edx,ebx,1 ; edx contains q }
db $66,$25,$ff,$ff,$ff,$7f { and eax,07fffffffh ; eax contains r }
les di,r
db $66
stosw { stosd ; mov eax to r }
db $66
mov ax,dx { mov eax,edx }
les di,q
db $66
stosw { stosd ; mov eax to q }
db $66
pop dx { = pop edx }
db $66
pop bx { = pop ebx }
db $66
pop ax { = pop eax }
pop es
pop di
end;{ PROC MD2P31 }
{$ENDIF}
FUNCTION UNIRANDOM;
{ Cf. Payne,W.H., Rabung, J.R., Bogyo, T.P. :
"Coding the Lehmer pseudo-random number generator.
Comm. ACM 12, 85-86 (1969) }
CONST MODULUS : longint = 2147483647; { = 2^31-1 }
FACTOR : longint = 397204094; { primitive root, e.g. used by SAS }
VAR Q,R,S : longint;
BEGIN { UNIRANDOM }
IF TEST8086>1 THEN
BEGIN
{ Priniple
seed := ( seed*FACTOR ) MOD MODULUS;
but can not use it due to overflow }
MD2P31 ( seed, factor, Q, R );
S := modulus - R;
IF S > Q
THEN seed := Q + R
ELSE seed := Q - S;
{ Single precision version. For more details on the division cf. to
Fishman, G.S., Moore, L.R. :
A statistical evaluation of multiplicative congruential random
number generators with modulus 2^31-1 .
Jour. Amer. Stat. Ass. 77, 129-136 *1982) }
UNIRANDOM := seed/MODULUS;
END
ELSE
{see comment at the beginning of unit, RandSeed from unit SYSTEM}
BEGIN
RandSeed := seed; Unirandom := random; seed := RandSeed;
END;
END; { FUNC UNIRANDOM }
{--------------------------------------------}
FUNCTION NORRANDOM ( VAR seed : longint ) : real;
CONST TwoPi = 2*Pi;
VAR s, r1, r2, z1, z2 : real;
BEGIN { NORRANDOM }
REPEAT r1 := UniRandom ( seed ) UNTIL r1>0;
r2 := UniRandom ( seed );
s := -2*ln(r1);
IF s<0 THEN s := 0; { round off error could give negative value s }
NORRANDOM := SQRT ( s ) * cos ( TwoPi*r2 );
END; { FUNC NORRANDOM }
{--------------------------------------------}
FUNCTION NorRandom2 ( VAR seed : longint;
VAR v2 : real ) : real;
{ Polar method due to Box, Mueller and Marsaglia,
cf. D.E. Knuth, The art of computer porgramming, vol 2, 117
first call: v2 must be > 1000 }
VAR s, t, v1 : real;
BEGIN { NORRANDOM2 }
IF v2<1000
THEN
BEGIN
NORRANDOM2 := v2; v2 := 2000
END
ELSE
BEGIN
REPEAT
v1 := UniRandom ( seed ); v1 := v1+v1-1;
v2 := UniRandom ( seed ); v2 := v2+v2-1;
s := Sqr (v1) + sqr(v2);
UNTIL (s<1) AND (s>0);
T := Sqrt ( -2*Ln(s)/s );
NORRANDOM2 := V1*T;
V2 := V2*T;
END;
END; { FUNC NORRANDOM2 }
{--------------------------------------------}
FUNCTION NORRANDOM3 ( mean, std : real;
VAR seed : longint ) : real;
BEGIN { NORRANDOM3 }
NORRANDOM3 := std*NORRANDOM (seed) + mean;
END; { FUNC NORRANDOM3 }
{--------------------------------------------}
FUNCTION EXPRANDOM ( mean : real;
VAR seed : longint ) : real;
BEGIN { EXPRANDOM }
EXPRANDOM := ( -ln(UNIRANDOM(seed))*mean );
END; { FUNC EXPRANDOM }
{--------------------------------------------}
FUNCTION LOGNORRANDOM ( mean, std : real;
VAR seed : longint ) : real;
VAR m2, s2, sum, mu, sigma : real;
BEGIN { LOGNORRANDOM }
m2 := SQR(mean); s2 := SQR(std); sum := m2+s2;
mu := 0.5*ln(SQR(m2)/sum);
sigma := SQRT (ln(sum/m2));
LOGNORRANDOM := exp (sigma * NORRANDOM (seed) + mu);
END; { FUNC LOGNORRANDOM }
{--------------------------------------------}
FUNCTION POIRANDOM ( lambda : real;
VAR seed : longint ) : longint;
VAR i : longint;
p,q,r,z : real;
BEGIN { POIRANDOM }
{ cf. H.E. Schaffer, Generator of random numbers satisfying the
Poisson distribution, Comm. ACM 13, 49 (1970)
G.S. Fishman, Sampling from the Poisson distribution on
a computer, Computing 17, 147-156 (1976) }
IF lambda<50
THEN
BEGIN
z := exp ( -lambda ); p := 1; i := -1;
REPEAT
i := i+1;
r := UNIRANDOM ( seed );
p := p*r;
UNTIL p<=z;
POIRANDOM := i;
END
ELSE
BEGIN
i := Round ( NORRANDOM3 ( lambda, SQRT(lambda), seed ) );
IF i<0 THEN i :=0;
POIRANDOM := i;
END;
END; { PROC POIRANDOM }
{--------------------------------------------}
PROCEDURE RanSphere ( VAR seed : longint;
VAR x,y,z : real );
{ R.E. Knop, Random Vectors Uniform in Solid Angle,
CACM 13 (1970), 326 }
VAR s, t, v1, v2 : real;
BEGIN { RanSphere }
REPEAT
v1 := UniRandom ( seed ); v1 := v1+v1-1;
v2 := UniRandom ( seed ); v2 := v2+v2-1;
s := Sqr (v1) + sqr(v2);
UNTIL (s<1) AND (s>0);
T := 2*Sqrt ( 1-S );
x := T*V1; y := T*V2; z := S+S-1;
END; { PROC RanSphere }
end. {UNIT RANGEN}
{--------------------- DEMO ---------------------- }
{ ----------- CUT ---------------- }
{$A+,B-,D+,E+,F-,G-,I+,L+,N-,O-,P-,Q-,R+,S+,T-,V+,X+,Y+}
{$M 16384,0,655360}
program testrand;
uses rangen;
{generates two independent streams of uniform random numbers}
VAR r1, r2 : real;
seed1, seed2 : longint;
i : integer;
BEGIN
write ('enter seed1, seed2 = ' ); readln ( seed1, seed2 );
for i := 1 TO 20 DO
BEGIN
r1 := unirandom(seed1);
r2 := unirandom(seed2);
writeln ( seed1:12, r1:12:8, seed2:12, r2:12:8 );
END;
readln;
END.
{ --------------------------------- ASM UNIT THAT CAN BE USED IN THIS UNIT ------------- }
{ ----------- CUT ---------------- }
; MD2P31.ASM, for Turbo Assembler
.MODEL TPASCAL ; 16-bit segments
.386C ; non-privileged 386 instructions
CODE SEGMENT BYTE PUBLIC
ASSUME cs:CODE,ds:NOTHING
; PASCAL declaration
; PROCEDURE MD2P31 (A,B : longint;
; VAR Q,R : longint); far; external;
; calculate product a*b, represent product as q*2^31+r, 0<=r<2^31
; return q and r
; Parameters (+2 because of push bp)
R EQU DWORD PTR ss:[bp+6]
Q EQU DWORD PTR ss:[bp+10]
B EQU DWORD PTR ss:[bp+14]
A EQU DWORD PTR ss:[bp+18]
MD2P31 PROC FAR
PUBLIC MD2P31
push bp
mov bp,sp ;get pointer into stack
push es
push di ;manual says we don't need to save
push eax ;those registers, but safety first!
push ebx
push edx
mov eax,a
mov ebx,b
imul ebx ; a*b in edx:eax
mov ebx,eax
SHLD edx,ebx,1 ; edx contains q
and eax,07fffffffh ; eax contains r
les di,r
stosd
mov eax,edx
les di,q
stosd
pop edx
pop ebx
pop eax
pop di
pop es
pop bp
retf 16 ;parameters take 16 bytes
MD2P31 ENDP
CODE ENDS
END
[Back to MATH SWAG index] [Back to Main SWAG index] [Original]