QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
lecuyeruniformrng.cpp
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) 2000, 2001, 2002, 2003 RiskMap srl
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
22
23namespace QuantLib {
24
25 const long LecuyerUniformRng::m1 = 2147483563L;
26 const long LecuyerUniformRng::a1 = 40014L;
27 const long LecuyerUniformRng::q1 = 53668L;
28 const long LecuyerUniformRng::r1 = 12211L;
29
30 const long LecuyerUniformRng::m2 = 2147483399L;
31 const long LecuyerUniformRng::a2 = 40692L;
32 const long LecuyerUniformRng::q2 = 52774L;
33 const long LecuyerUniformRng::r2 = 3791L;
34
35 const int LecuyerUniformRng::bufferSize = 32;
36
37 // int(1+m1/bufferSize) = int(1+(m1-1)/bufferSize)
38 const long LecuyerUniformRng::bufferNormalizer = 67108862L;
39
40 const long double LecuyerUniformRng::maxRandom = 1.0-QL_EPSILON;
41
43 : buffer(LecuyerUniformRng::bufferSize) {
44 // Need to prevent seed=0, so use seed=0 to have a "random" seed
45 temp2 = temp1 = (seed != 0 ? seed : SeedGenerator::instance().get());
46 // Load the shuffle table (after 8 warm-ups)
47 for (int j=bufferSize+7; j>=0; j--) {
48 long k = temp1/q1;
49 temp1 = a1*(temp1-k*q1)-k*r1;
50 if (temp1 < 0)
51 temp1 += m1;
52 if (j < bufferSize)
53 buffer[j] = temp1;
54 }
55 y = buffer[0];
56 }
57
59 long k = temp1/q1;
60 // Compute temp1=(a1*temp1) % m1
61 // without overflows (Schrage's method)
62 temp1 = a1*(temp1-k*q1)-k*r1;
63 if (temp1 < 0)
64 temp1 += m1;
65 k = temp2/q2;
66 // Compute temp2=(a2*temp2) % m2
67 // without overflows (Schrage's method)
68 temp2 = a2*(temp2-k*q2)-k*r2;
69 if (temp2 < 0)
70 temp2 += m2;
71 // Will be in the range 0..bufferSize-1
72 int j = y/bufferNormalizer;
73 // Here temp1 is shuffled, temp1 and temp2 are
74 // combined to generate output
75 y = buffer[j]-temp2;
76 buffer[j] = temp1;
77 if (y < 1)
78 y += m1-1;
79 double result = y/double(m1);
80 // users don't expect endpoint values
81 if (result > maxRandom)
82 result = (double) maxRandom;
83 return {result, 1.0};
84 }
85
86}
Uniform random number generator.
static const long double maxRandom
static const long bufferNormalizer
static SeedGenerator & instance()
access to the unique instance
Definition: singleton.hpp:104
#define QL_EPSILON
Definition: qldefines.hpp:178
L'Ecuyer uniform random number generator.
Definition: any.hpp:35
Random seed generator.
weighted sample
Definition: sample.hpp:35