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.qk15; 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 qk15(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 dqk15 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a1a2 28 //***keywords 15-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 15-point 57 // kronrod rule (resk) obtained by optimal addition 58 // of abscissae to the7-point gauss rule(resg). 59 // 60 // abserr - double precision 61 // estimate of the modulus of the absolute error, 62 // which should not exceed abs(i-result) 63 // 64 // resabs - double precision 65 // approximation to the integral j 66 // 67 // resasc - double precision 68 // approximation to the integral of abs(f-i/(b-a)) 69 // over (a,b) 70 // 71 //***references (none) 72 //***routines called d1mach 73 //***end prologue dqk15 74 // 75 Real absc,centr,dhlgth, 76 epmach,fc,fsum,fval1,fval2,hlgth, 77 resg,resk,reskh,uflow; 78 Real[7] fv1_, fv2_; 79 int j=1,jtw=1,jtwm1=1; 80 // 81 // the abscissae and weights are given for the interval (-1,1). 82 // because of symmetry only the positive abscissae and their 83 // corresponding weights are given. 84 // 85 // xgk - abscissae of the 15-point kronrod rule 86 // xgk(2), xgk(4), ... abscissae of the 7-point 87 // gauss rule 88 // xgk(1), xgk(3), ... abscissae which are optimally 89 // added to the 7-point gauss rule 90 // 91 // wgk - weights of the 15-point kronrod rule 92 // 93 // wg - weights of the 7-point gauss rule 94 // 95 // 96 // gauss quadrature weights and kronron quadrature abscissae and weights 97 // as evaluated with 80 decimal digit arithmetic by l. w. fullerton, 98 // bell labs, nov. 1981. 99 // 100 static immutable Real[4] wg_ = [ 101 0.1294849661_6886969327_0611432679_082, 102 0.2797053914_8927666790_1467771423_780, 103 0.3818300505_0511894495_0369775488_975, 104 0.4179591836_7346938775_5102040816_327]; 105 // 106 static immutable Real[8] xgk_ = [ 107 0.9914553711_2081263920_6854697526_329, 108 0.9491079123_4275852452_6189684047_851, 109 0.8648644233_5976907278_9712788640_926, 110 0.7415311855_9939443986_3864773280_788, 111 0.5860872354_6769113029_4144838258_730, 112 0.4058451513_7739716690_6606412076_961, 113 0.2077849550_0789846760_0689403773_245, 114 0.0000000000_0000000000_0000000000_000]; 115 // 116 static immutable Real[8] wgk_ = [ 117 0.0229353220_1052922496_3732008058_970, 118 0.0630920926_2997855329_0700663189_204, 119 0.1047900103_2225018383_9876322541_518, 120 0.1406532597_1552591874_5189590510_238, 121 0.1690047266_3926790282_6583426598_550, 122 0.1903505780_6478540991_3256402421_014, 123 0.2044329400_7529889241_4161999234_649, 124 0.2094821410_8472782801_2999174891_714]; 125 // 126 auto fv1 = dimension(fv1_.ptr, 7); 127 auto fv2 = dimension(fv2_.ptr, 7); 128 auto wg = dimension(wg_.ptr, 4); 129 auto wgk = dimension(wgk_.ptr, 8); 130 auto xgk = dimension(xgk_.ptr, 8); 131 // 132 // 133 // list of major variables 134 // ----------------------- 135 // 136 // centr - mid point of the interval 137 // hlgth - half-length of the interval 138 // absc - abscissa 139 // fval* - function value 140 // resg - result of the 7-point gauss formula 141 // resk - result of the 15-point kronrod formula 142 // reskh - approximation to the mean value of f over (a,b), 143 // i.e. to i/(b-a) 144 // 145 // machine dependent constants 146 // --------------------------- 147 // 148 // epmach is the largest relative spacing. 149 // uflow is the smallest positive magnitude. 150 // 151 //***first executable statement dqk15 152 epmach = Real.epsilon; 153 uflow = Real.min_normal; 154 // 155 centr = 0.5*(a+b); 156 hlgth = 0.5*(b-a); 157 dhlgth = fabs(hlgth); 158 // 159 // compute the 15-point kronrod approximation to 160 // the integral, and estimate the absolute error. 161 // 162 fc = f(centr); 163 resg = fc*wg[4]; 164 resk = fc*wgk[8]; 165 resabs = fabs(resk); 166 for (j=1; j<=3; j++) { //do 10 j=1,3 167 jtw = j*2; 168 absc = hlgth*xgk[jtw]; 169 fval1 = f(centr-absc); 170 fval2 = f(centr+absc); 171 fv1[jtw] = fval1; 172 fv2[jtw] = fval2; 173 fsum = fval1+fval2; 174 resg = resg+wg[j]*fsum; 175 resk = resk+wgk[jtw]*fsum; 176 resabs = resabs+wgk[jtw]*(fabs(fval1)+fabs(fval2)); 177 l10: ;} 178 for (j=1; j<=4; j++) { //do 15 j = 1,4 179 jtwm1 = j*2-1; 180 absc = hlgth*xgk[jtwm1]; 181 fval1 = f(centr-absc); 182 fval2 = f(centr+absc); 183 fv1[jtwm1] = fval1; 184 fv2[jtwm1] = fval2; 185 fsum = fval1+fval2; 186 resk = resk+wgk[jtwm1]*fsum; 187 resabs = resabs+wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); 188 l15: ;} 189 reskh = resk*0.5; 190 resasc = wgk[8]*fabs(fc-reskh); 191 for (j=1; j<=7; j++) { //do 20 j=1,7 192 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 193 l20: ;} 194 result = resk*hlgth; 195 resabs = resabs*dhlgth; 196 resasc = resasc*dhlgth; 197 abserr = fabs((resk-resg)*hlgth); 198 if(resasc != 0.0 && abserr != 0.0) 199 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 200 if(resabs > uflow/(0.5e2*epmach)) abserr = max 201 ((epmach*0.5e2)*resabs,abserr); 202 return; 203 } 204 205 version(unittest) import scid.core.testing; 206 unittest 207 { 208 alias qk15!(float, float delegate(float)) fqk15; 209 alias qk15!(double, double delegate(double)) dqk15; 210 alias qk15!(double, double function(double)) dfqk15; 211 alias qk15!(real, real delegate(real)) rqk15; 212 213 double f(double x) { return x^^3; } 214 double result, abserr, resabs, resasc; 215 qk15(&f, 0.0, 1.0, result, abserr, resabs, resasc); 216 assert (isAccurate(result, abserr, 0.25, 1e-6)); 217 }