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.qk21; 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 qk21(Real, Func)(Func f, Real a, Real b, out Real result, 22 out Real abserr, out Real resabs, out Real resasc) 23 { 24 //***begin prologue dqk21 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a1a2 28 //***keywords 21-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 driver 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 21-point 57 // kronrod rule (resk) obtained by optimal addition 58 // of abscissae to the 10-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 dqk21 74 // 75 Real absc,centr,dhlgth, 76 epmach,fc,fsum,fval1,fval2,hlgth, 77 resg,resk,reskh,uflow; 78 Real[10] fv1_, fv2_; 79 int j,jtw,jtwm1; 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 21-point kronrod rule 86 // xgk(2), xgk(4), ... abscissae of the 10-point 87 // gauss rule 88 // xgk(1), xgk(3), ... abscissae which are optimally 89 // added to the 10-point gauss rule 90 // 91 // wgk - weights of the 21-point kronrod rule 92 // 93 // wg - weights of the 10-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[5] wg_ = [ 101 0.0666713443_0868813759_3568809893_332, 102 0.1494513491_5058059314_5776339657_697, 103 0.2190863625_1598204399_5534934228_163, 104 0.2692667193_0999635509_1226921569_469, 105 0.2955242247_1475287017_3892994651_338]; 106 // 107 static immutable Real[11] xgk_ = [ 108 0.9956571630_2580808073_5527280689_003, 109 0.9739065285_1717172007_7964012084_452, 110 0.9301574913_5570822600_1207180059_508, 111 0.8650633666_8898451073_2096688423_493, 112 0.7808177265_8641689706_3717578345_042, 113 0.6794095682_9902440623_4327365114_874, 114 0.5627571346_6860468333_9000099272_694, 115 0.4333953941_2924719079_9265943165_784, 116 0.2943928627_0146019813_1126603103_866, 117 0.1488743389_8163121088_4826001129_720, 118 0.0000000000_0000000000_0000000000_000]; 119 // 120 static immutable Real[11] wgk_ = [ 121 0.0116946388_6737187427_8064396062_192, 122 0.0325581623_0796472747_8818972459_390, 123 0.0547558965_7435199603_1381300244_580, 124 0.0750396748_1091995276_7043140916_190, 125 0.0931254545_8369760553_5065465083_366, 126 0.1093871588_0229764189_9210590325_805, 127 0.1234919762_6206585107_7958109831_074, 128 0.1347092173_1147332592_8054001771_707, 129 0.1427759385_7706008079_7094273138_717, 130 0.1477391049_0133849137_4841515972_068, 131 0.1494455540_0291690566_4936468389_821]; 132 // 133 auto fv1 = dimension(fv1_.ptr, 10); 134 auto fv2 = dimension(fv2_.ptr, 10); 135 auto wg = dimension(wg_.ptr, 5); 136 auto wgk = dimension(wgk_.ptr, 11); 137 auto xgk = dimension(xgk_.ptr, 11); 138 // 139 // 140 // list of major variables 141 // ----------------------- 142 // 143 // centr - mid point of the interval 144 // hlgth - half-length of the interval 145 // absc - abscissa 146 // fval* - function value 147 // resg - result of the 10-point gauss formula 148 // resk - result of the 21-point kronrod formula 149 // reskh - approximation to the mean value of f over (a,b), 150 // i.e. to i/(b-a) 151 // 152 // 153 // machine dependent constants 154 // --------------------------- 155 // 156 // epmach is the largest relative spacing. 157 // uflow is the smallest positive magnitude. 158 // 159 //***first executable statement dqk21 160 epmach = Real.epsilon; 161 uflow = Real.min_normal; 162 // 163 centr = 0.5*(a+b); 164 hlgth = 0.5*(b-a); 165 dhlgth = fabs(hlgth); 166 // 167 // compute the 21-point kronrod approximation to 168 // the integral, and estimate the absolute error. 169 // 170 resg = 0.0; 171 fc = f(centr); 172 resk = wgk[11]*fc; 173 resabs = fabs(resk); 174 for (j=1; j<=5; j++) { //do 10 j = 1,5 175 jtw = 2*j; 176 absc = hlgth*xgk[jtw]; 177 fval1 = f(centr-absc); 178 fval2 = f(centr+absc); 179 fv1[jtw] = fval1; 180 fv2[jtw] = fval2; 181 fsum = fval1+fval2; 182 resg = resg+wg[j]*fsum; 183 resk = resk+wgk[jtw]*fsum; 184 resabs = resabs+wgk[jtw]*(fabs(fval1)+fabs(fval2)); 185 l10: ;} 186 for (j=1; j<=5; j++) { //do 15 j = 1,5 187 jtwm1 = 2*j-1; 188 absc = hlgth*xgk[jtwm1]; 189 fval1 = f(centr-absc); 190 fval2 = f(centr+absc); 191 fv1[jtwm1] = fval1; 192 fv2[jtwm1] = fval2; 193 fsum = fval1+fval2; 194 resk = resk+wgk[jtwm1]*fsum; 195 resabs = resabs+wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); 196 l15: ;} 197 reskh = resk*0.5; 198 resasc = wgk[11]*fabs(fc-reskh); 199 for (j=1; j<=10; j++) { //do 20 j=1,10 200 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 201 l20: ;} 202 result = resk*hlgth; 203 resabs = resabs*dhlgth; 204 resasc = resasc*dhlgth; 205 abserr = fabs((resk-resg)*hlgth); 206 if(resasc != 0.0 && abserr != 0.0) 207 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 208 if(resabs > uflow/(0.5e2*epmach)) abserr = max 209 ((epmach*0.5e2)*resabs,abserr); 210 return; 211 } 212 213 version(unittest) import scid.core.testing; 214 unittest 215 { 216 alias qk21!(float, float delegate(float)) fqk21; 217 alias qk21!(double, double delegate(double)) dqk21; 218 alias qk21!(double, double function(double)) dfqk21; 219 alias qk21!(real, real delegate(real)) rqk21; 220 221 double f(double x) { return x^^3; } 222 double result, abserr, resabs, resasc; 223 qk21(&f, 0.0, 1.0, result, abserr, resabs, resasc); 224 assert (isAccurate(result, abserr, 0.25, 1e-6)); 225 }