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.qk51; 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 qk51(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 dqk51 25 //***date written 800101 (yymmdd) 26 //***revision date 830518 (yymmdd) 27 //***category no. h2a1a2 28 //***keywords 51-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 subroutine 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 51-point 57 // kronrod rule (resk) obtained by optimal addition 58 // of abscissae to the 25-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 dqk51 74 // 75 Real absc,centr,dhlgth, 76 epmach,fc,fsum,fval1,fval2,hlgth, 77 resg,resk,reskh,uflow; 78 Real[25] 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 51-point kronrod rule 86 // xgk(2), xgk(4), ... abscissae of the 25-point 87 // gauss rule 88 // xgk(1), xgk(3), ... abscissae which are optimally 89 // added to the 25-point gauss rule 90 // 91 // wgk - weights of the 51-point kronrod rule 92 // 93 // wg - weights of the 25-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[13] wg_ = [ 101 0.0113937985_0102628794_7902964113_235, 102 0.0263549866_1503213726_1901815295_299, 103 0.0409391567_0130631265_5623487711_646, 104 0.0549046959_7583519192_5936891540_473, 105 0.0680383338_1235691720_7187185656_708, 106 0.0801407003_3500101801_3234959669_111, 107 0.0910282619_8296364981_1497220702_892, 108 0.1005359490_6705064420_2206890392_686, 109 0.1085196244_7426365311_6093957050_117, 110 0.1148582591_4571164833_9325545869_556, 111 0.1194557635_3578477222_8178126512_901, 112 0.1222424429_9031004168_8959518945_852, 113 0.1231760537_2671545120_3902873079_050]; 114 // 115 static immutable Real[26] xgk_ = [ 116 0.9992621049_9260983419_3457486540_341, 117 0.9955569697_9049809790_8784946893_902, 118 0.9880357945_3407724763_7331014577_406, 119 0.9766639214_5951751149_8315386479_594, 120 0.9616149864_2584251241_8130033660_167, 121 0.9429745712_2897433941_4011169658_471, 122 0.9207471152_8170156174_6346084546_331, 123 0.8949919978_7827536885_1042006782_805, 124 0.8658470652_9327559544_8996969588_340, 125 0.8334426287_6083400142_1021108693_570, 126 0.7978737979_9850005941_0410904994_307, 127 0.7592592630_3735763057_7282865204_361, 128 0.7177664068_1308438818_6654079773_298, 129 0.6735663684_7346836448_5120633247_622, 130 0.6268100990_1031741278_8122681624_518, 131 0.5776629302_4122296772_3689841612_654, 132 0.5263252843_3471918259_9623778158_010, 133 0.4730027314_4571496052_2182115009_192, 134 0.4178853821_9303774885_1814394594_572, 135 0.3611723058_0938783773_5821730127_641, 136 0.3030895389_3110783016_7478909980_339, 137 0.2438668837_2098843204_5190362797_452, 138 0.1837189394_2104889201_5969888759_528, 139 0.1228646926_1071039638_7359818808_037, 140 0.0615444830_0568507888_6546392366_797, 141 0.0000000000_0000000000_0000000000_000]; 142 // 143 static immutable Real[26] wgk_ = [ 144 0.0019873838_9233031592_6507851882_843, 145 0.0055619321_3535671375_8040236901_066, 146 0.0094739733_8617415160_7207710523_655, 147 0.0132362291_9557167481_3656405846_976, 148 0.0168478177_0912829823_1516667536_336, 149 0.0204353711_4588283545_6568292235_939, 150 0.0240099456_0695321622_0092489164_881, 151 0.0274753175_8785173780_2948455517_811, 152 0.0307923001_6738748889_1109020215_229, 153 0.0340021302_7432933783_6748795229_551, 154 0.0371162714_8341554356_0330625367_620, 155 0.0400838255_0403238207_4839284467_076, 156 0.0428728450_2017004947_6895792439_495, 157 0.0455029130_4992178890_9870584752_660, 158 0.0479825371_3883671390_6392255756_915, 159 0.0502776790_8071567196_3325259433_440, 160 0.0523628858_0640747586_4366712137_873, 161 0.0542511298_8854549014_4543370459_876, 162 0.0559508112_2041231730_8240686382_747, 163 0.0574371163_6156783285_3582693939_506, 164 0.0586896800_2239420796_1974175856_788, 165 0.0597203403_2417405997_9099291932_562, 166 0.0605394553_7604586294_5360267517_565, 167 0.0611285097_1705304830_5859030416_293, 168 0.0614711898_7142531666_1544131965_264, 169 // note: wgk (26) was calculated from the values of wgk(1..25) 170 0.0615808180_6783293507_8759824240_066]; 171 // 172 auto fv1 = dimension(fv1_.ptr, 25); 173 auto fv2 = dimension(fv2_.ptr, 25); 174 auto xgk = dimension(xgk_.ptr, 26); 175 auto wgk = dimension(wgk_.ptr, 26); 176 auto wg = dimension(wg_.ptr, 13); 177 // 178 // 179 // list of major variables 180 // ----------------------- 181 // 182 // centr - mid point of the interval 183 // hlgth - half-length of the interval 184 // absc - abscissa 185 // fval* - function value 186 // resg - result of the 25-point gauss formula 187 // resk - result of the 51-point kronrod formula 188 // reskh - approximation to the mean value of f over (a,b), 189 // i.e. to i/(b-a) 190 // 191 // machine dependent constants 192 // --------------------------- 193 // 194 // epmach is the largest relative spacing. 195 // uflow is the smallest positive magnitude. 196 // 197 //***first executable statement dqk51 198 epmach = Real.epsilon; 199 uflow = Real.min_normal; 200 // 201 centr = 0.5*(a+b); 202 hlgth = 0.5*(b-a); 203 dhlgth = fabs(hlgth); 204 // 205 // compute the 51-point kronrod approximation to 206 // the integral, and estimate the absolute error. 207 // 208 fc = f(centr); 209 resg = wg[13]*fc; 210 resk = wgk[26]*fc; 211 resabs = fabs(resk); 212 for (j=1; j<=12; j++) { //do 10 j=1,12 213 jtw = j*2; 214 absc = hlgth*xgk[jtw]; 215 fval1 = f(centr-absc); 216 fval2 = f(centr+absc); 217 fv1[jtw] = fval1; 218 fv2[jtw] = fval2; 219 fsum = fval1+fval2; 220 resg = resg+wg[j]*fsum; 221 resk = resk+wgk[jtw]*fsum; 222 resabs = resabs+wgk[jtw]*(fabs(fval1)+fabs(fval2)); 223 l10: ;} 224 for (j=1; j<=13; j++) { //do 15 j = 1,13 225 jtwm1 = j*2-1; 226 absc = hlgth*xgk[jtwm1]; 227 fval1 = f(centr-absc); 228 fval2 = f(centr+absc); 229 fv1[jtwm1] = fval1; 230 fv2[jtwm1] = fval2; 231 fsum = fval1+fval2; 232 resk = resk+wgk[jtwm1]*fsum; 233 resabs = resabs+wgk[jtwm1]*(fabs(fval1)+fabs(fval2)); 234 l15: ;} 235 reskh = resk*0.5; 236 resasc = wgk[26]*fabs(fc-reskh); 237 for (j=1; j<=25; j++) { //do 20 j=1,25 238 resasc = resasc+wgk[j]*(fabs(fv1[j]-reskh)+fabs(fv2[j]-reskh)); 239 l20: ;} 240 result = resk*hlgth; 241 resabs = resabs*dhlgth; 242 resasc = resasc*dhlgth; 243 abserr = fabs((resk-resg)*hlgth); 244 if(resasc != 0.0 && abserr != 0.0) 245 abserr = resasc*min(0.1e1,((cast(Real)0.2e3)*abserr/resasc)^^(cast(Real)1.5)); 246 if(resabs > uflow/(0.5e2*epmach)) abserr = max 247 ((epmach*0.5e2)*resabs,abserr); 248 return; 249 } 250 251 version(unittest) import scid.core.testing; 252 unittest 253 { 254 alias qk51!(float, float delegate(float)) fqk51; 255 alias qk51!(double, double delegate(double)) dqk51; 256 alias qk51!(double, double function(double)) dfqk51; 257 alias qk51!(real, real delegate(real)) rqk51; 258 259 double f(double x) { return x^^3; } 260 double result, abserr, resabs, resasc; 261 qk51(&f, 0.0, 1.0, result, abserr, resabs, resasc); 262 assert (isAccurate(result, abserr, 0.25, 1e-6)); 263 }