QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
knuthuniformrng.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 int KnuthUniformRng::KK = 100;
26 const int KnuthUniformRng::LL = 37;
27 const int KnuthUniformRng::TT = 70;
28 const int KnuthUniformRng::QUALITY = 1009;
29
31 : ranf_arr_buf(QUALITY), ran_u(QUALITY) {
33 ranf_start(seed != 0 ? seed : SeedGenerator::instance().get());
34 }
35
37 int t,s,j;
38 std::vector<double> u(KK+KK-1),ul(KK+KK-1);
39 double ulp=(1.0/(1L<<30))/(1L<<22); // 2 to the -52
40 double ss=2.0*ulp*((seed&0x3fffffff)+2);
41
42 for (j=0;j<KK;j++) {
43 u[j]=ss; ul[j]=0.0; // bootstrap the buffer
44 ss+=ss; if (ss>=1.0) ss-=1.0-2*ulp; // cyclic shift of 51 bits
45 }
46 for (;j<KK+KK-1;j++) u[j]=ul[j]=0.0;
47 u[1]+=ulp;ul[1]=ulp; // make u[1] (and only u[1]) "odd"
48 s=seed&0x3fffffff;
49 t=TT-1;
50 while (t != 0) {
51 for (j=KK-1;j>0;--j) ul[j+j]=ul[j],u[j+j]=u[j]; // "square"
52 for (j=KK+KK-2;j>KK-LL;j-=2)
53 ul[KK+KK-1-j]=0.0,u[KK+KK-1-j]=u[j]-ul[j];
54 for (j=KK+KK-2;j>=KK;--j)
55 if (ul[j] != 0.0) {
56 ul[j - (KK - LL)] = ulp - ul[j - (KK - LL)],
57 u[j - (KK - LL)] = mod_sum(u[j - (KK - LL)], u[j]);
58 ul[j - KK] = ulp - ul[j - KK], u[j - KK] = mod_sum(u[j - KK], u[j]);
59 }
60 if (is_odd(s)) { // "multiply by z"
61 for (j=KK;j>0;--j) ul[j]=ul[j-1],u[j]=u[j-1];
62 ul[0]=ul[KK],u[0]=u[KK]; // shift the buffer cyclically
63 if (ul[KK] != 0.0)
64 ul[LL] = ulp - ul[LL], u[LL] = mod_sum(u[LL], u[KK]);
65 }
66 if (s != 0)
67 s >>= 1;
68 else
69 t--;
70 }
71 for (j=0;j<LL;j++) ran_u[j+KK-LL]=u[j];
72 for (;j<KK;j++) ran_u[j-LL]=u[j];
73 }
74
75 void KnuthUniformRng::ranf_array(std::vector<double>& aa,
76 int n) const {
77 int i,j;
78 for (j=0;j<KK;j++) aa[j]=ran_u[j];
79 for (;j<n;j++) aa[j]=mod_sum(aa[j-KK],aa[j-LL]);
80 for (i=0;i<LL;i++,j++) ran_u[i]=mod_sum(aa[j-KK],aa[j-LL]);
81 for (;i<KK;i++,j++) ran_u[i]=mod_sum(aa[j-KK],ran_u[i-LL]);
82 }
83
86 ranf_arr_ptr = 1;
88 return ranf_arr_buf[0];
89 }
90
91}
void ranf_array(std::vector< double > &aa, int n) const
double mod_sum(double x, double y) const
std::vector< double > ran_u
std::vector< double > ranf_arr_buf
static SeedGenerator & instance()
access to the unique instance
Definition: singleton.hpp:104
const DefaultType & t
Knuth uniform random number generator.
Definition: any.hpp:35
Random seed generator.