1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2010 Kakhkhor Abdijalilov
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include <ql/experimental/math/zigguratrng.hpp>
21#include <ql/math/distributions/normaldistribution.hpp>
22#include <cmath>
23
24namespace QuantLib {
25
26 namespace {
27
28 // tail probability
29 const Real p_ = 2.880541027242713E-004;
30 const Real q_ = 1.0 - p_;
31
32 /* The tabulated values were calculated following Marsaglia
33 and Tsang (2000). */
34
35 // values of exp(-0.5*x*x)
36 const Real f_ [128] = {
37 1.000000000000000E+000, 9.635996931557717E-001,
38 9.362826817083744E-001, 9.130436479920410E-001,
39 8.922816508023054E-001, 8.732430489268560E-001,
40 8.555006078850665E-001, 8.387836053106493E-001,
41 8.229072113952640E-001, 8.077382946961230E-001,
42 7.931770117838610E-001, 7.791460859417049E-001,
43 7.655841739092376E-001, 7.524415591857053E-001,
44 7.396772436833397E-001, 7.272569183545073E-001,
45 7.151515074204785E-001, 7.033360990258188E-001,
46 6.917891434460373E-001, 6.804918410064157E-001,
47 6.694276673577075E-001, 6.585820000586550E-001,
48 6.479418211185520E-001, 6.374954773431460E-001,
49 6.272324852578157E-001, 6.171433708265636E-001,
50 6.072195366326060E-001, 5.974531509518134E-001,
51 5.878370544418217E-001, 5.783646811267034E-001,
52 5.690299910747226E-001, 5.598274127106959E-001,
53 5.507517931210564E-001, 5.417983550317252E-001,
54 5.329626593899887E-001, 5.242405726789938E-001,
55 5.156282382498731E-001, 5.071220510813057E-001,
56 4.987186354765854E-001, 4.904148252893227E-001,
57 4.822076463348397E-001, 4.740943006982505E-001,
58 4.660721526945719E-001, 4.581387162728729E-001,
59 4.502916436869279E-001, 4.425287152802475E-001,
60 4.348478302546628E-001, 4.272469983095633E-001,
61 4.197243320540391E-001, 4.122780401070255E-001,
62 4.049064208114891E-001, 3.976078564980433E-001,
63 3.903808082413902E-001, 3.832238110598844E-001,
64 3.761354695144552E-001, 3.691144536682758E-001,
65 3.621594953730338E-001, 3.552693848515477E-001,
66 3.484429675498729E-001, 3.416791412350141E-001,
67 3.349768533169716E-001, 3.283350983761528E-001,
68 3.217529158792090E-001, 3.152293880681579E-001,
69 3.087636380092523E-001, 3.023548277894802E-001,
70 2.960021568498564E-001, 2.897048604458110E-001,
71 2.834622082260129E-001, 2.772735029218981E-001,
72 2.711380791410257E-001, 2.650553022581624E-001,
73 2.590245673987112E-001, 2.530452985097663E-001,
74 2.471169475146971E-001, 2.412389935477517E-001,
75 2.354109422657280E-001, 2.296323252343031E-001,
76 2.239026993871343E-001, 2.182216465563709E-001,
77 2.125887730737364E-001, 2.070037094418741E-001,
78 2.014661100762035E-001, 1.959756531181106E-001,
79 1.905320403209139E-001, 1.851349970107136E-001,
80 1.797842721249623E-001, 1.744796383324025E-001,
81 1.692208922389250E-001, 1.640078546849280E-001,
82 1.588403711409353E-001, 1.537183122095867E-001,
83 1.486415742436971E-001, 1.436100800919331E-001,
84 1.386237799858510E-001, 1.336826525846477E-001,
85 1.287867061971040E-001, 1.239359802039816E-001,
86 1.191305467087186E-001, 1.143705124498883E-001,
87 1.096560210158178E-001, 1.049872554103546E-001,
88 1.003644410295456E-001, 9.578784912257826E-002,
89 9.125780082763474E-002, 8.677467189554304E-002,
90 8.233889824295743E-002, 7.795098251465470E-002,
91 7.361150188475492E-002, 6.932111739418027E-002,
92 6.508058521363191E-002, 6.089077034856640E-002,
93 5.675266348153862E-002, 5.266740190350321E-002,
94 4.863629586028410E-002, 4.466086220087247E-002,
95 4.074286807479065E-002, 3.688438878696881E-002,
96 3.308788614650520E-002, 2.935631744025387E-002,
97 2.569329193614964E-002, 2.210330461611161E-002,
98 1.859210273716583E-002, 1.516729801067205E-002,
99 1.183947865798232E-002, 8.624484412930473E-003,
100 5.548995220816476E-003, 2.669629083902507E-003
101 };
102
103 // acceptance thresholds 2^24*x[i]/x[i+1]. k_[0] is special
104 const Size k_[128] = {
105 15555141, 0, 12590647, 14272656,
106 14988942, 15384587, 15635012, 15807564,
107 15933580, 16029597, 16105158, 16166150,
108 16216402, 16258511, 16294298, 16325081,
109 16351834, 16375294, 16396029, 16414482,
110 16431005, 16445883, 16459346, 16471581,
111 16482747, 16492974, 16502372, 16511034,
112 16519042, 16526462, 16533356, 16539772,
113 16545758, 16551351, 16556587, 16561496,
114 16566104, 16570437, 16574515, 16578357,
115 16581980, 16585401, 16588633, 16591688,
116 16594579, 16597314, 16599905, 16602358,
117 16604682, 16606885, 16608972, 16610949,
118 16612822, 16614597, 16616276, 16617865,
119 16619367, 16620786, 16622125, 16623387,
120 16624575, 16625690, 16626735, 16627713,
121 16628624, 16629470, 16630253, 16630974,
122 16631634, 16632233, 16632773, 16633254,
123 16633677, 16634041, 16634346, 16634593,
124 16634781, 16634910, 16634979, 16634987,
125 16634934, 16634817, 16634637, 16634390,
126 16634075, 16633689, 16633231, 16632698,
127 16632085, 16631390, 16630609, 16629737,
128 16628768, 16627698, 16626520, 16625226,
129 16623808, 16622257, 16620563, 16618714,
130 16616696, 16614494, 16612091, 16609465,
131 16606593, 16603449, 16599999, 16596206,
132 16592025, 16587402, 16582273, 16576559,
133 16570163, 16562965, 16554812, 16545511,
134 16534809, 16522368, 16507733, 16490265,
135 16469045, 16442690, 16409026, 16364394,
136 16302111, 16208408, 16049219, 15707338
137 };
138
139 // values of 2^{-24}*x[i]. w_[0] is special.
140 const double w_[128] = {
141 2.213171867573477E-007, 1.623158840564778E-008,
142 2.162882274558596E-008, 2.542424120326624E-008,
143 2.845751269184242E-008, 3.103351823837397E-008,
144 3.330064883086164E-008, 3.534334554922425E-008,
145 3.721467240506913E-008, 3.895036212891571E-008,
146 4.057573787247544E-008, 4.210946627340346E-008,
147 4.356574479471913E-008, 4.495565083232566E-008,
148 4.628801273561392E-008, 4.756999377168848E-008,
149 4.880749623079987E-008, 5.000544871575862E-008,
150 5.116801519263080E-008, 5.229875022755345E-008,
151 5.340071633852936E-008, 5.447657412343023E-008,
152 5.552865246542405E-008, 5.655900391923845E-008,
153 5.756944891143612E-008, 5.856161138431779E-008,
154 5.953694781545649E-008, 6.049677105184184E-008,
155 6.144227004387700E-008, 6.237452630714050E-008,
156 6.329452775023089E-008, 6.420318036567782E-008,
157 6.510131817439508E-008, 6.598971173307500E-008,
158 6.686907545162751E-008, 6.774007391947947E-008,
159 6.860332740181531E-008, 6.945941663712532E-008,
160 7.030888704386109E-008, 7.115225242518010E-008,
161 7.198999824564194E-008, 7.282258454149729E-008,
162 7.365044851627824E-008, 7.447400686528278E-008,
163 7.529365786588351E-008, 7.610978326509584E-008,
164 7.692274999129007E-008, 7.773291171314836E-008,
165 7.854061026581177E-008, 7.934617696152180E-008,
166 8.014993379984568E-008, 8.095219459071287E-008,
167 8.175326600192373E-008, 8.255344854147119E-008,
168 8.335303748390705E-008, 8.415232374905104E-008,
169 8.495159474056128E-008, 8.575113515123489E-008,
170 8.655122774137352E-008, 8.735215409611426E-008,
171 8.815419536728245E-008, 8.895763300505963E-008,
172 8.976274948457178E-008, 9.056982903238356E-008,
173 9.137915835783214E-008, 9.219102739414587E-008,
174 9.300573005436895E-008, 9.382356500725440E-008,
175 9.464483647849558E-008, 9.546985508294559E-008,
176 9.629893869382930E-008, 9.713241335539087E-008,
177 9.797061424595009E-008, 9.881388669897357E-008,
178 9.966258729051657E-008, 1.005170850022725E-007,
179 1.013777624705017E-007, 1.022450173323223E-007,
180 1.031192636822607E-007, 1.040009336536155E-007,
181 1.048904791411299E-007, 1.057883736837368E-007,
182 1.066951145288121E-007, 1.076112249025135E-007,
183 1.085372565144899E-007, 1.094737923296323E-007,
184 1.104214496447496E-007, 1.113808835142578E-007,
185 1.123527905763905E-007, 1.133379133403490E-007,
186 1.143370450055439E-007, 1.153510348970830E-007,
187 1.163807946174674E-007, 1.174273050337859E-007,
188 1.184916242434419E-007, 1.195748966907839E-007,
189 1.206783636434635E-007, 1.218033752829236E-007,
190 1.229514047207811E-007, 1.241240643255547E-007,
191 1.253231248369812E-007, 1.265505378645533E-007,
192 1.278084625218070E-007, 1.290992971506620E-007,
193 1.304257173581136E-007, 1.317907219454484E-007,
194 1.331976887933646E-007, 1.346504434266883E-007,
195 1.361533438964878E-007, 1.377113869008423E-007,
196 1.393303418955523E-007, 1.410169225999109E-007,
197 1.427790092234294E-007, 1.446259406525023E-007,
198 1.465689049606532E-007, 1.486214710528821E-007,
199 1.508003278008381E-007, 1.531263366890930E-007,
200 1.556260733859904E-007, 1.583341605221148E-007,
201 1.612969382476045E-007, 1.645785196056458E-007,
202 1.682713836756925E-007, 1.725163463961286E-007,
203 1.775441320326934E-007, 1.837747608550914E-007,
204 1.921108355867039E-007, 2.051961336074264E-007
205 };
206
207 }
208
209 ZigguratRng::ZigguratRng(unsigned long seed)
210 : mt32_(seed) {}
211
212 Real ZigguratRng::nextGaussian() const {
213 static const int c[2] = {-1, 1};
214 Real x;
215
216 for (;;) {
217 unsigned long j = mt32_.nextInt32(); // generate 32 bits of randomness
218 int f = j & 1; // 1 bit to choose a tails
219 j >>= 1;
220 unsigned long i = j & 0x7f; // 7 bits to choose a strip
221 j >>= 7; // the last 24 bits for accepttion/rejection
222 x = (c[f]*static_cast<long>(j))*w_[i]; // x is uniform
223 // within the i-th strip
224 if (j < k_[i]) // if true, accept x
225 break;
226
227 // handle rejections
228 if (i!=0) { // upper strips
229 if ((f_[i-1]-f_[i])*mt32_.nextReal() + f_[i] < std::exp(x: -0.5*x*x))
230 break;
231 } else { // base strip, sample from the tail
232 x = c[f]*InverseCumulativeNormal::standard_value(
233 x: p_*mt32_.nextReal()+q_);
234 break;
235 }
236 }
237
238 return x;
239 }
240
241}
242

source code of quantlib/ql/experimental/math/zigguratrng.cpp