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.qk31; 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 qk31(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 dqk31 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a1a2 28 //***keywords 31-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 31-point 57 // gauss-kronrod rule (resk), obtained by optimal 58 // addition of abscissae to the 15-point gauss 59 // rule (resg). 60 // 61 // abserr - double precison 62 // estimate of the modulus of the modulus, 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 integral of abs(f-i/(b-a)) 70 // over (a,b) 71 // 72 //***references (none) 73 //***routines called d1mach 74 //***end prologue dqk31 75 Real absc,centr,dhlgth, 76 epmach,fc,fsum,fval1,fval2,hlgth, 77 resg,resk,reskh,uflow; 78 Real[15] 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 31-point kronrod rule 86 // xgk(2), xgk(4), ... abscissae of the 15-point 87 // gauss rule 88 // xgk(1), xgk(3), ... abscissae which are optimally 89 // added to the 15-point gauss rule 90 // 91 // wgk - weights of the 31-point kronrod rule 92 // 93 // wg - weights of the 15-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 immutable static Real[8] wg_ = [ 101 0.0307532419_9611726835_4628393577_204, 102 0.0703660474_8810812470_9267416450_667, 103 0.1071592204_6717193501_1869546685_869, 104 0.1395706779_2615431444_7804794511_028, 105 0.1662692058_1699393355_3200860481_209, 106 0.1861610000_1556221102_6800561866_423, 107 0.1984314853_2711157645_6118326443_839, 108 0.2025782419_2556127288_0620199967_519]; 109 // 110 immutable static Real[16] xgk_ = [ 111 0.9980022986_9339706028_5172840152_271, 112 0.9879925180_2048542848_9565718586_613, 113 0.9677390756_7913913425_7347978784_337, 114 0.9372733924_0070590430_7758947710_209, 115 0.8972645323_4408190088_2509656454_496, 116 0.8482065834_1042721620_0648320774_217, 117 0.7904185014_4246593296_7649294817_947, 118 0.7244177313_6017004741_6186054613_938, 119 0.6509967412_9741697053_3735895313_275, 120 0.5709721726_0853884753_7226737253_911, 121 0.4850818636_4023968069_3655740232_351, 122 0.3941513470_7756336989_7207370981_045, 123 0.2991800071_5316881216_6780024266_389, 124 0.2011940939_9743452230_0628303394_596, 125 0.1011420669_1871749902_7074231447_392, 126 0.0000000000_0000000000_0000000000_000]; 127 // 128 immutable static Real[16] wgk_ = [ 129 0.0053774798_7292334898_7792051430_128, 130 0.0150079473_2931612253_8374763075_807, 131 0.0254608473_2671532018_6874001019_653, 132 0.0353463607_9137584622_2037948478_360, 133 0.0445897513_2476487660_8227299373_280, 134 0.0534815246_9092808726_5343147239_430, 135 0.0620095678_0067064028_5139230960_803, 136 0.0698541213_1872825870_9520077099_147, 137 0.0768496807_5772037889_4432777482_659, 138 0.0830805028_2313302103_8289247286_104, 139 0.0885644430_5621177064_7275443693_774, 140 0.0931265981_7082532122_5486872747_346, 141 0.0966427269_8362367850_5179907627_589, 142 0.0991735987_2179195933_2393173484_603, 143 0.1007698455_2387559504_4946662617_570, 144 0.1013300070_1479154901_7374792767_493]; 145 // 146 auto fv1 = dimension(fv1_.ptr, 15); 147 auto fv2 = dimension(fv2_.ptr, 15); 148 auto wg = dimension(wg_.ptr, 8); 149 auto wgk = dimension(wgk_.ptr, 16); 150 auto xgk = dimension(xgk_.ptr, 16); 151 // 152 // 153 // list of major variables 154 // ----------------------- 155 // centr - mid point of the interval 156 // hlgth - half-length of the interval 157 // absc - abscissa 158 // fval* - function value 159 // resg - result of the 15-point gauss formula 160 // resk - result of the 31-point kronrod formula 161 // reskh - approximation to the mean value of f over (a,b), 162 // i.e. to i/(b-a) 163 // 164 // machine dependent constants 165 // --------------------------- 166 // epmach is the largest relative spacing. 167 // uflow is the smallest positive magnitude. 168 //***first executable statement dqk31 169 epmach = Real.epsilon; 170 uflow = Real.min_normal; 171 // 172 centr = 0.5*(a+b); 173 hlgth = 0.5*(b-a); 174 dhlgth = fabs(hlgth); 175 // 176 // compute the 31-point kronrod approximation to 177 // the integral, and estimate the absolute error. 178 // 179 fc = f(centr); 180 resg = wg[8]*fc; 181 resk = wgk[16]*fc; 182 resabs = fabs(resk); 183 for (j=1; j<=7; j++) { //do 10 j=1,7 184 jtw = j*2; 185 absc = hlgth*xgk[jtw]; 186 fval1 = f(centr-absc); 187 fval2 = f(centr+absc); 188 fv1[jtw] = fval1; 189 fv2[jtw] = fval2; 190 fsum = fval1+fval2; 191 resg = resg+wg[j]*fsum; 192 resk = resk+wgk[jtw]*fsum; 193 resabs = resabs+wgk[jtw]*(fabs(fval1)+fabs(fval2)); 194 l10: ;} 195 for (j=1; j<=8; j++) { //do 15 j = 1,8 196 jtwm1 = j*2-1; 197 absc = hlgth*xgk[jtwm1]; 198 fval1 = f(centr-absc); 199 fval2 = f(centr+absc); 200 fv1[jtwm1] = fval1; 201 fv2[jtwm1] = fval2; 202 fsum = fval1+fval2; 203 resk = resk+wgk[jtwm1]*fsum; 204 resabs = resabs+wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); 205 l15: ;} 206 reskh = resk*0.5; 207 resasc = wgk[16]*fabs(fc-reskh); 208 for (j=1; j<=15; j++) { //do 20 j=1,15 209 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 210 l20: ;} 211 result = resk*hlgth; 212 resabs = resabs*dhlgth; 213 resasc = resasc*dhlgth; 214 abserr = fabs((resk-resg)*hlgth); 215 if(resasc != 0.0 && abserr != 0.0) 216 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 217 if(resabs > uflow/(0.5e2*epmach)) abserr = max 218 ((epmach*0.5e2)*resabs,abserr); 219 return; 220 } 221 222 version(unittest) import scid.core.testing; 223 unittest 224 { 225 alias qk31!(float, float delegate(float)) fqk31; 226 alias qk31!(double, double delegate(double)) dqk31; 227 alias qk31!(double, double function(double)) dfqk31; 228 alias qk31!(real, real delegate(real)) rqk31; 229 230 double f(double x) { return x^^3; } 231 double result, abserr, resabs, resasc; 232 qk31(&f, 0.0, 1.0, result, abserr, resabs, resasc); 233 assert (isAccurate(result, abserr, 0.25, 1e-6)); 234 }