1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/quadpack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.quadpack.qk15i; 9 10 11 import std.algorithm: max, min; 12 import std.math; 13 14 import scid.core.fortran; 15 import scid.core.meta; 16 17 18 19 20 /// 21 void qk15i(Real, Func)(Func f, Real boun, int inf, Real a, Real b, 22 out Real result, out Real abserr, out Real resabs, out Real resasc) 23 { 24 //***begin prologue dqk15i 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a3a2,h2a4a2 28 //***keywords 15-point transformed 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 the original (infinite integration range is mapped 32 // onto the interval (0,1) and (a,b) is a part of (0,1). 33 // it is the purpose to compute 34 // i = integral of transformed integrand over (a,b), 35 // j = integral of abs(transformed integrand) over (a,b). 36 //***description 37 // 38 // integration rule 39 // standard fortran subroutine 40 // double precision version 41 // 42 // parameters 43 // on entry 44 // f - double precision 45 // fuction subprogram defining the integrand 46 // function f(x). the actual name for f needs to be 47 // declared e x t e r n a l in the calling program. 48 // 49 // boun - double precision 50 // finite bound of original integration 51 // range (set to zero if inf = +2) 52 // 53 // inf - integer 54 // if inf = -1, the original interval is 55 // (-infinity,bound), 56 // if inf = +1, the original interval is 57 // (bound,+infinity), 58 // if inf = +2, the original interval is 59 // (-infinity,+infinity) and 60 // the integral is computed as the sum of two 61 // integrals, one over (-infinity,0) and one over 62 // (0,+infinity). 63 // 64 // a - double precision 65 // lower limit for integration over subrange 66 // of (0,1) 67 // 68 // b - double precision 69 // upper limit for integration over subrange 70 // of (0,1) 71 // 72 // on return 73 // result - double precision 74 // approximation to the integral i 75 // result is computed by applying the 15-point 76 // kronrod rule(resk) obtained by optimal addition 77 // of abscissae to the 7-point gauss rule(resg). 78 // 79 // abserr - double precision 80 // estimate of the modulus of the absolute error, 81 // which should equal or exceed abs(i-result) 82 // 83 // resabs - double precision 84 // approximation to the integral j 85 // 86 // resasc - double precision 87 // approximation to the integral of 88 // abs((transformed integrand)-i/(b-a)) over (a,b) 89 // 90 //***references (none) 91 //***routines called d1mach 92 //***end prologue dqk15i 93 // 94 Real absc,absc1,absc2,centr,dabs,dinf, 95 epmach,fc,fsum,fval1,fval2,hlgth, 96 resg,resk,reskh,tabsc1,tabsc2,uflow; 97 Real[7] fv1_, fv2_; 98 int j; 99 // 100 // the abscissae and weights are supplied for the interval 101 // (-1,1). because of symmetry only the positive abscissae and 102 // their corresponding weights are given. 103 // 104 // xgk - abscissae of the 15-point kronrod rule 105 // xgk(2), xgk(4), ... abscissae of the 7-point 106 // gauss rule 107 // xgk(1), xgk(3), ... abscissae which are optimally 108 // added to the 7-point gauss rule 109 // 110 // wgk - weights of the 15-point kronrod rule 111 // 112 // wg - weights of the 7-point gauss rule, corresponding 113 // to the abscissae xgk(2), xgk(4), ... 114 // wg(1), wg(3), ... are set to zero. 115 // 116 static immutable Real[8] wg_ = [ 117 0.0, 118 0.1294849661_6886969327_0611432679_082, 119 0.0, 120 0.2797053914_8927666790_1467771423_780, 121 0.0, 122 0.3818300505_0511894495_0369775488_975, 123 0.0, 124 0.4179591836_7346938775_5102040816_327 125 ]; 126 // 127 static immutable Real[8] xgk_ = [ 128 0.9914553711_2081263920_6854697526_329, 129 0.9491079123_4275852452_6189684047_851, 130 0.8648644233_5976907278_9712788640_926, 131 0.7415311855_9939443986_3864773280_788, 132 0.5860872354_6769113029_4144838258_730, 133 0.4058451513_7739716690_6606412076_961, 134 0.2077849550_0789846760_0689403773_245, 135 0.0000000000_0000000000_0000000000_000 136 ]; 137 // 138 static immutable Real[8] wgk_ = [ 139 0.0229353220_1052922496_3732008058_970, 140 0.0630920926_2997855329_0700663189_204, 141 0.1047900103_2225018383_9876322541_518, 142 0.1406532597_1552591874_5189590510_238, 143 0.1690047266_3926790282_6583426598_550, 144 0.1903505780_6478540991_3256402421_014, 145 0.2044329400_7529889241_4161999234_649, 146 0.2094821410_8472782801_2999174891_714 147 ]; 148 // 149 auto fv1 = dimension(fv1_.ptr, 7); 150 auto fv2 = dimension(fv2_.ptr, 7); 151 auto xgk = dimension(xgk_.ptr, 8); 152 auto wgk = dimension(wgk_.ptr, 8); 153 auto wg = dimension(wg_.ptr, 8); 154 // 155 // 156 // list of major variables 157 // ----------------------- 158 // 159 // centr - mid point of the interval 160 // hlgth - half-length of the interval 161 // absc* - abscissa 162 // tabsc* - transformed abscissa 163 // fval* - function value 164 // resg - result of the 7-point gauss formula 165 // resk - result of the 15-point kronrod formula 166 // reskh - approximation to the mean value of the transformed 167 // integrand over (a,b), i.e. to i/(b-a) 168 // 169 // machine dependent constants 170 // --------------------------- 171 // 172 // epmach is the largest relative spacing. 173 // uflow is the smallest positive magnitude. 174 // 175 //***first executable statement dqk15i 176 epmach = Real.epsilon; 177 uflow = Real.min_normal; 178 dinf = min(1,inf); 179 // 180 centr = 0.5*(a+b); 181 hlgth = 0.5*(b-a); 182 tabsc1 = boun+dinf*(0.1e1-centr)/centr; 183 fval1 = f(tabsc1); 184 if(inf == 2) fval1 = fval1+f(-tabsc1); 185 fc = (fval1/centr)/centr; 186 // 187 // compute the 15-point kronrod approximation to 188 // the integral, and estimate the error. 189 // 190 resg = wg[8]*fc; 191 resk = wgk[8]*fc; 192 resabs = fabs(resk); 193 for (j=1; j<=7; j++) { // end: 10 194 absc = hlgth*xgk[j]; 195 absc1 = centr-absc; 196 absc2 = centr+absc; 197 tabsc1 = boun+dinf*(0.1e1-absc1)/absc1; 198 tabsc2 = boun+dinf*(0.1e1-absc2)/absc2; 199 fval1 = f(tabsc1); 200 fval2 = f(tabsc2); 201 if(inf == 2) fval1 = fval1+f(-tabsc1); 202 if(inf == 2) fval2 = fval2+f(-tabsc2); 203 fval1 = (fval1/absc1)/absc1; 204 fval2 = (fval2/absc2)/absc2; 205 fv1[j] = fval1; 206 fv2[j] = fval2; 207 fsum = fval1+fval2; 208 resg = resg+wg[j]*fsum; 209 resk = resk+wgk[j]*fsum; 210 resabs = resabs+wgk[j]*(fabs(fval1)+fabs(fval2)); 211 l10: ;} 212 reskh = resk*0.5; 213 resasc = wgk[8]*fabs(fc-reskh); 214 for (j=1; j<=7; j++) { // end: 20 215 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 216 l20: ;} 217 result = resk*hlgth; 218 resasc = resasc*hlgth; 219 resabs = resabs*hlgth; 220 abserr = fabs((resk-resg)*hlgth); 221 if(resasc != 0.0 && abserr != 0.0) abserr = resasc* 222 min(0.1e1, ((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 223 if(resabs > uflow/(0.5e2*epmach)) abserr = max 224 ((epmach*0.5e2)*resabs,abserr); 225 return; 226 } 227 228 229 unittest 230 { 231 alias qk15i!(float, float delegate(float)) fqk15i; 232 alias qk15i!(double, double delegate(double)) dqk15i; 233 alias qk15i!(double, double function(double)) dfqk15i; 234 alias qk15i!(real, real delegate(real)) rqk15i; 235 } 236