27 const Real resultAsc) {
28 err = std::fabs(err) ;
29 if (resultAsc != 0 && err != 0){
30 Real scale = std::pow((200 * err / resultAsc), 1.5) ;
32 err = resultAsc * scale ;
50 static const Real x1[5] = {
51 0.973906528517171720077964012084452,
52 0.865063366688984510732096688423493,
53 0.679409568299024406234327365114874,
54 0.433395394129247190799265943165784,
55 0.148874338981631210884826001129720
59 static const Real w10[5] = {
60 0.066671344308688137593568809893332,
61 0.149451349150580593145776339657697,
62 0.219086362515982043995534934228163,
63 0.269266719309996355091226921569469,
64 0.295524224714752870173892994651338
68 static const Real x2[5] = {
69 0.995657163025808080735527280689003,
70 0.930157491355708226001207180059508,
71 0.780817726586416897063717578345042,
72 0.562757134668604683339000099272694,
73 0.294392862701460198131126603103866
77 static const Real w21a[5] = {
78 0.032558162307964727478818972459390,
79 0.075039674810919952767043140916190,
80 0.109387158802297641899210590325805,
81 0.134709217311473325928054001771707,
82 0.147739104901338491374841515972068
86 static const Real w21b[6] = {
87 0.011694638867371874278064396062192,
88 0.054755896574351996031381300244580,
89 0.093125454583697605535065465083366,
90 0.123491976262065851077958109831074,
91 0.142775938577060080797094273138717,
92 0.149445554002916905664936468389821
96 static const Real x3[11] = {
97 0.999333360901932081394099323919911,
98 0.987433402908088869795961478381209,
99 0.954807934814266299257919200290473,
100 0.900148695748328293625099494069092,
101 0.825198314983114150847066732588520,
102 0.732148388989304982612354848755461,
103 0.622847970537725238641159120344323,
104 0.499479574071056499952214885499755,
105 0.364901661346580768043989548502644,
106 0.222254919776601296498260928066212,
107 0.074650617461383322043914435796506
111 static const Real w43a[10] = {
112 0.016296734289666564924281974617663,
113 0.037522876120869501461613795898115,
114 0.054694902058255442147212685465005,
115 0.067355414609478086075553166302174,
116 0.073870199632393953432140695251367,
117 0.005768556059769796184184327908655,
118 0.027371890593248842081276069289151,
119 0.046560826910428830743339154433824,
120 0.061744995201442564496240336030883,
121 0.071387267268693397768559114425516
125 static const Real w43b[12] = {
126 0.001844477640212414100389106552965,
127 0.010798689585891651740465406741293,
128 0.021895363867795428102523123075149,
129 0.032597463975345689443882222526137,
130 0.042163137935191811847627924327955,
131 0.050741939600184577780189020092084,
132 0.058379395542619248375475369330206,
133 0.064746404951445885544689259517511,
134 0.069566197912356484528633315038405,
135 0.072824441471833208150939535192842,
136 0.074507751014175118273571813842889,
137 0.074722147517403005594425168280423
141 static const Real x4[22] = {
142 0.999902977262729234490529830591582,
143 0.997989895986678745427496322365960,
144 0.992175497860687222808523352251425,
145 0.981358163572712773571916941623894,
146 0.965057623858384619128284110607926,
147 0.943167613133670596816416634507426,
148 0.915806414685507209591826430720050,
149 0.883221657771316501372117548744163,
150 0.845710748462415666605902011504855,
151 0.803557658035230982788739474980964,
152 0.757005730685495558328942793432020,
153 0.706273209787321819824094274740840,
154 0.651589466501177922534422205016736,
155 0.593223374057961088875273770349144,
156 0.531493605970831932285268948562671,
157 0.466763623042022844871966781659270,
158 0.399424847859218804732101665817923,
159 0.329874877106188288265053371824597,
160 0.258503559202161551802280975429025,
161 0.185695396568346652015917141167606,
162 0.111842213179907468172398359241362,
163 0.037352123394619870814998165437704
167 static const Real w87a[21] = {
168 0.008148377384149172900002878448190,
169 0.018761438201562822243935059003794,
170 0.027347451050052286161582829741283,
171 0.033677707311637930046581056957588,
172 0.036935099820427907614589586742499,
173 0.002884872430211530501334156248695,
174 0.013685946022712701888950035273128,
175 0.023280413502888311123409291030404,
176 0.030872497611713358675466394126442,
177 0.035693633639418770719351355457044,
178 0.000915283345202241360843392549948,
179 0.005399280219300471367738743391053,
180 0.010947679601118931134327826856808,
181 0.016298731696787335262665703223280,
182 0.021081568889203835112433060188190,
183 0.025370969769253827243467999831710,
184 0.029189697756475752501446154084920,
185 0.032373202467202789685788194889595,
186 0.034783098950365142750781997949596,
187 0.036412220731351787562801163687577,
188 0.037253875503047708539592001191226
192 static const Real w87b[23] = {
193 0.000274145563762072350016527092881,
194 0.001807124155057942948341311753254,
195 0.004096869282759164864458070683480,
196 0.006758290051847378699816577897424,
197 0.009549957672201646536053581325377,
198 0.012329447652244853694626639963780,
199 0.015010447346388952376697286041943,
200 0.017548967986243191099665352925900,
201 0.019938037786440888202278192730714,
202 0.022194935961012286796332102959499,
203 0.024339147126000805470360647041454,
204 0.026374505414839207241503786552615,
205 0.028286910788771200659968002987960,
206 0.030052581128092695322521110347341,
207 0.031646751371439929404586051078883,
208 0.033050413419978503290785944862689,
209 0.034255099704226061787082821046821,
210 0.035262412660156681033782717998428,
211 0.036076989622888701185500318003895,
212 0.036698604498456094498018047441094,
213 0.037120549269832576114119958413599,
214 0.037334228751935040321235449094698,
215 0.037361073762679023410321241766599
229 Real relativeAccuracy)
230 :
Integrator(absoluteAccuracy, maxEvaluations),
231 relativeAccuracy_(relativeAccuracy) {}
239 Real fv1[5], fv2[5], fv3[5], fv4[5];
241 Real res10, res21, res43, res87;
249 const Real halfLength = 0.5 * (
b - a);
250 const Real center = 0.5 * (
b + a);
251 const Real fCenter =
f(center);
256 res21 = w21b[5] * fCenter;
257 resAbs = w21b[5] * std::fabs(fCenter);
259 for (k = 0; k < 5; k++) {
260 Real abscissa = halfLength * x1[k];
261 Real fval1 =
f(center + abscissa);
262 Real fval2 =
f(center - abscissa);
263 Real fval = fval1 + fval2;
264 res10 += w10[k] * fval;
265 res21 += w21a[k] * fval;
266 resAbs += w21a[k] * (std::fabs(fval1) + std::fabs(fval2));
272 for (k = 0; k < 5; k++) {
273 Real abscissa = halfLength * x2[k];
274 Real fval1 =
f(center + abscissa);
275 Real fval2 =
f(center - abscissa);
276 Real fval = fval1 + fval2;
277 res21 += w21b[k] * fval;
278 resAbs += w21b[k] * (std::fabs(fval1) + std::fabs(fval2));
279 savfun[k + 5] = fval;
284 result = res21 * halfLength;
285 resAbs *= halfLength ;
286 Real mean = 0.5 * res21;
287 resasc = w21b[5] * std::fabs(fCenter - mean);
289 for (k = 0; k < 5; k++)
290 resasc += (w21a[k] * (std::fabs(fv1[k] - mean)
291 + std::fabs(fv2[k] - mean))
292 + w21b[k] * (std::fabs(fv3[k] - mean)
293 + std::fabs(fv4[k] - mean)));
295 err = rescaleError ((res21 - res10) * halfLength, resAbs, resasc) ;
296 resasc *= halfLength ;
307 res43 = w43b[11] * fCenter;
309 for (k = 0; k < 10; k++)
310 res43 += savfun[k] * w43a[k];
312 for (k = 0; k < 11; k++){
313 Real abscissa = halfLength * x3[k];
314 Real fval = (
f(center + abscissa)
315 +
f(center - abscissa));
316 res43 += fval * w43b[k];
317 savfun[k + 10] = fval;
322 result = res43 * halfLength;
323 err = rescaleError ((res43 - res21) * halfLength, resAbs, resasc);
333 res87 = w87b[22] * fCenter;
335 for (k = 0; k < 21; k++)
336 res87 += savfun[k] * w87a[k];
338 for (k = 0; k < 22; k++){
339 Real abscissa = halfLength * x4[k];
340 res87 += w87b[k] * (
f(center + abscissa)
341 +
f(center - abscissa));
345 result = res87 * halfLength ;
346 err = rescaleError ((res87 - res43) * halfLength, resAbs, resasc);
362 static const Real g7w[] = { 0.417959183673469,
367 static const Real k15w[] = { 0.209482141084728,
377 static const Real k15t[] = { 0.000000000000000,
390 Real tolerance)
const {
392 Real halflength = (
b - a) / 2;
393 Real center = (a +
b) / 2;
405 for (j = 1, j2 = 2; j < 4; j++, j2 += 2) {
406 t = halflength * k15t[j2];
407 fsum =
f(center -
t) +
f(center +
t);
409 k15 += fsum * k15w[j2];
413 for (j2 = 1; j2 < 8; j2 += 2) {
414 t = halflength * k15t[j2];
415 fsum =
f(center -
t) +
f(center +
t);
416 k15 += fsum * k15w[j2];
420 g7 = halflength * g7;
421 k15 = halflength * k15;
429 if (std::fabs(k15 - g7) < tolerance) {
434 "maximum number of function evaluations "
444 :
Integrator(absoluteAccuracy, maxEvaluations) {
447 ") not allowed. It must be >= 15");
Real integrate(const ext::function< Real(Real)> &f, Real a, Real b) const override
GaussKronrodAdaptive(Real tolerance, Size maxFunctionEvaluations=Null< Size >())
Real integrateRecursively(const ext::function< Real(Real)> &f, Real a, Real b, Real tolerance) const
Real integrate(const ext::function< Real(Real)> &f, Real a, Real b) const override
void setRelativeAccuracy(Real)
Real relativeAccuracy() const
GaussKronrodNonAdaptive(Real absoluteAccuracy, Size maxEvaluations, Real relativeAccuracy)
Real absoluteAccuracy() const
Size numberOfEvaluations() const
Size maxEvaluations() const
void increaseNumberOfEvaluations(Size increase) const
void setAbsoluteError(Real error) const
void setNumberOfEvaluations(Size evaluations) const
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
ext::function< Real(Real)> b
#define QL_MIN_POSITIVE_REAL
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
Integral of a 1-dimensional function using the Gauss-Kronrod method.