QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
zigguratgaussianrng.hpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2024 Ralf Konrad Eckel
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/*! \file zigguratgaussianrng.hpp
21 \brief Ziggurat Gaussian random-number generator
22*/
23
24#ifndef quantlib_ziggurat_gaussian_rng_h
25#define quantlib_ziggurat_gaussian_rng_h
26
28#include <cstdint>
29
30namespace QuantLib {
31
32 //! Gaussian random number generator
33 /*! It uses the Ziggurat transformation to return a
34 normal distributed Gaussian deviate with average 0.0 and
35 standard deviation of 1.0, from a random integer
36 in the [0,0xffffffffffffffffULL]-interval like.
37
38 For a more detailed description see the article
39 "An Improved Ziggurat Method to Generate Normal Random Samples"
40 by Jurgen A. Doornik
41 (https://www.doornik.com/research/ziggurat.pdf).
42
43 The code here is inspired by the rust implementation in
44 https://github.com/rust-random/rand/blob/d42daabf65a3ceaf58c2eefc7eb477c4d5a9b4ba/rand_distr/src/normal.rs
45 and
46 https://github.com/rust-random/rand/blob/d42daabf65a3ceaf58c2eefc7eb477c4d5a9b4ba/rand_distr/src/utils.rs.
47
48 Class RNG must implement the following interface:
49 \code
50 Real nextReal() const;
51 std::uint64_t nextInt64() const;
52 \endcode
53 Currently, Xoshiro256StarStarUniformRng is the only RNG supporting this.
54 */
55 template <class RNG>
57 public:
59
60 explicit ZigguratGaussianRng(const RNG& uint64Generator)
61 : uint64Generator_(uint64Generator) {}
62
63 //! returns a sample from a Gaussian distribution
64 sample_type next() const { return {nextReal(), 1.0}; }
65
66 //! return a random number from a Gaussian distribution
67 Real nextReal() const;
68
69 private:
71
72 typedef Real ZigguratTable[257];
73
74 Real pdf(Real x) const;
75
76 //! compute a random number in the tail by hand
77 Real zeroCase(Real u) const;
78
79 Real normR() const;
80 Real normX(int i) const;
81 Real normF(int i) const;
82 };
83
84 template <class RNG>
86 while (true) {
87 // As an optimisation we re-implement the conversion
88 // to a double in the interval (-1,1).
89 // From the remaining 12 most significant bits we use 8 to construct `i`.
90 //
91 // This saves us generating a whole extra random number, while the added
92 // precision of using 64 bits for double does not buy us much.
93 std::uint64_t randomU64 = uint64Generator_.nextInt64();
94 auto u = 2.0 * (Real(randomU64 >> 11) + 0.5) * (1.0 / Real(1ULL << 53)) - 1.0;
95 auto i = (int)(randomU64 & 0xff);
96
97 auto x = u * normX(i);
98
99 if (std::abs(x) < normX(i + 1)) {
100 return x;
101 }
102 if (i == 0) {
103 // compute a random number in the tail by hand
104 return zeroCase(u);
105 }
106 if (normF(i + 1) + (normF(i) - normF(i + 1) * uint64Generator_.nextReal()) < pdf(x)) {
107 return x;
108 }
109 }
110 }
111
112 template <class RNG>
114 // compute a random number in the tail by hand
115 Real x, y;
116 do {
117 x = std::log(uint64Generator_.nextReal()) / normR();
118 y = std::log(uint64Generator_.nextReal());
119 } while (-2.0 * y < x * x);
120
121 return (u < 0.0) ? x - normR() : normR() - x;
122 }
123
124 template <class RNG>
126 return std::exp(-x * x / 2.0);
127 }
128
129 template <class RNG>
131 return 3.654152885361008796;
132 }
133
134 template <class RNG>
136 static const ZigguratTable normX = {
137 3.910757959537090045, 3.654152885361008796, 3.449278298560964462, 3.320244733839166074,
138 3.224575052047029100, 3.147889289517149969, 3.083526132001233044, 3.027837791768635434,
139 2.978603279880844834, 2.934366867207854224, 2.894121053612348060, 2.857138730872132548,
140 2.822877396825325125, 2.790921174000785765, 2.760944005278822555, 2.732685359042827056,
141 2.705933656121858100, 2.680514643284522158, 2.656283037575502437, 2.633116393630324570,
142 2.610910518487548515, 2.589575986706995181, 2.569035452680536569, 2.549221550323460761,
143 2.530075232158516929, 2.511544441625342294, 2.493583041269680667, 2.476149939669143318,
144 2.459208374333311298, 2.442725318198956774, 2.426670984935725972, 2.411018413899685520,
145 2.395743119780480601, 2.380822795170626005, 2.366237056715818632, 2.351967227377659952,
146 2.337996148795031370, 2.324308018869623016, 2.310888250599850036, 2.297723348901329565,
147 2.284800802722946056, 2.272108990226823888, 2.259637095172217780, 2.247375032945807760,
148 2.235313384928327984, 2.223443340090905718, 2.211756642882544366, 2.200245546609647995,
149 2.188902771624720689, 2.177721467738641614, 2.166695180352645966, 2.155817819875063268,
150 2.145083634046203613, 2.134487182844320152, 2.124023315687815661, 2.113687150684933957,
151 2.103474055713146829, 2.093379631137050279, 2.083399693996551783, 2.073530263516978778,
152 2.063767547809956415, 2.054107931648864849, 2.044547965215732788, 2.035084353727808715,
153 2.025713947862032960, 2.016433734904371722, 2.007240830558684852, 1.998132471356564244,
154 1.989106007615571325, 1.980158896898598364, 1.971288697931769640, 1.962493064942461896,
155 1.953769742382734043, 1.945116560006753925, 1.936531428273758904, 1.928012334050718257,
156 1.919557336591228847, 1.911164563769282232, 1.902832208548446369, 1.894558525668710081,
157 1.886341828534776388, 1.878180486290977669, 1.870072921069236838, 1.862017605397632281,
158 1.854013059758148119, 1.846057850283119750, 1.838150586580728607, 1.830289919680666566,
159 1.822474540091783224, 1.814703175964167636, 1.806974591348693426, 1.799287584547580199,
160 1.791640986550010028, 1.784033659547276329, 1.776464495522344977, 1.768932414909077933,
161 1.761436365316706665, 1.753975320315455111, 1.746548278279492994, 1.739154261283669012,
162 1.731792314050707216, 1.724461502945775715, 1.717160915015540690, 1.709889657069006086,
163 1.702646854797613907, 1.695431651932238548, 1.688243209434858727, 1.681080704722823338,
164 1.673943330923760353, 1.666830296159286684, 1.659740822855789499, 1.652674147080648526,
165 1.645629517902360339, 1.638606196773111146, 1.631603456932422036, 1.624620582830568427,
166 1.617656869570534228, 1.610711622367333673, 1.603784156023583041, 1.596873794420261339,
167 1.589979870021648534, 1.583101723393471438, 1.576238702733332886, 1.569390163412534456,
168 1.562555467528439657, 1.555733983466554893, 1.548925085471535512, 1.542128153226347553,
169 1.535342571438843118, 1.528567729435024614, 1.521803020758293101, 1.515047842773992404,
170 1.508301596278571965, 1.501563685112706548, 1.494833515777718391, 1.488110497054654369,
171 1.481394039625375747, 1.474683555695025516, 1.467978458615230908, 1.461278162507407830,
172 1.454582081885523293, 1.447889631277669675, 1.441200224845798017, 1.434513276002946425,
173 1.427828197027290358, 1.421144398672323117, 1.414461289772464658, 1.407778276843371534,
174 1.401094763676202559, 1.394410150925071257, 1.387723835686884621, 1.381035211072741964,
175 1.374343665770030531, 1.367648583594317957, 1.360949343030101844, 1.354245316759430606,
176 1.347535871177359290, 1.340820365893152122, 1.334098153216083604, 1.327368577624624679,
177 1.320630975217730096, 1.313884673146868964, 1.307128989027353860, 1.300363230327433728,
178 1.293586693733517645, 1.286798664489786415, 1.279998415710333237, 1.273185207661843732,
179 1.266358287014688333, 1.259516886060144225, 1.252660221891297887, 1.245787495544997903,
180 1.238897891102027415, 1.231990574742445110, 1.225064693752808020, 1.218119375481726552,
181 1.211153726239911244, 1.204166830140560140, 1.197157747875585931, 1.190125515422801650,
182 1.183069142678760732, 1.175987612011489825, 1.168879876726833800, 1.161744859441574240,
183 1.154581450355851802, 1.147388505416733873, 1.140164844363995789, 1.132909248648336975,
184 1.125620459211294389, 1.118297174115062909, 1.110938046009249502, 1.103541679420268151,
185 1.096106627847603487, 1.088631390649514197, 1.081114409698889389, 1.073554065787871714,
186 1.065948674757506653, 1.058296483326006454, 1.050595664586207123, 1.042844313139370538,
187 1.035040439828605274, 1.027181966030751292, 1.019266717460529215, 1.011292417434978441,
188 1.003256679539591412, 0.995156999629943084, 0.986990747093846266, 0.978755155288937750,
189 0.970447311058864615, 0.962064143217605250, 0.953602409875572654, 0.945058684462571130,
190 0.936429340280896860, 0.927710533396234771, 0.918898183643734989, 0.909987953490768997,
191 0.900975224455174528, 0.891855070726792376, 0.882622229578910122, 0.873271068082494550,
192 0.863795545546826915, 0.854189171001560554, 0.844444954902423661, 0.834555354079518752,
193 0.824512208745288633, 0.814306670128064347, 0.803929116982664893, 0.793369058833152785,
194 0.782615023299588763, 0.771654424216739354, 0.760473406422083165, 0.749056662009581653,
195 0.737387211425838629, 0.725446140901303549, 0.713212285182022732, 0.700661841097584448,
196 0.687767892786257717, 0.674499822827436479, 0.660822574234205984, 0.646695714884388928,
197 0.632072236375024632, 0.616896989996235545, 0.601104617743940417, 0.584616766093722262,
198 0.567338257040473026, 0.549151702313026790, 0.529909720646495108, 0.509423329585933393,
199 0.487443966121754335, 0.463634336771763245, 0.437518402186662658, 0.408389134588000746,
200 0.375121332850465727, 0.335737519180459465, 0.286174591747260509, 0.215241895913273806,
201 0.000000000000000000};
202 return normX[i];
203 }
204
205 template <class RNG>
207 static const ZigguratGaussianRng<RNG>::ZigguratTable normF = {
208 0.000477467764586655, 0.001260285930498598, 0.002609072746106363, 0.004037972593371872,
209 0.005522403299264754, 0.007050875471392110, 0.008616582769422917, 0.010214971439731100,
210 0.011842757857943104, 0.013497450601780807, 0.015177088307982072, 0.016880083152595839,
211 0.018605121275783350, 0.020351096230109354, 0.022117062707379922, 0.023902203305873237,
212 0.025705804008632656, 0.027527235669693315, 0.029365939758230111, 0.031221417192023690,
213 0.033093219458688698, 0.034980941461833073, 0.036884215688691151, 0.038802707404656918,
214 0.040736110656078753, 0.042684144916619378, 0.044646552251446536, 0.046623094902089664,
215 0.048613553216035145, 0.050617723861121788, 0.052635418276973649, 0.054666461325077916,
216 0.056710690106399467, 0.058767952921137984, 0.060838108349751806, 0.062921024437977854,
217 0.065016577971470438, 0.067124653828023989, 0.069245144397250269, 0.071377949059141965,
218 0.073522973714240991, 0.075680130359194964, 0.077849336702372207, 0.080030515814947509,
219 0.082223595813495684, 0.084428509570654661, 0.086645194450867782, 0.088873592068594229,
220 0.091113648066700734, 0.093365311913026619, 0.095628536713353335, 0.097903279039215627,
221 0.100189498769172020, 0.102487158942306270, 0.104796225622867056, 0.107116667775072880,
222 0.109448457147210021, 0.111791568164245583, 0.114145977828255210, 0.116511665626037014,
223 0.118888613443345698, 0.121276805485235437, 0.123676228202051403, 0.126086870220650349,
224 0.128508722280473636, 0.130941777174128166, 0.133386029692162844, 0.135841476571757352,
225 0.138308116449064322, 0.140785949814968309, 0.143274978974047118, 0.145775208006537926,
226 0.148286642733128721, 0.150809290682410169, 0.153343161060837674, 0.155888264725064563,
227 0.158444614156520225, 0.161012223438117663, 0.163591108232982951, 0.166181285765110071,
228 0.168782774801850333, 0.171395595638155623, 0.174019770082499359, 0.176655321444406654,
229 0.179302274523530397, 0.181960655600216487, 0.184630492427504539, 0.187311814224516926,
230 0.190004651671193070, 0.192709036904328807, 0.195425003514885592, 0.198152586546538112,
231 0.200891822495431333, 0.203642749311121501, 0.206405406398679298, 0.209179834621935651,
232 0.211966076307852941, 0.214764175252008499, 0.217574176725178370, 0.220396127481011589,
233 0.223230075764789593, 0.226076071323264877, 0.228934165415577484, 0.231804410825248525,
234 0.234686861873252689, 0.237581574432173676, 0.240488605941449107, 0.243408015423711988,
235 0.246339863502238771, 0.249284212419516704, 0.252241126056943765, 0.255210669955677150,
236 0.258192911338648023, 0.261187919133763713, 0.264195763998317568, 0.267216518344631837,
237 0.270250256366959984, 0.273297054069675804, 0.276356989296781264, 0.279430141762765316,
238 0.282516593084849388, 0.285616426816658109, 0.288729728483353931, 0.291856585618280984,
239 0.294997087801162572, 0.298151326697901342, 0.301319396102034120, 0.304501391977896274,
240 0.307697412505553769, 0.310907558127563710, 0.314131931597630143, 0.317370638031222396,
241 0.320623784958230129, 0.323891482377732021, 0.327173842814958593, 0.330470981380537099,
242 0.333783015832108509, 0.337110066638412809, 0.340452257045945450, 0.343809713148291340,
243 0.347182563958251478, 0.350570941482881204, 0.353974980801569250, 0.357394820147290515,
244 0.360830600991175754, 0.364282468130549597, 0.367750569780596226, 0.371235057669821344,
245 0.374736087139491414, 0.378253817247238111, 0.381788410875031348, 0.385340034841733958,
246 0.388908860020464597, 0.392495061461010764, 0.396098818517547080, 0.399720314981931668,
247 0.403359739222868885, 0.407017284331247953, 0.410693148271983222, 0.414387534042706784,
248 0.418100649839684591, 0.421832709231353298, 0.425583931339900579, 0.429354541031341519,
249 0.433144769114574058, 0.436954852549929273, 0.440785034667769915, 0.444635565397727750,
250 0.448506701509214067, 0.452398706863882505, 0.456311852680773566, 0.460246417814923481,
251 0.464202689050278838, 0.468180961407822172, 0.472181538469883255, 0.476204732721683788,
252 0.480250865911249714, 0.484320269428911598, 0.488413284707712059, 0.492530263646148658,
253 0.496671569054796314, 0.500837575128482149, 0.505028667945828791, 0.509245245998136142,
254 0.513487720749743026, 0.517756517232200619, 0.522052074674794864, 0.526374847174186700,
255 0.530725304406193921, 0.535103932383019565, 0.539511234259544614, 0.543947731192649941,
256 0.548413963257921133, 0.552910490428519918, 0.557437893621486324, 0.561996775817277916,
257 0.566587763258951771, 0.571211506738074970, 0.575868682975210544, 0.580559996103683473,
258 0.585286179266300333, 0.590047996335791969, 0.594846243770991268, 0.599681752622167719,
259 0.604555390700549533, 0.609468064928895381, 0.614420723892076803, 0.619414360609039205,
260 0.624450015550274240, 0.629528779928128279, 0.634651799290960050, 0.639820277456438991,
261 0.645035480824251883, 0.650298743114294586, 0.655611470583224665, 0.660975147780241357,
262 0.666391343912380640, 0.671861719900766374, 0.677388036222513090, 0.682972161648791376,
263 0.688616083008527058, 0.694321916130032579, 0.700091918140490099, 0.705928501336797409,
264 0.711834248882358467, 0.717811932634901395, 0.723864533472881599, 0.729995264565802437,
265 0.736207598131266683, 0.742505296344636245, 0.748892447223726720, 0.755373506511754500,
266 0.761953346841546475, 0.768637315803334831, 0.775431304986138326, 0.782341832659861902,
267 0.789376143571198563, 0.796542330428254619, 0.803849483176389490, 0.811307874318219935,
268 0.818929191609414797, 0.826726833952094231, 0.834716292992930375, 0.842915653118441077,
269 0.851346258465123684, 0.860033621203008636, 0.869008688043793165, 0.878309655816146839,
270 0.887984660763399880, 0.898095921906304051, 0.908726440060562912, 0.919991505048360247,
271 0.932060075968990209, 0.945198953453078028, 0.959879091812415930, 0.977101701282731328,
272 1.000000000000000000};
273 return normF[i];
274 }
275}
276
277#endif
Gaussian random number generator.
sample_type next() const
returns a sample from a Gaussian distribution
Real nextReal() const
return a random number from a Gaussian distribution
ZigguratGaussianRng(const RNG &uint64Generator)
Real zeroCase(Real u) const
compute a random number in the tail by hand
QL_REAL Real
real number
Definition: types.hpp:50
Definition: any.hpp:35
weighted sample
weighted sample
Definition: sample.hpp:35