1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/quadpack. 3 // An idiomatic D port can be found in scid.internal.calculus.integrate_qng. 4 5 /** This module is deprecated in favour of scid.internal.calculus.integrate_qng. 6 7 Authors: Lars Tandle Kyllingstad 8 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 9 License: Boost License 1.0 10 */ 11 module scid.ports.quadpack.qng; 12 13 14 deprecated: 15 16 17 import std.algorithm: max, min; 18 import std.conv; 19 import std.math; 20 import std.traits; 21 22 import scid.core.fortran; 23 import scid.types; 24 25 version(unittest) 26 { 27 import scid.core.testing; 28 } 29 30 31 32 33 /// 34 void qng(Real, Func)(Func f, Real a, Real b, Real epsabs, Real epsrel, 35 out Real result, out Real abserr, out int neval, out int ier) 36 { 37 //***begin prologue dqng 38 //***date written 800101 (yymmdd) 39 //***revision date 810101 (yymmdd) 40 //***category no. h2a1a1 41 //***keywords automatic integrator, smooth integrand, 42 // non-adaptive, gauss-kronrod(patterson) 43 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 44 // de doncker,elise,appl math & progr. div. - k.u.leuven 45 // kahaner,david,nbs - modified (2/82) 46 //***purpose the routine calculates an approximation result to a 47 // given definite integral i = integral of f over (a,b), 48 // hopefully satisfying following claim for accuracy 49 // abs(i-result).le.max(epsabs,epsrel*abs(i)). 50 //***description 51 // 52 // non-adaptive integration 53 // standard fortran subroutine 54 // double precision version 55 // 56 // f - double precision 57 // function subprogram defining the integrand function 58 // f(x). the actual name for f needs to be declared 59 // e x t e r n a l in the driver program. 60 // 61 // a - double precision 62 // lower limit of integration 63 // 64 // b - double precision 65 // upper limit of integration 66 // 67 // epsabs - double precision 68 // absolute accuracy requested 69 // epsrel - double precision 70 // relative accuracy requested 71 // if epsabs.le.0 72 // and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 73 // the routine will end with ier = 6. 74 // 75 // on return 76 // result - double precision 77 // approximation to the integral i 78 // result is obtained by applying the 21-point 79 // gauss-kronrod rule (res21) obtained by optimal 80 // addition of abscissae to the 10-point gauss rule 81 // (res10), or by applying the 43-point rule (res43) 82 // obtained by optimal addition of abscissae to the 83 // 21-point gauss-kronrod rule, or by applying the 84 // 87-point rule (res87) obtained by optimal addition 85 // of abscissae to the 43-point rule. 86 // 87 // abserr - double precision 88 // estimate of the modulus of the absolute error, 89 // which should equal or exceed abs(i-result) 90 // 91 // neval - integer 92 // number of integrand evaluations 93 // 94 // ier - ier = 0 normal and reliable termination of the 95 // routine. it is assumed that the requested 96 // accuracy has been achieved. 97 // ier.gt.0 abnormal termination of the routine. it is 98 // assumed that the requested accuracy has 99 // not been achieved. 100 // error messages 101 // ier = 1 the maximum number of steps has been 102 // executed. the integral is probably too 103 // difficult to be calculated by dqng. 104 // = 6 the input is invalid, because 105 // epsabs.le.0 and 106 // epsrel.lt.max(50*rel.mach.acc.,0.5d-28). 107 // result, abserr and neval are set to zero. 108 // 109 //***references (none) 110 //***routines called d1mach,xerror 111 //***end prologue dqng 112 // 113 Real absc,centr,dhlgth, 114 epmach,fcentr,fval,fval1,fval2, 115 hlgth,res10,res21,res43,res87,resabs,resasc, 116 reskh,uflow; 117 Real[5] fv1_, fv2_, fv3_, fv4_; 118 Real[21] savfun_; 119 int ipx=1,k=1,l=1; 120 // 121 // the following data statements contain the 122 // abscissae and weights of the integration rules used. 123 // 124 // x1 abscissae common to the 10-, 21-, 43- and 87- 125 // point rule 126 // x2 abscissae common to the 21-, 43- and 87-point rule 127 // x3 abscissae common to the 43- and 87-point rule 128 // x4 abscissae of the 87-point rule 129 // w10 weights of the 10-point formula 130 // w21a weights of the 21-point formula for abscissae x1 131 // w21b weights of the 21-point formula for abscissae x2 132 // w43a weights of the 43-point formula for abscissae x1, x3 133 // w43b weights of the 43-point formula for abscissae x3 134 // w87a weights of the 87-point formula for abscissae x1, 135 // x2, x3 136 // w87b weights of the 87-point formula for abscissae x4 137 // 138 // 139 // gauss-kronrod-patterson quadrature coefficients for use in 140 // quadpack routine qng. these coefficients were calculated with 141 // 101 decimal digit arithmetic by l. w. fullerton, bell labs, nov 1981. 142 // 143 static immutable Real[5] x1_ = [ 144 0.9739065285_1717172007_7964012084_452, 145 0.8650633666_8898451073_2096688423_493, 146 0.6794095682_9902440623_4327365114_874, 147 0.4333953941_2924719079_9265943165_784, 148 0.1488743389_8163121088_4826001129_720]; 149 static immutable Real[5] w10_ = [ 150 0.0666713443_0868813759_3568809893_332, 151 0.1494513491_5058059314_5776339657_697, 152 0.2190863625_1598204399_5534934228_163, 153 0.2692667193_0999635509_1226921569_469, 154 0.2955242247_1475287017_3892994651_338]; 155 // 156 static immutable Real[5] x2_ = [ 157 0.9956571630_2580808073_5527280689_003, 158 0.9301574913_5570822600_1207180059_508, 159 0.7808177265_8641689706_3717578345_042, 160 0.5627571346_6860468333_9000099272_694, 161 0.2943928627_0146019813_1126603103_866]; 162 static immutable Real[5] w21a_ = [ 163 0.0325581623_0796472747_8818972459_390, 164 0.0750396748_1091995276_7043140916_190, 165 0.1093871588_0229764189_9210590325_805, 166 0.1347092173_1147332592_8054001771_707, 167 0.1477391049_0133849137_4841515972_068]; 168 static immutable Real[6] w21b_ = [ 169 0.0116946388_6737187427_8064396062_192, 170 0.0547558965_7435199603_1381300244_580, 171 0.0931254545_8369760553_5065465083_366, 172 0.1234919762_6206585107_7958109831_074, 173 0.1427759385_7706008079_7094273138_717, 174 0.1494455540_0291690566_4936468389_821]; 175 // 176 static immutable Real[11] x3_ = [ 177 0.9993333609_0193208139_4099323919_911, 178 0.9874334029_0808886979_5961478381_209, 179 0.9548079348_1426629925_7919200290_473, 180 0.9001486957_4832829362_5099494069_092, 181 0.8251983149_8311415084_7066732588_520, 182 0.7321483889_8930498261_2354848755_461, 183 0.6228479705_3772523864_1159120344_323, 184 0.4994795740_7105649995_2214885499_755, 185 0.3649016613_4658076804_3989548502_644, 186 0.2222549197_7660129649_8260928066_212, 187 0.0746506174_6138332204_3914435796_506]; 188 static immutable Real[10] w43a_ = [ 189 0.0162967342_8966656492_4281974617_663, 190 0.0375228761_2086950146_1613795898_115, 191 0.0546949020_5825544214_7212685465_005, 192 0.0673554146_0947808607_5553166302_174, 193 0.0738701996_3239395343_2140695251_367, 194 0.0057685560_5976979618_4184327908_655, 195 0.0273718905_9324884208_1276069289_151, 196 0.0465608269_1042883074_3339154433_824, 197 0.0617449952_0144256449_6240336030_883, 198 0.0713872672_6869339776_8559114425_516]; 199 static immutable Real[12] w43b_ = [ 200 0.0018444776_4021241410_0389106552_965, 201 0.0107986895_8589165174_0465406741_293, 202 0.0218953638_6779542810_2523123075_149, 203 0.0325974639_7534568944_3882222526_137, 204 0.0421631379_3519181184_7627924327_955, 205 0.0507419396_0018457778_0189020092_084, 206 0.0583793955_4261924837_5475369330_206, 207 0.0647464049_5144588554_4689259517_511, 208 0.0695661979_1235648452_8633315038_405, 209 0.0728244414_7183320815_0939535192_842, 210 0.0745077510_1417511827_3571813842_889, 211 0.0747221475_1740300559_4425168280_423]; 212 // 213 static immutable Real[22] x4_ = [ 214 0.9999029772_6272923449_0529830591_582, 215 0.9979898959_8667874542_7496322365_960, 216 0.9921754978_6068722280_8523352251_425, 217 0.9813581635_7271277357_1916941623_894, 218 0.9650576238_5838461912_8284110607_926, 219 0.9431676131_3367059681_6416634507_426, 220 0.9158064146_8550720959_1826430720_050, 221 0.8832216577_7131650137_2117548744_163, 222 0.8457107484_6241566660_5902011504_855, 223 0.8035576580_3523098278_8739474980_964, 224 0.7570057306_8549555832_8942793432_020, 225 0.7062732097_8732181982_4094274740_840, 226 0.6515894665_0117792253_4422205016_736, 227 0.5932233740_5796108887_5273770349_144, 228 0.5314936059_7083193228_5268948562_671, 229 0.4667636230_4202284487_1966781659_270, 230 0.3994248478_5921880473_2101665817_923, 231 0.3298748771_0618828826_5053371824_597, 232 0.2585035592_0216155180_2280975429_025, 233 0.1856953965_6834665201_5917141167_606, 234 0.1118422131_7990746817_2398359241_362, 235 0.0373521233_9461987081_4998165437_704]; 236 static immutable Real[21] w87a_ = [ 237 0.0081483773_8414917290_0002878448_190, 238 0.0187614382_0156282224_3935059003_794, 239 0.0273474510_5005228616_1582829741_283, 240 0.0336777073_1163793004_6581056957_588, 241 0.0369350998_2042790761_4589586742_499, 242 0.0028848724_3021153050_1334156248_695, 243 0.0136859460_2271270188_8950035273_128, 244 0.0232804135_0288831112_3409291030_404, 245 0.0308724976_1171335867_5466394126_442, 246 0.0356936336_3941877071_9351355457_044, 247 0.0009152833_4520224136_0843392549_948, 248 0.0053992802_1930047136_7738743391_053, 249 0.0109476796_0111893113_4327826856_808, 250 0.0162987316_9678733526_2665703223_280, 251 0.0210815688_8920383511_2433060188_190, 252 0.0253709697_6925382724_3467999831_710, 253 0.0291896977_5647575250_1446154084_920, 254 0.0323732024_6720278968_5788194889_595, 255 0.0347830989_5036514275_0781997949_596, 256 0.0364122207_3135178756_2801163687_577, 257 0.0372538755_0304770853_9592001191_226]; 258 static immutable Real[23] w87b_ = [ 259 0.0002741455_6376207235_0016527092_881, 260 0.0018071241_5505794294_8341311753_254, 261 0.0040968692_8275916486_4458070683_480, 262 0.0067582900_5184737869_9816577897_424, 263 0.0095499576_7220164653_6053581325_377, 264 0.0123294476_5224485369_4626639963_780, 265 0.0150104473_4638895237_6697286041_943, 266 0.0175489679_8624319109_9665352925_900, 267 0.0199380377_8644088820_2278192730_714, 268 0.0221949359_6101228679_6332102959_499, 269 0.0243391471_2600080547_0360647041_454, 270 0.0263745054_1483920724_1503786552_615, 271 0.0282869107_8877120065_9968002987_960, 272 0.0300525811_2809269532_2521110347_341, 273 0.0316467513_7143992940_4586051078_883, 274 0.0330504134_1997850329_0785944862_689, 275 0.0342550997_0422606178_7082821046_821, 276 0.0352624126_6015668103_3782717998_428, 277 0.0360769896_2288870118_5500318003_895, 278 0.0366986044_9845609449_8018047441_094, 279 0.0371205492_6983257611_4119958413_599, 280 0.0373342287_5193504032_1235449094_698, 281 0.0373610737_6267902341_0321241766_599]; 282 // 283 auto fv1 = dimension(fv1_.ptr, 5); 284 auto fv2 = dimension(fv2_.ptr, 5); 285 auto fv3 = dimension(fv3_.ptr, 5); 286 auto fv4 = dimension(fv4_.ptr, 5); 287 auto x1 = dimension(x1_.ptr, 5); 288 auto x2 = dimension(x2_.ptr, 5); 289 auto x3 = dimension(x3_.ptr, 11); 290 auto x4 = dimension(x4_.ptr, 22); 291 auto w10 = dimension(w10_.ptr, 5); 292 auto w21a = dimension(w21a_.ptr, 5); 293 auto w21b = dimension(w21b_.ptr, 6); 294 auto w43a = dimension(w43a_.ptr, 10); 295 auto w43b = dimension(w43b_.ptr, 12); 296 auto w87a = dimension(w87a_.ptr, 21); 297 auto w87b = dimension(w87b_.ptr, 23); 298 auto savfun = dimension(savfun_.ptr, 21); 299 // 300 // list of major variables 301 // ----------------------- 302 // 303 // centr - mid point of the integration interval 304 // hlgth - half-length of the integration interval 305 // fcentr - function value at mid point 306 // absc - abscissa 307 // fval - function value 308 // savfun - array of function values which have already been 309 // computed 310 // res10 - 10-point gauss result 311 // res21 - 21-point kronrod result 312 // res43 - 43-point result 313 // res87 - 87-point result 314 // resabs - approximation to the integral of abs(f) 315 // resasc - approximation to the integral of abs(f-i/(b-a)) 316 // 317 // machine dependent constants 318 // --------------------------- 319 // 320 // epmach is the largest relative spacing. 321 // uflow is the smallest positive magnitude. 322 // 323 //***first executable statement dqng 324 epmach = Real.epsilon; 325 uflow = Real.min_normal; 326 // 327 // test on validity of parameters 328 // ------------------------------ 329 // 330 result = 0.0; 331 abserr = 0.0; 332 neval = 0; 333 ier = 6; 334 if(epsabs <= 0.0 && epsrel < max(0.5e2*epmach,0.5e-28)) 335 goto l80; 336 hlgth = 0.5*(b-a); 337 dhlgth = fabs(hlgth); 338 centr = 0.5*(b+a); 339 fcentr = f(centr); 340 neval = 21; 341 ier = 1; 342 // 343 // compute the integral using the 10- and 21-point formula. 344 // 345 for (l=1; l<=3; l++) { //do 70 l = 1,3 346 if (l == 1) goto l5; 347 else if (l == 2) goto l25; 348 else if (l == 3) goto l45; 349 //goto l(5,25,45),l 350 l5: res10 = 0.0; 351 res21 = w21b[6]*fcentr; 352 resabs = w21b[6]*fabs(fcentr); 353 for (k=1; k<=5; k++) { //do 10 k=1,5 354 absc = hlgth*x1[k]; 355 fval1 = f(centr+absc); 356 fval2 = f(centr-absc); 357 fval = fval1+fval2; 358 res10 = res10+w10[k]*fval; 359 res21 = res21+w21a[k]*fval; 360 resabs = resabs+w21a[k]*(fabs(fval1)+fabs(fval2)); 361 savfun[k] = fval; 362 fv1[k] = fval1; 363 fv2[k] = fval2; 364 l10: ;} 365 ipx = 5; 366 for (k=1; k<=5; k++) { //do 15 k=1,5 367 ipx = ipx+1; 368 absc = hlgth*x2[k]; 369 fval1 = f(centr+absc); 370 fval2 = f(centr-absc); 371 fval = fval1+fval2; 372 res21 = res21+w21b[k]*fval; 373 resabs = resabs+w21b[k]*(fabs(fval1)+fabs(fval2)); 374 savfun[ipx] = fval; 375 fv3[k] = fval1; 376 fv4[k] = fval2; 377 l15: ;} 378 // 379 // test for convergence. 380 // 381 result = res21*hlgth; 382 resabs = resabs*dhlgth; 383 reskh = 0.5*res21; 384 resasc = w21b[6]*fabs(fcentr-reskh); 385 for (k=1; k<=5; k++) { //do 20 k = 1,5 386 resasc = resasc+w21a[k]*(fabs(fv1[k]-reskh)+fabs(fv2[k]-reskh)) 387 +w21b[k]*(fabs(fv3[k]-reskh)+fabs(fv4[k]-reskh)); 388 l20: ;} 389 abserr = fabs((res21-res10)*hlgth); 390 resasc = resasc*dhlgth; 391 goto l65; 392 // 393 // compute the integral using the 43-point formula. 394 // 395 l25: res43 = w43b[12]*fcentr; 396 neval = 43; 397 for (k=1; k<=10; k++) { //do 30 k=1,10 398 res43 = res43+savfun[k]*w43a[k]; 399 l30: ;} 400 for (k=1; k<=11; k++) { //do 40 k=1,11 401 ipx = ipx+1; 402 absc = hlgth*x3[k]; 403 fval = f(absc+centr)+f(centr-absc); 404 res43 = res43+fval*w43b[k]; 405 savfun[ipx] = fval; 406 l40: ;} 407 // 408 // test for convergence. 409 // 410 result = res43*hlgth; 411 abserr = fabs((res43-res21)*hlgth); 412 goto l65; 413 // 414 // compute the integral using the 87-point formula. 415 // 416 l45: res87 = w87b[23]*fcentr; 417 neval = 87; 418 for (k=1; k<=21; k++) { //do 50 k=1,21 419 res87 = res87+savfun[k]*w87a[k]; 420 l50: ;} 421 for (k=1; k<=22; k++) { // do 60 k=1,22 422 absc = hlgth*x4[k]; 423 res87 = res87+w87b[k]*(f(absc+centr)+f(centr-absc)); 424 l60: ;} 425 result = res87*hlgth; 426 abserr = fabs((res87-res43)*hlgth); 427 l65: if(resasc != 0.0 && abserr != 0.0) 428 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 429 if (resabs > uflow/(0.5e2*epmach)) abserr = max 430 ((epmach*0.5e2)*resabs,abserr); 431 if (abserr <= max(epsabs,epsrel*fabs(result))) ier = 0; 432 // ***jump out of do-loop 433 if (ier == 0) goto l999; 434 l70: ;} 435 l80: throw new Exception("abnormal return from qng: "~to!string(ier)); 436 l999: return; 437 } 438 439 unittest 440 { 441 alias qng!(float, float delegate(float)) fqng; 442 alias qng!(double, double delegate(double)) dqng; 443 alias qng!(double, double function(double)) dfqng; 444 alias qng!(real, real delegate(real)) rqng; 445 446 float f(float x) { return x <= 0.0 ? 0.0 : sqrt(x)*log(x); } 447 448 enum : float 449 { 450 a = 0.0, 451 b = 1.0, 452 epsabs = 0.0, 453 epsrel = 0.001 454 } 455 float result, abserr; 456 int neval, ier; 457 458 qng(&f, a, b, epsabs, epsrel, result, abserr, neval, ier); 459 460 assert (isAccurate(result, abserr, -4.0f/9.0f, epsrel, epsabs)); 461 }