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.qc25s; 9 10 11 import std.math: fabs, log; 12 13 import scid.core.fortran; 14 import scid.ports.quadpack.qcheb; 15 import scid.ports.quadpack.qk15w; 16 import scid.ports.quadpack.qwgts; 17 18 19 20 21 /// 22 void qc25s(Real, Func)(Func f, Real a, Real b, Real bl, Real br, Real alfa, 23 Real beta, Real* ri_, Real* rj_, Real* rg_, Real* rh_, out Real result, 24 out Real abserr, out Real resasc, int integr, out int nev) 25 { 26 //***begin prologue dqc25s 27 //***date written 810101 (yymmdd) 28 //***revision date 830518 (yymmdd) 29 //***category no. h2a2a2 30 //***keywords 25-point clenshaw-curtis integration 31 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 32 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 33 //***purpose to compute i = integral of f*w over (bl,br), with error 34 // estimate, where the weight function w has a singular 35 // behaviour of algebraico-logarithmic type at the points 36 // a and/or b. (bl,br) is a part of (a,b). 37 //***description 38 // 39 // integration rules for integrands having algebraico-logarithmic 40 // end point singularities 41 // standard fortran subroutine 42 // double precision version 43 // 44 // parameters 45 // f - double precision 46 // function subprogram defining the integrand 47 // f(x). the actual name for f needs to be declared 48 // e x t e r n a l in the driver program. 49 // 50 // a - double precision 51 // left end point of the original interval 52 // 53 // b - double precision 54 // right end point of the original interval, b.gt.a 55 // 56 // bl - double precision 57 // lower limit of integration, bl.ge.a 58 // 59 // br - double precision 60 // upper limit of integration, br.le.b 61 // 62 // alfa - double precision 63 // parameter in the weight function 64 // 65 // beta - double precision 66 // parameter in the weight function 67 // 68 // ri,rj,rg,rh - double precision 69 // modified chebyshev moments for the application 70 // of the generalized clenshaw-curtis 71 // method (computed in subroutine dqmomo) 72 // 73 // result - double precision 74 // approximation to the integral 75 // result is computed by using a generalized 76 // clenshaw-curtis method if b1 = a or br = b. 77 // in all other cases the 15-point kronrod 78 // rule is applied, obtained by optimal addition of 79 // abscissae to the 7-point gauss rule. 80 // 81 // abserr - double precision 82 // estimate of the modulus of the absolute error, 83 // which should equal or exceed abs(i-result) 84 // 85 // resasc - double precision 86 // approximation to the integral of abs(f*w-i/(b-a)) 87 // 88 // integr - integer 89 // which determines the weight function 90 // = 1 w(x) = (x-a)**alfa*(b-x)**beta 91 // = 2 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a) 92 // = 3 w(x) = (x-a)**alfa*(b-x)**beta*log(b-x) 93 // = 4 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)* 94 // log(b-x) 95 // 96 // nev - integer 97 // number of integrand evaluations 98 //***references (none) 99 //***routines called dqcheb,dqk15w 100 //***end prologue dqc25s 101 // 102 Real centr,dc,factor,fix,hlgth,resabs,res12, 103 res24,u; 104 Real[13] cheb12_; 105 Real[25] cheb24_, fval_; 106 int i=1,isym=1; 107 // 108 // the vector x contains the values cos(k*pi/24) 109 // k = 1, ..., 11, to be used for the computation of the 110 // chebyshev series expansion of f. 111 // 112 static immutable Real[11] x_ = [ 113 0.9914448613_7381041114_4557526928_563, 114 0.9659258262_8906828674_9743199728_897, 115 0.9238795325_1128675612_8183189396_788, 116 0.8660254037_8443864676_3723170752_936, 117 0.7933533402_9123516457_9776961501_299, 118 0.7071067811_8654752440_0844362104_849, 119 0.6087614290_0872063941_6097542898_164, 120 0.5000000000_0000000000_0000000000_000, 121 0.3826834323_6508977172_8459984030_399, 122 0.2588190451_0252076234_8898837624_048, 123 0.1305261922_2005159154_8406227895_489]; 124 // 125 auto cheb12 = dimension(cheb12_.ptr, 13); 126 auto cheb24 = dimension(cheb24_.ptr, 25); 127 auto fval = dimension(fval_.ptr, 25); 128 auto rg = dimension(rg_, 25); 129 auto rh = dimension(rh_, 25); 130 auto ri = dimension(ri_, 25); 131 auto rj = dimension(rj_, 25); 132 auto x = dimension(x_.ptr, 11); 133 // 134 // list of major variables 135 // ----------------------- 136 // 137 // fval - value of the function f at the points 138 // (br-bl)*0.5*cos(k*pi/24)+(br+bl)*0.5 139 // k = 0, ..., 24 140 // cheb12 - coefficients of the chebyshev series expansion 141 // of degree 12, for the function f, in the 142 // interval (bl,br) 143 // cheb24 - coefficients of the chebyshev series expansion 144 // of degree 24, for the function f, in the 145 // interval (bl,br) 146 // res12 - approximation to the integral obtained from cheb12 147 // res24 - approximation to the integral obtained from cheb24 148 // dqwgts - external function subprogram defining 149 // the four possible weight functions 150 // hlgth - half-length of the interval (bl,br) 151 // centr - mid point of the interval (bl,br) 152 // 153 //***first executable statement dqc25s 154 nev = 25; 155 if(bl == a && (alfa != 0.0 || integr == 2 || integr == 4)) 156 goto l10; 157 if(br == b && (beta != 0.0 || integr == 3 || integr == 4)) 158 goto l140; 159 // 160 // if a > bl and b < br, apply the 15-point gauss-kronrod 161 // scheme. 162 // 163 // 164 qk15w!(Real,Func)(f,&qwgts!Real,a,b,alfa,beta,integr,bl,br, 165 result,abserr,resabs,resasc); 166 nev = 15; 167 goto l270; 168 // 169 // this part of the program is executed only if a = bl. 170 // ---------------------------------------------------- 171 // 172 // compute the chebyshev series expansion of the 173 // following function 174 // f1 = (0.5*(b+b-br-a)-0.5*(br-a)*x)**beta 175 // *f(0.5*(br-a)*x+0.5*(br+a)) 176 // 177 l10: hlgth = 0.5*(br-bl); 178 centr = 0.5*(br+bl); 179 fix = b-centr; 180 fval[1] = 0.5*f(hlgth+centr)*((fix-hlgth)^^beta); 181 fval[13] = f(centr)*(fix^^beta); 182 fval[25] = 0.5*f(centr-hlgth)*((fix+hlgth)^^beta); 183 for (i=2; i<=12; i++) { //do 20 i=2,12 184 u = hlgth*x[i-1]; 185 isym = 26-i; 186 fval[i] = f(u+centr)*((fix-u)^^beta); 187 fval[isym] = f(centr-u)*((fix+u)^^beta); 188 l20: ;} 189 factor = hlgth^^(alfa+(cast(Real)0.1e1)); 190 result = 0.0; 191 abserr = 0.0; 192 res12 = 0.0; 193 res24 = 0.0; 194 if(integr > 2) goto l70; 195 qcheb!Real(x.ptr,fval.ptr,cheb12.ptr,cheb24.ptr); 196 // 197 // integr = 1 (or 2) 198 // 199 for (i=1; i<=13; i++) { //do 30 i=1,13 200 res12 = res12+cheb12[i]*ri[i]; 201 res24 = res24+cheb24[i]*ri[i]; 202 l30: ;} 203 for (i=14; i<=25; i++) { //do 40 i=14,25 204 res24 = res24+cheb24[i]*ri[i]; 205 l40: ;} 206 if(integr == 1) goto l130; 207 // 208 // integr = 2 209 // 210 dc = log(br-bl); 211 result = res24*dc; 212 abserr = fabs((res24-res12)*dc); 213 res12 = 0.0; 214 res24 = 0.0; 215 for (i=1; i<=13; i++) { //do 50 i=1,13 216 res12 = res12+cheb12[i]*rg[i]; 217 res24 = res12+cheb24[i]*rg[i]; 218 l50: ;} 219 for (i=14; i<=25; i++) { //do 60 i=14,25 220 res24 = res24+cheb24[i]*rg[i]; 221 l60: ;} 222 goto l130; 223 // 224 // compute the chebyshev series expansion of the 225 // following function 226 // f4 = f1*log(0.5*(b+b-br-a)-0.5*(br-a)*x) 227 // 228 l70: fval[1] = fval[1]*log(fix-hlgth); 229 fval[13] = fval[13]*log(fix); 230 fval[25] = fval[25]*log(fix+hlgth); 231 for (i=2; i<=12; i++) { //do 80 i=2,12 232 u = hlgth*x[i-1]; 233 isym = 26-i; 234 fval[i] = fval[i]*log(fix-u); 235 fval[isym] = fval[isym]*log(fix+u); 236 l80: ;} 237 qcheb!Real(x.ptr,fval.ptr,cheb12.ptr,cheb24.ptr); 238 // 239 // integr = 3 (or 4) 240 // 241 for (i=1; i<=13; i++) { //do 90 i=1,13 242 res12 = res12+cheb12[i]*ri[i]; 243 res24 = res24+cheb24[i]*ri[i]; 244 l90: ;} 245 for (i=14; i<=25; i++) { //do 100 i=14,25 246 res24 = res24+cheb24[i]*ri[i]; 247 l100: ;} 248 if(integr == 3) goto l130; 249 // 250 // integr = 4 251 // 252 dc = log(br-bl); 253 result = res24*dc; 254 abserr = fabs((res24-res12)*dc); 255 res12 = 0.0; 256 res24 = 0.0; 257 for (i=1; i<=13; i++) { //do 110 i=1,13 258 res12 = res12+cheb12[i]*rg[i]; 259 res24 = res24+cheb24[i]*rg[i]; 260 l110: ;} 261 for (i=14; i<=25; i++) { //do 120 i=14,25 262 res24 = res24+cheb24[i]*rg[i]; 263 l120: ;} 264 l130: result = (result+res24)*factor; 265 abserr = (abserr+fabs(res24-res12))*factor; 266 goto l270; 267 // 268 // this part of the program is executed only if b = br. 269 // ---------------------------------------------------- 270 // 271 // compute the chebyshev series expansion of the 272 // following function 273 // f2 = (0.5*(b+bl-a-a)+0.5*(b-bl)*x)**alfa 274 // *f(0.5*(b-bl)*x+0.5*(b+bl)) 275 // 276 l140: hlgth = 0.5*(br-bl); 277 centr = 0.5*(br+bl); 278 fix = centr-a; 279 fval[1] = 0.5*f(hlgth+centr)*((fix+hlgth)^^alfa); 280 fval[13] = f(centr)*(fix^^alfa); 281 fval[25] = 0.5*f(centr-hlgth)*((fix-hlgth)^^alfa); 282 for (i=2; i<=12; i++) { //do 150 i=2,12 283 u = hlgth*x[i-1]; 284 isym = 26-i; 285 fval[i] = f(u+centr)*((fix+u)^^alfa); 286 fval[isym] = f(centr-u)*((fix-u)^^alfa); 287 l150: ;} 288 factor = hlgth^^(beta+(cast(Real)0.1e1)); 289 result = 0.0; 290 abserr = 0.0; 291 res12 = 0.0; 292 res24 = 0.0; 293 if(integr == 2 || integr == 4) goto l200; 294 // 295 // integr = 1 (or 3) 296 // 297 qcheb!Real(x.ptr,fval.ptr,cheb12.ptr,cheb24.ptr); 298 for (i=1; i<=13; i++) { //do 160 i=1,13 299 res12 = res12+cheb12[i]*rj[i]; 300 res24 = res24+cheb24[i]*rj[i]; 301 l160: ;} 302 for (i=14; i<=25; i++) { //do 170 i=14,25 303 res24 = res24+cheb24[i]*rj[i]; 304 l170: ;} 305 if(integr == 1) goto l260; 306 // 307 // integr = 3 308 // 309 dc = log(br-bl); 310 result = res24*dc; 311 abserr = fabs((res24-res12)*dc); 312 res12 = 0.0; 313 res24 = 0.0; 314 for (i=1; i<=13; i++) { //do 180 i=1,13 315 res12 = res12+cheb12[i]*rh[i]; 316 res24 = res24+cheb24[i]*rh[i]; 317 l180: ;} 318 for (i=14; i<=25; i++) { //do 190 i=14,25 319 res24 = res24+cheb24[i]*rh[i]; 320 l190: ;} 321 goto l260; 322 // 323 // compute the chebyshev series expansion of the 324 // following function 325 // f3 = f2*log(0.5*(b-bl)*x+0.5*(b+bl-a-a)) 326 // 327 l200: fval[1] = fval[1]*log(hlgth+fix); 328 fval[13] = fval[13]*log(fix); 329 fval[25] = fval[25]*log(fix-hlgth); 330 for (i=2; i<=12; i++) { //do 210 i=2,12 331 u = hlgth*x[i-1]; 332 isym = 26-i; 333 fval[i] = fval[i]*log(u+fix); 334 fval[isym] = fval[isym]*log(fix-u); 335 l210: ;} 336 qcheb!Real(x.ptr,fval.ptr,cheb12.ptr,cheb24.ptr); 337 // 338 // integr = 2 (or 4) 339 // 340 for (i=1; i<=13; i++) { //do 220 i=1,13 341 res12 = res12+cheb12[i]*rj[i]; 342 res24 = res24+cheb24[i]*rj[i]; 343 l220: ;} 344 for (i=14; i<=25; i++) { //do 230 i=14,25 345 res24 = res24+cheb24[i]*rj[i]; 346 l230: ;} 347 if(integr == 2) goto l260; 348 dc = log(br-bl); 349 result = res24*dc; 350 abserr = fabs((res24-res12)*dc); 351 res12 = 0.0; 352 res24 = 0.0; 353 // 354 // integr = 4 355 // 356 for (i=1; i<=13; i++) { //do 240 i=1,13 357 res12 = res12+cheb12[i]*rh[i]; 358 res24 = res24+cheb24[i]*rh[i]; 359 l240: ;} 360 for (i=14; i<=25; i++) { //do 250 i=14,25 361 res24 = res24+cheb24[i]*rh[i]; 362 l250: ;} 363 l260: result = (result+res24)*factor; 364 abserr = (abserr+fabs(res24-res12))*factor; 365 l270: return; 366 } 367 368 369 unittest 370 { 371 alias qc25s!(float, float delegate(float)) fqc25s; 372 alias qc25s!(double, double delegate(double)) dqc25s; 373 alias qc25s!(double, double function(double)) dfqc25s; 374 alias qc25s!(real, real delegate(real)) rqc25s; 375 }