/
random.ps
87 lines (76 loc) · 2.82 KB
/
random.ps
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Postscript is Copyright (c) 2017, Peter J Billam %
% Permission is granted to any individual or institution to use, copy, %
% modify or redistribute this software, so long as it is not resold for %
% profit, and provided this notice is retained. It is provided "as is", %
% without any express or implied warranty. http://www.pjb.com.au %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/gauss_rand_r 0.5 1024.0 div 1024.0 div 1024.0 div def % rand scaling factor
usertime realtime add srand
% The algorithm by Richard P. Brent (1973), Computer Centre, ANU, Canberra
% Published (1981) as CACM Algorithm 448
% http://wwwmaths.anu.edu.au/~brent/pd/rpb023i.pdf
% was replaced in 20170510, fixing some obscure bugs ...
% http://www.design.caltech.edu/erik/Misc/Gaussian.html
% The polar form of the Box-Muller transformation is faster and more robust:
% float x1, x2, w, y1, y2;
% do {
% x1 = 2.0 * ranf() - 1.0;
% x2 = 2.0 * ranf() - 1.0;
% w = x1 * x1 + x2 * x2;
% } while ( w >= 1.0 );
% w = sqrt( (-2.0 * log( w ) ) / w );
% y1 = x1 * w;
% y2 = x2 * w;
% where ranf() obtains a random number uniformly distributed in [0,1]
/gauss_rand_already false def
/gauss_rand {
gauss_rand_already {
/gauss_rand_already false def
gauss_rand_y2 % alternate between this and gauss_rand_y1
} {
{
/gauss_rand_x1 rand gauss_rand_r mul 2.0 mul 1.0 sub def
/gauss_rand_x2 rand gauss_rand_r mul 2.0 mul 1.0 sub def
/gauss_rand_w
gauss_rand_x1 gauss_rand_x1 mul gauss_rand_x2 gauss_rand_x2 mul add def
gauss_rand_w 1.0 le {exit} if
} loop
/gauss_rand_w gauss_rand_w ln -2.0 mul gauss_rand_w div sqrt def
/gauss_rand_y1 gauss_rand_x1 gauss_rand_w mul def
/gauss_rand_y2 gauss_rand_x2 gauss_rand_w mul def
/gauss_rand_already true def
gauss_rand_y1
} ifelse
} def
gauss_rand pop % initialise, and discard the resulting zero
/grand { % usage: mean stddev grand
gauss_rand mul add
} bind def
/irand { % usage: n irand -> a random integer 0 .. n-1
rand exch round cvi mod
} bind def
/urand { % usage: urand -> a random number between 0 and 1
rand gauss_rand_r mul
} bind def
/randomget { % usage: [ an array ] randomget
dup length rand exch mod get
} bind def
/randomgetn { % usage: [ an array ] n randomgetn
10 dict begin [ /n /arr_in ] { exch def } forall
n arr_in length ge {
arr_in
} {
/already n dict def [
n {
{
/i rand arr_in length mod def
/candidate arr_in i get def
already i known not { already i true put candidate exit } if
} loop
} repeat
] } ifelse
} bind def
/rayleigh_rand { % usage: sigma rayleigh_rand
1.0 urand sub ln -2 mul sqrt mul
} bind def