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.qk61; 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 qk61(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 dqk61 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a1a2 28 //***keywords 61-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 dabs(f) over (a,b) 34 //***description 35 // 36 // integration rule 37 // standard fortran subroutine 38 // double precision version 39 // 40 // 41 // parameters 42 // on entry 43 // f - double precision 44 // function subprogram defining the integrand 45 // function f(x). the actual name for f needs to be 46 // declared e x t e r n a l in the calling program. 47 // 48 // a - double precision 49 // lower limit of integration 50 // 51 // b - double precision 52 // upper limit of integration 53 // 54 // on return 55 // result - double precision 56 // approximation to the integral i 57 // result is computed by applying the 61-point 58 // kronrod rule (resk) obtained by optimal addition of 59 // abscissae to the 30-point gauss rule (resg). 60 // 61 // abserr - double precision 62 // estimate of the modulus of the absolute error, 63 // which should equal or exceed dabs(i-result) 64 // 65 // resabs - double precision 66 // approximation to the integral j 67 // 68 // resasc - double precision 69 // approximation to the integral of dabs(f-i/(b-a)) 70 // 71 // 72 //***references (none) 73 //***routines called d1mach 74 //***end prologue dqk61 75 // 76 Real dabsc,centr,dhlgth, 77 epmach,fc,fsum,fval1,fval2,hlgth, 78 resg,resk,reskh,uflow; 79 Real[30] fv1_, fv2_; 80 int j,jtw,jtwm1; 81 // 82 // the abscissae and weights are given for the 83 // interval (-1,1). because of symmetry only the positive 84 // abscissae and their corresponding weights are given. 85 // 86 // xgk - abscissae of the 61-point kronrod rule 87 // xgk(2), xgk(4) ... abscissae of the 30-point 88 // gauss rule 89 // xgk(1), xgk(3) ... optimally added abscissae 90 // to the 30-point gauss rule 91 // 92 // wgk - weights of the 61-point kronrod rule 93 // 94 // wg - weigths of the 30-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[15] wg_ = [ 102 0.0079681924_9616660561_5465883474_674, 103 0.0184664683_1109095914_2302131912_047, 104 0.0287847078_8332336934_9719179611_292, 105 0.0387991925_6962704959_6801936446_348, 106 0.0484026728_3059405290_2938140422_808, 107 0.0574931562_1761906648_1721689402_056, 108 0.0659742298_8218049512_8128515115_962, 109 0.0737559747_3770520626_8243850022_191, 110 0.0807558952_2942021535_4694938460_530, 111 0.0868997872_0108297980_2387530715_126, 112 0.0921225222_3778612871_7632707087_619, 113 0.0963687371_7464425963_9468626351_810, 114 0.0995934205_8679526706_2780282103_569, 115 0.1017623897_4840550459_6428952168_554, 116 0.1028526528_9355884034_1285636705_415]; 117 // 118 static immutable Real[31] xgk_ = [ 119 0.9994844100_5049063757_1325895705_811, 120 0.9968934840_7464954027_1630050918_695, 121 0.9916309968_7040459485_8628366109_486, 122 0.9836681232_7974720997_0032581605_663, 123 0.9731163225_0112626837_4693868423_707, 124 0.9600218649_6830751221_6871025581_798, 125 0.9443744447_4855997941_5831324037_439, 126 0.9262000474_2927432587_9324277080_474, 127 0.9055733076_9990779854_6522558925_958, 128 0.8825605357_9205268154_3116462530_226, 129 0.8572052335_4606109895_8658510658_944, 130 0.8295657623_8276839744_2898119732_502, 131 0.7997278358_2183908301_3668942322_683, 132 0.7677774321_0482619491_7977340974_503, 133 0.7337900624_5322680472_6171131369_528, 134 0.6978504947_9331579693_2292388026_640, 135 0.6600610641_2662696137_0053668149_271, 136 0.6205261829_8924286114_0477556431_189, 137 0.5793452358_2636169175_6024932172_540, 138 0.5366241481_4201989926_4169793311_073, 139 0.4924804678_6177857499_3693061207_709, 140 0.4470337695_3808917678_0609900322_854, 141 0.4004012548_3039439253_5476211542_661, 142 0.3527047255_3087811347_1037207089_374, 143 0.3040732022_7362507737_2677107199_257, 144 0.2546369261_6788984643_9805129817_805, 145 0.2045251166_8230989143_8957671002_025, 146 0.1538699136_0858354696_3794672743_256, 147 0.1028069379_6673703014_7096751318_001, 148 0.0514718425_5531769583_3025213166_723, 149 0.0000000000_0000000000_0000000000_000]; 150 // 151 static immutable Real[31] wgk_ = [ 152 0.0013890136_9867700762_4551591226_760, 153 0.0038904611_2709988405_1267201844_516, 154 0.0066307039_1593129217_3319826369_750, 155 0.0092732796_5951776342_8441146892_024, 156 0.0118230152_5349634174_2232898853_251, 157 0.0143697295_0704580481_2451432443_580, 158 0.0169208891_8905327262_7572289420_322, 159 0.0194141411_9394238117_3408951050_128, 160 0.0218280358_2160919229_7167485738_339, 161 0.0241911620_7808060136_5686370725_232, 162 0.0265099548_8233310161_0601709335_075, 163 0.0287540487_6504129284_3978785354_334, 164 0.0309072575_6238776247_2884252943_092, 165 0.0329814470_5748372603_1814191016_854, 166 0.0349793380_2806002413_7499670731_468, 167 0.0368823646_5182122922_3911065617_136, 168 0.0386789456_2472759295_0348651532_281, 169 0.0403745389_5153595911_1995279752_468, 170 0.0419698102_1516424614_7147541285_970, 171 0.0434525397_0135606931_6831728117_073, 172 0.0448148001_3316266319_2355551616_723, 173 0.0460592382_7100698811_6271735559_374, 174 0.0471855465_6929915394_5261478181_099, 175 0.0481858617_5708712914_0779492298_305, 176 0.0490554345_5502977888_7528165367_238, 177 0.0497956834_2707420635_7811569379_942, 178 0.0504059214_0278234684_0893085653_585, 179 0.0508817958_9874960649_2297473049_805, 180 0.0512215478_4925877217_0656282604_944, 181 0.0514261285_3745902593_3862879215_781, 182 0.0514947294_2945156755_8340433647_099]; 183 // 184 auto fv1 = dimension(fv1_.ptr, 30); 185 auto fv2 = dimension(fv2_.ptr, 30); 186 auto xgk = dimension(xgk_.ptr, 31); 187 auto wgk = dimension(wgk_.ptr, 31); 188 auto wg = dimension(wg_.ptr, 15); 189 // 190 // list of major variables 191 // ----------------------- 192 // 193 // centr - mid point of the interval 194 // hlgth - half-length of the interval 195 // dabsc - abscissa 196 // fval* - function value 197 // resg - result of the 30-point gauss rule 198 // resk - result of the 61-point kronrod rule 199 // reskh - approximation to the mean value of f 200 // over (a,b), i.e. to i/(b-a) 201 // 202 // machine dependent constants 203 // --------------------------- 204 // 205 // epmach is the largest relative spacing. 206 // uflow is the smallest positive magnitude. 207 // 208 epmach = Real.epsilon; 209 uflow = Real.min_normal; 210 // 211 centr = 0.5*(b+a); 212 hlgth = 0.5*(b-a); 213 dhlgth = fabs(hlgth); 214 // 215 // compute the 61-point kronrod approximation to the 216 // integral, and estimate the absolute error. 217 // 218 //***first executable statement dqk61 219 resg = 0.0; 220 fc = f(centr); 221 resk = wgk[31]*fc; 222 resabs = fabs(resk); 223 for (j=1; j<=15; j++) { //do 10 j=1,15 224 jtw = j*2; 225 dabsc = hlgth*xgk[jtw]; 226 fval1 = f(centr-dabsc); 227 fval2 = f(centr+dabsc); 228 fv1[jtw] = fval1; 229 fv2[jtw] = fval2; 230 fsum = fval1+fval2; 231 resg = resg+wg[j]*fsum; 232 resk = resk+wgk[jtw]*fsum; 233 resabs = resabs+wgk[jtw]*(fabs(fval1)+fabs(fval2)); 234 l10: ;} 235 for (j=1; j<=15; j++) { //do 15 j=1,15 236 jtwm1 = j*2-1; 237 dabsc = hlgth*xgk[jtwm1]; 238 fval1 = f(centr-dabsc); 239 fval2 = f(centr+dabsc); 240 fv1[jtwm1] = fval1; 241 fv2[jtwm1] = fval2; 242 fsum = fval1+fval2; 243 resk = resk+wgk[jtwm1]*fsum; 244 resabs = resabs+wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); 245 l15: ;} 246 reskh = resk*0.5; 247 resasc = wgk[31]*fabs(fc-reskh); 248 for (j=1; j<=30; j++) { //do 20 j=1,30 249 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 250 l20: ;} 251 result = resk*hlgth; 252 resabs = resabs*dhlgth; 253 resasc = resasc*dhlgth; 254 abserr = fabs((resk-resg)*hlgth); 255 if(resasc != 0.0 && abserr != 0.0) 256 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 257 if(resabs > uflow/(0.5e2*epmach)) abserr = max 258 ((epmach*0.5e2)*resabs,abserr); 259 return; 260 } 261 262 version(unittest) import scid.core.testing; 263 unittest 264 { 265 alias qk61!(float, float delegate(float)) fqk61; 266 alias qk61!(double, double delegate(double)) dqk61; 267 alias qk61!(double, double function(double)) dfqk61; 268 alias qk61!(real, real delegate(real)) rqk61; 269 270 double f(double x) { return x^^3; } 271 double result, abserr, resabs, resasc; 272 qk61(&f, 0.0, 1.0, result, abserr, resabs, resasc); 273 assert (isAccurate(result, abserr, 0.25, 1e-6)); 274 }