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
88
89
90
91
92
93
94
|
/*
Copyright (C) 2000, 2001, 2002 RiskMap srl
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
QuantLib is free software: you can redistribute it and/or modify it under the
terms of the QuantLib license. You should have received a copy of the
license along with this program; if not, please email ferdinando@ametrano.net
The license is also available online at http://quantlib.org/html/license.html
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
*/
/*! \file knuthuniformrng.cpp
\brief Knuth uniform random number generator
\fullpath
ql/RandomNumbers/%knuthuniformrng.cpp
*/
// $Id: knuthuniformrng.cpp,v 1.4 2002/01/16 14:41:27 nando Exp $
#include <ql/RandomNumbers/knuthuniformrng.hpp>
namespace QuantLib {
namespace RandomNumbers {
const int KnuthUniformRng::KK = 100;
const int KnuthUniformRng::LL = 37;
const int KnuthUniformRng::TT = 70;
const int KnuthUniformRng::QUALITY = 1009;
KnuthUniformRng::KnuthUniformRng(long seed)
: ranf_arr_buf(QUALITY), ran_u(QUALITY) {
ranf_arr_ptr = ranf_arr_sentinel = ranf_arr_buf.end();
ranf_start(seed != 0 ? seed : long(QL_TIME(0)));
}
void KnuthUniformRng::ranf_start(long seed) {
int t,s,j;
std::vector<double> u(KK+KK-1),ul(KK+KK-1);
double ulp=(1.0/(1L<<30))/(1L<<22); // 2 to the -52
double ss=2.0*ulp*((seed&0x3fffffff)+2);
for (j=0;j<KK;j++) {
u[j]=ss; ul[j]=0.0; // bootstrap the buffer
ss+=ss; if (ss>=1.0) ss-=1.0-2*ulp; // cyclic shift of 51 bits
}
for (;j<KK+KK-1;j++) u[j]=ul[j]=0.0;
u[1]+=ulp;ul[1]=ulp; // make u[1] (and only u[1]) "odd"
s=seed&0x3fffffff;
t=TT-1; while (t) {
for (j=KK-1;j>0;j--) ul[j+j]=ul[j],u[j+j]=u[j]; // "square"
for (j=KK+KK-2;j>KK-LL;j-=2)
ul[KK+KK-1-j]=0.0,u[KK+KK-1-j]=u[j]-ul[j];
for (j=KK+KK-2;j>=KK;j--) if(ul[j]) {
ul[j-(KK-LL)]=ulp-ul[j-(KK-LL)],
u[j-(KK-LL)]=mod_sum(u[j-(KK-LL)],u[j]);
ul[j-KK]=ulp-ul[j-KK],u[j-KK]=mod_sum(u[j-KK],u[j]);
}
if (is_odd(s)) { // "multiply by z"
for (j=KK;j>0;j--) ul[j]=ul[j-1],u[j]=u[j-1];
ul[0]=ul[KK],u[0]=u[KK]; // shift the buffer cyclically
if (ul[KK]) ul[LL]=ulp-ul[LL],u[LL]=mod_sum(u[LL],u[KK]);
}
if (s) s>>=1; else t--;
}
for (j=0;j<LL;j++) ran_u[j+KK-LL]=u[j];
for (;j<KK;j++) ran_u[j-LL]=u[j];
}
void KnuthUniformRng::ranf_array(std::vector<double>& aa,
int n) const {
int i,j;
for (j=0;j<KK;j++) aa[j]=ran_u[j];
for (;j<n;j++) aa[j]=mod_sum(aa[j-KK],aa[j-LL]);
for (i=0;i<LL;i++,j++) ran_u[i]=mod_sum(aa[j-KK],aa[j-LL]);
for (;i<KK;i++,j++) ran_u[i]=mod_sum(aa[j-KK],ran_u[i-LL]);
}
double KnuthUniformRng::ranf_arr_cycle() const {
ranf_array(ranf_arr_buf,QUALITY);
ranf_arr_ptr=ranf_arr_buf.begin()+1;
return ranf_arr_buf[0];
}
}
}
|
