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_qk. 4 5 /** Authors: Lars Tandle Kyllingstad 6 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 7 License: Boost License 1.0 8 */ 9 module scid.ports.quadpack.qk41; 10 11 12 import std.algorithm: max, min; 13 import std.math; 14 15 import scid.core.fortran; 16 17 18 19 20 /// 21 void qk41(Real, Func)(Func f, Real a, Real b, out Real result, out Real abserr, 22 out Real resabs, out Real resasc) 23 { 24 //***begin prologue dqk41 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a1a2 28 //***keywords 41-point gauss-kronrod rules 29 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 30 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 31 //***purpose to compute i = integral of f over (a,b), with error 32 // estimate 33 // j = integral of abs(f) over (a,b) 34 //***description 35 // 36 // integration rules 37 // standard fortran subroutine 38 // double precision version 39 // 40 // parameters 41 // on entry 42 // f - double precision 43 // function subprogram defining the integrand 44 // function f(x). the actual name for f needs to be 45 // declared e x t e r n a l in the calling program. 46 // 47 // a - double precision 48 // lower limit of integration 49 // 50 // b - double precision 51 // upper limit of integration 52 // 53 // on return 54 // result - double precision 55 // approximation to the integral i 56 // result is computed by applying the 41-point 57 // gauss-kronrod rule (resk) obtained by optimal 58 // addition of abscissae to the 20-point gauss 59 // rule (resg). 60 // 61 // abserr - double precision 62 // estimate of the modulus of the absolute error, 63 // which should not exceed abs(i-result) 64 // 65 // resabs - double precision 66 // approximation to the integral j 67 // 68 // resasc - double precision 69 // approximation to the integal of abs(f-i/(b-a)) 70 // over (a,b) 71 // 72 //***references (none) 73 //***routines called d1mach 74 //***end prologue dqk41 75 // 76 Real absc,centr,dhlgth, 77 epmach,fc,fsum,fval1,fval2,hlgth, 78 resg,resk,reskh,uflow; 79 Real[20] fv1_, fv2_; 80 int j,jtw,jtwm1; 81 // 82 // the abscissae and weights are given for the interval (-1,1). 83 // because of symmetry only the positive abscissae and their 84 // corresponding weights are given. 85 // 86 // xgk - abscissae of the 41-point gauss-kronrod rule 87 // xgk(2), xgk(4), ... abscissae of the 20-point 88 // gauss rule 89 // xgk(1), xgk(3), ... abscissae which are optimally 90 // added to the 20-point gauss rule 91 // 92 // wgk - weights of the 41-point gauss-kronrod rule 93 // 94 // wg - weights of the 20-point gauss rule 95 // 96 // 97 // gauss quadrature weights and kronron quadrature abscissae and weights 98 // as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 99 // bell labs, nov. 1981. 100 // 101 static immutable Real[10] wg_ = [ 102 0.0176140071_3915211831_1861962351_853, 103 0.0406014298_0038694133_1039952274_932, 104 0.0626720483_3410906356_9506535187_042, 105 0.0832767415_7670474872_4758143222_046, 106 0.1019301198_1724043503_6750135480_350, 107 0.1181945319_6151841731_2377377711_382, 108 0.1316886384_4917662689_8494499748_163, 109 0.1420961093_1838205132_9298325067_165, 110 0.1491729864_7260374678_7828737001_969, 111 0.1527533871_3072585069_8084331955_098]; 112 // 113 static immutable Real[21] xgk_ = [ 114 0.9988590315_8827766383_8315576545_863, 115 0.9931285991_8509492478_6122388471_320, 116 0.9815078774_5025025919_3342994720_217, 117 0.9639719272_7791379126_7666131197_277, 118 0.9408226338_3175475351_9982722212_443, 119 0.9122344282_5132590586_7752441203_298, 120 0.8782768112_5228197607_7442995113_078, 121 0.8391169718_2221882339_4529061701_521, 122 0.7950414288_3755119835_0638833272_788, 123 0.7463319064_6015079261_4305070355_642, 124 0.6932376563_3475138480_5490711845_932, 125 0.6360536807_2651502545_2836696226_286, 126 0.5751404468_1971031534_2946036586_425, 127 0.5108670019_5082709800_4364050955_251, 128 0.4435931752_3872510319_9992213492_640, 129 0.3737060887_1541956067_2548177024_927, 130 0.3016278681_1491300432_0555356858_592, 131 0.2277858511_4164507808_0496195368_575, 132 0.1526054652_4092267550_5220241022_678, 133 0.0765265211_3349733375_4640409398_838, 134 0.0000000000_0000000000_0000000000_000]; 135 // 136 static immutable Real[21] wgk_ = [ 137 0.0030735837_1852053150_1218293246_031, 138 0.0086002698_5564294219_8661787950_102, 139 0.0146261692_5697125298_3787960308_868, 140 0.0203883734_6126652359_8010231432_755, 141 0.0258821336_0495115883_4505067096_153, 142 0.0312873067_7703279895_8543119323_801, 143 0.0366001697_5820079803_0557240707_211, 144 0.0416688733_2797368626_3788305936_895, 145 0.0464348218_6749767472_0231880926_108, 146 0.0509445739_2372869193_2707670050_345, 147 0.0551951053_4828599474_4832372419_777, 148 0.0591114008_8063957237_4967220648_594, 149 0.0626532375_5478116802_5870122174_255, 150 0.0658345971_3361842211_1563556969_398, 151 0.0686486729_2852161934_5623411885_368, 152 0.0710544235_5344406830_5790361723_210, 153 0.0730306903_3278666749_5189417658_913, 154 0.0745828754_0049918898_6581418362_488, 155 0.0757044976_8455667465_9542775376_617, 156 0.0763778676_7208073670_5502835038_061, 157 0.0766007119_1799965644_5049901530_102]; 158 // 159 auto fv1 = dimension(fv1_.ptr, 20); 160 auto fv2 = dimension(fv2_.ptr, 20); 161 auto xgk = dimension(xgk_.ptr, 21); 162 auto wgk = dimension(wgk_.ptr, 21); 163 auto wg = dimension(wg_.ptr, 10); 164 // 165 // 166 // list of major variables 167 // ----------------------- 168 // 169 // centr - mid point of the interval 170 // hlgth - half-length of the interval 171 // absc - abscissa 172 // fval* - function value 173 // resg - result of the 20-point gauss formula 174 // resk - result of the 41-point kronrod formula 175 // reskh - approximation to mean value of f over (a,b), i.e. 176 // to i/(b-a) 177 // 178 // machine dependent constants 179 // --------------------------- 180 // 181 // epmach is the largest relative spacing. 182 // uflow is the smallest positive magnitude. 183 // 184 //***first executable statement dqk41 185 epmach = Real.epsilon; 186 uflow = Real.min_normal; 187 // 188 centr = 0.5*(a+b); 189 hlgth = 0.5*(b-a); 190 dhlgth = fabs(hlgth); 191 // 192 // compute the 41-point gauss-kronrod approximation to 193 // the integral, and estimate the absolute error. 194 // 195 resg = 0.0; 196 fc = f(centr); 197 resk = wgk[21]*fc; 198 resabs = fabs(resk); 199 for (j=1; j<=10; j++) { //do 10 j=1,10 200 jtw = j*2; 201 absc = hlgth*xgk[jtw]; 202 fval1 = f(centr-absc); 203 fval2 = f(centr+absc); 204 fv1[jtw] = fval1; 205 fv2[jtw] = fval2; 206 fsum = fval1+fval2; 207 resg = resg+wg[j]*fsum; 208 resk = resk+wgk[jtw]*fsum; 209 resabs = resabs+wgk[jtw]*(fabs(fval1)+fabs(fval2)); 210 l10: ;} 211 for (j=1; j<=10; j++) { //do 15 j = 1,10 212 jtwm1 = j*2-1; 213 absc = hlgth*xgk[jtwm1]; 214 fval1 = f(centr-absc); 215 fval2 = f(centr+absc); 216 fv1[jtwm1] = fval1; 217 fv2[jtwm1] = fval2; 218 fsum = fval1+fval2; 219 resk = resk+wgk[jtwm1]*fsum; 220 resabs = resabs+wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); 221 l15: ;} 222 reskh = resk*0.5; 223 resasc = wgk[21]*fabs(fc-reskh); 224 for (j=1; j<=20; j++) { //do 20 j=1,20 225 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 226 l20: ;} 227 result = resk*hlgth; 228 resabs = resabs*dhlgth; 229 resasc = resasc*dhlgth; 230 abserr = fabs((resk-resg)*hlgth); 231 if(resasc != 0.0 && abserr != 0.) 232 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 233 if(resabs > uflow/(0.5e2*epmach)) abserr = max 234 ((epmach*0.5e2)*resabs,abserr); 235 return; 236 } 237 238 version(unittest) import scid.core.testing; 239 unittest 240 { 241 alias qk41!(float, float delegate(float)) fqk41; 242 alias qk41!(double, double delegate(double)) dqk41; 243 alias qk41!(double, double function(double)) dfqk41; 244 alias qk41!(real, real delegate(real)) rqk41; 245 246 double f(double x) { return x^^3; } 247 double result, abserr, resabs, resasc; 248 qk41(&f, 0.0, 1.0, result, abserr, resabs, resasc); 249 assert (isAccurate(result, abserr, 0.25, 1e-6)); 250 }