1 /** D port of intde2.c from Takuya OOURA's DE-Quadrature package 2 ($(LINK http://www.kurims.kyoto-u.ac.jp/~ooura/)). 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.intde.intde2; 9 10 11 import std.math; 12 13 14 15 16 /** I = integral of f(x) over (a,b) 17 18 Examples: 19 --- 20 auto aw = new real[8000]; 21 intdeini(aw.length, tiny, eps, aw.ptr); // initialization of aw 22 ... 23 intde(f, a, b, aw.ptr, &i, &err); 24 --- 25 26 Params: 27 lenaw = length of aw 28 tiny = minimum value that 1/tiny does not 29 overflow 30 eps = relative error requested 31 aw = points and weights of the quadrature 32 formula, aw[0...lenaw-1] 33 f = integrand f(x) 34 a = lower limit of integration 35 b = upper limit of integration 36 i = approximation to the integral 37 err = estimate of the absolute error 38 39 Remarks: 40 <pre> 41 initial parameters 42 lenaw > 1000, 43 IEEE double : 44 lenaw = 8000; 45 tiny = 1.0e-307; 46 function 47 f(x) needs to be analytic over (a,b). 48 relative error 49 eps is relative error requested excluding 50 cancellation of significant digits. 51 i.e. eps means : (absolute error) / 52 (integral_a^b |f(x)| dx). 53 eps does not mean : (absolute error) / I. 54 error message 55 err >= 0 : normal termination. 56 err < 0 : abnormal termination. 57 i.e. convergent error is detected : 58 1. f(x) or (d/dx)^n f(x) has 59 discontinuous points or sharp 60 peaks over (a,b). 61 you must divide the interval 62 (a,b) at this points. 63 2. relative error of f(x) is 64 greater than eps. 65 3. f(x) has oscillatory factor 66 and frequency of the oscillation 67 is very high. 68 </pre> 69 */ 70 void intdeini(Real)(size_t lenaw, Real tiny, Real eps, Real *aw) 71 { 72 assert (lenaw > 1000); 73 /* ---- adjustable parameter ---- */ 74 enum : Real { EFS = 0.1, HOFF = 8.5 } 75 /* ------------------------------ */ 76 int nk, k, j; 77 Real tinyln, epsln, h0, ehp, ehm, h, t, ep, em, xw, wg; 78 79 tinyln = -log(tiny); 80 epsln = 1 - log(EFS * eps); 81 h0 = HOFF / epsln; 82 ehp = exp(h0); 83 ehm = 1 / ehp; 84 aw[2] = eps; 85 aw[3] = exp(-ehm * epsln); 86 aw[4] = sqrt(EFS * eps); 87 enum NOFF = 5; 88 aw[NOFF] = 0.5; 89 aw[NOFF + 1] = h0; 90 aw[NOFF + 2] = PI_4 * h0; 91 h = 2; 92 nk = 0; 93 k = NOFF + 3; 94 do { 95 t = h * 0.5; 96 do { 97 em = exp(h0 * t); 98 ep = PI_2 * em; 99 em = PI_2 / em; 100 j = k; 101 do { 102 xw = 1 / (1 + exp(ep - em)); 103 wg = xw * (1 - xw) * h0; 104 aw[j] = xw; 105 aw[j + 1] = wg * 4; 106 aw[j + 2] = wg * (ep + em); 107 ep *= ehp; 108 em *= ehm; 109 j += 3; 110 } while (ep < tinyln && j <= lenaw - 3); 111 t += h; 112 k += nk; 113 } while (t < 1); 114 h *= 0.5; 115 if (nk == 0) { 116 if (j > lenaw - 6) j -= 3; 117 nk = j - NOFF; 118 k += nk; 119 aw[1] = nk; 120 } 121 } while (2 * k - NOFF - 3 <= lenaw); 122 aw[0] = k - 3; 123 } 124 125 126 /// ditto 127 void intde(Real, Func)(Func f, Real a, Real b, Real *aw, 128 Real *i, Real *err) 129 { 130 int lenawm, nk, k, j, jtmp, jm, m, klim; 131 Real epsh, ba, ir, xa, fa, fb, errt, errh, errd, h, iback, irback; 132 133 enum NOFF = 5; 134 lenawm = cast(int) (aw[0] + 0.5); 135 nk = cast(int) (aw[1] + 0.5); 136 epsh = aw[4]; 137 ba = b - a; 138 *i = f((a + b) * aw[NOFF]); 139 ir = *i * aw[NOFF + 1]; 140 *i *= aw[NOFF + 2]; 141 *err = fabs(*i); 142 k = nk + NOFF; 143 j = NOFF; 144 do { 145 j += 3; 146 xa = ba * aw[j]; 147 fa = f(a + xa); 148 fb = f(b - xa); 149 ir += (fa + fb) * aw[j + 1]; 150 fa *= aw[j + 2]; 151 fb *= aw[j + 2]; 152 *i += fa + fb; 153 *err += fabs(fa) + fabs(fb); 154 } while (aw[j] > epsh && j < k); 155 errt = *err * aw[3]; 156 errh = *err * epsh; 157 errd = 1 + 2 * errh; 158 jtmp = j; 159 while (fabs(fa) > errt && j < k) { 160 j += 3; 161 fa = f(a + ba * aw[j]); 162 ir += fa * aw[j + 1]; 163 fa *= aw[j + 2]; 164 *i += fa; 165 } 166 jm = j; 167 j = jtmp; 168 while (fabs(fb) > errt && j < k) { 169 j += 3; 170 fb = f(b - ba * aw[j]); 171 ir += fb * aw[j + 1]; 172 fb *= aw[j + 2]; 173 *i += fb; 174 } 175 if (j < jm) jm = j; 176 jm -= NOFF + 3; 177 h = 1; 178 m = 1; 179 klim = k + nk; 180 while (errd > errh && klim <= lenawm) { 181 iback = *i; 182 irback = ir; 183 do { 184 jtmp = k + jm; 185 for (j = k + 3; j <= jtmp; j += 3) { 186 xa = ba * aw[j]; 187 fa = f(a + xa); 188 fb = f(b - xa); 189 ir += (fa + fb) * aw[j + 1]; 190 *i += (fa + fb) * aw[j + 2]; 191 } 192 k += nk; 193 j = jtmp; 194 do { 195 j += 3; 196 fa = f(a + ba * aw[j]); 197 ir += fa * aw[j + 1]; 198 fa *= aw[j + 2]; 199 *i += fa; 200 } while (fabs(fa) > errt && j < k); 201 j = jtmp; 202 do { 203 j += 3; 204 fb = f(b - ba * aw[j]); 205 ir += fb * aw[j + 1]; 206 fb *= aw[j + 2]; 207 *i += fb; 208 } while (fabs(fb) > errt && j < k); 209 } while (k < klim); 210 errd = h * (fabs(*i - 2 * iback) + fabs(ir - 2 * irback)); 211 h *= 0.5; 212 m *= 2; 213 klim = 2 * klim - NOFF; 214 } 215 *i *= h * ba; 216 if (errd > errh) { 217 *err = -errd * (m * fabs(ba)); 218 } else { 219 *err = *err * aw[2] * (m * fabs(ba)); 220 } 221 } 222 223 224 unittest 225 { 226 real f1(real x) { return 1/sqrt(x); } 227 real f2(real x) { return sqrt(4 - x*x); } 228 229 real[] aw = new real[8000]; 230 real i, err; 231 real tiny = 1e-307; 232 real eps = 1e-15; 233 234 intdeini(aw.length, tiny, eps, aw.ptr); 235 236 intde(&f1, 0.0L, 1.0L, aw.ptr, &i, &err); 237 assert (approxEqual(i, 2.0, eps)); 238 239 intde(&f2, 0.0L, 2.0L, aw.ptr, &i, &err); 240 assert (approxEqual(i, PI, eps)); 241 } 242 243 244 245 246 /** I = integral of f(x) over (a,infinity), 247 f(x) has not oscillatory factor. 248 249 Examples: 250 --- 251 auto aw = new real[8000]; 252 intdeiini(aw.length, tiny, eps, aw.ptr); // initialization of aw 253 ... 254 intdei(f, a, aw.ptr, &i, &err); 255 --- 256 257 Params: 258 lenaw = length of aw 259 tiny = minimum value that 1/tiny does not 260 overflow 261 eps = relative error requested 262 aw = points and weights of the quadrature 263 formula, aw[0...lenaw-1] 264 f = integrand f(x) 265 a = lower limit of integration 266 i = approximation to the integral 267 err = estimate of the absolute error 268 269 Remarks: 270 <pre> 271 initial parameters 272 lenaw > 1000, 273 IEEE double : 274 lenaw = 8000; 275 tiny = 1.0e-307; 276 function 277 f(x) needs to be analytic over (a,infinity). 278 relative error 279 eps is relative error requested excluding 280 cancellation of significant digits. 281 i.e. eps means : (absolute error) / 282 (integral_a^infinity |f(x)| dx). 283 eps does not mean : (absolute error) / I. 284 error message 285 err >= 0 : normal termination. 286 err < 0 : abnormal termination. 287 i.e. convergent error is detected : 288 1. f(x) or (d/dx)^n f(x) has 289 discontinuous points or sharp 290 peaks over (a,infinity). 291 you must divide the interval 292 (a,infinity) at this points. 293 2. relative error of f(x) is 294 greater than eps. 295 3. f(x) has oscillatory factor 296 and decay of f(x) is very slow 297 as x -> infinity. 298 </pre> 299 */ 300 void intdeiini(Real)(size_t lenaw, Real tiny, Real eps, Real *aw) 301 { 302 assert (lenaw > 1000); 303 /* ---- adjustable parameter ---- */ 304 enum : Real { EFS = 0.1, HOFF = 11.0 } 305 /* ------------------------------ */ 306 int nk, k, j; 307 Real tinyln, epsln, h0, ehp, ehm, h, t, ep, em, xp, xm, 308 wp, wm; 309 310 tinyln = -log(tiny); 311 epsln = 1 - log(EFS * eps); 312 h0 = HOFF / epsln; 313 ehp = exp(h0); 314 ehm = 1 / ehp; 315 aw[2] = eps; 316 aw[3] = exp(-ehm * epsln); 317 aw[4] = sqrt(EFS * eps); 318 enum NOFF = 5; 319 aw[NOFF] = 1; 320 aw[NOFF + 1] = 4 * h0; 321 aw[NOFF + 2] = PI_2 * h0; 322 h = 2; 323 nk = 0; 324 k = NOFF + 6; 325 do { 326 t = h * 0.5; 327 do { 328 em = exp(h0 * t); 329 ep = PI_4 * em; 330 em = PI_4 / em; 331 j = k; 332 do { 333 xp = exp(ep - em); 334 xm = 1 / xp; 335 wp = xp * ((ep + em) * h0); 336 wm = xm * ((ep + em) * h0); 337 aw[j] = xm; 338 aw[j + 1] = xp; 339 aw[j + 2] = xm * (4 * h0); 340 aw[j + 3] = xp * (4 * h0); 341 aw[j + 4] = wm; 342 aw[j + 5] = wp; 343 ep *= ehp; 344 em *= ehm; 345 j += 6; 346 } while (ep < tinyln && j <= lenaw - 6); 347 t += h; 348 k += nk; 349 } while (t < 1); 350 h *= 0.5; 351 if (nk == 0) { 352 if (j > lenaw - 12) j -= 6; 353 nk = j - NOFF; 354 k += nk; 355 aw[1] = nk; 356 } 357 } while (2 * k - NOFF - 6 <= lenaw); 358 aw[0] = k - 6; 359 } 360 361 362 /// ditto 363 void intdei(Real, Func)(Func f, Real a, Real *aw, Real *i, 364 Real *err) 365 { 366 int lenawm, nk, k, j, jtmp, jm, m, klim; 367 Real epsh, ir, fp, fm, errt, errh, errd, h, iback, irback; 368 369 enum NOFF = 5; 370 lenawm = cast(int) (aw[0] + 0.5); 371 nk = cast(int) (aw[1] + 0.5); 372 epsh = aw[4]; 373 *i = f(a + aw[NOFF]); 374 ir = *i * aw[NOFF + 1]; 375 *i *= aw[NOFF + 2]; 376 *err = fabs(*i); 377 k = nk + NOFF; 378 j = NOFF; 379 do { 380 j += 6; 381 fm = f(a + aw[j]); 382 fp = f(a + aw[j + 1]); 383 ir += fm * aw[j + 2] + fp * aw[j + 3]; 384 fm *= aw[j + 4]; 385 fp *= aw[j + 5]; 386 *i += fm + fp; 387 *err += fabs(fm) + fabs(fp); 388 } while (aw[j] > epsh && j < k); 389 errt = *err * aw[3]; 390 errh = *err * epsh; 391 errd = 1 + 2 * errh; 392 jtmp = j; 393 while (fabs(fm) > errt && j < k) { 394 j += 6; 395 fm = f(a + aw[j]); 396 ir += fm * aw[j + 2]; 397 fm *= aw[j + 4]; 398 *i += fm; 399 } 400 jm = j; 401 j = jtmp; 402 while (fabs(fp) > errt && j < k) { 403 j += 6; 404 fp = f(a + aw[j + 1]); 405 ir += fp * aw[j + 3]; 406 fp *= aw[j + 5]; 407 *i += fp; 408 } 409 if (j < jm) jm = j; 410 jm -= NOFF + 6; 411 h = 1; 412 m = 1; 413 klim = k + nk; 414 while (errd > errh && klim <= lenawm) { 415 iback = *i; 416 irback = ir; 417 do { 418 jtmp = k + jm; 419 for (j = k + 6; j <= jtmp; j += 6) { 420 fm = f(a + aw[j]); 421 fp = f(a + aw[j + 1]); 422 ir += fm * aw[j + 2] + fp * aw[j + 3]; 423 *i += fm * aw[j + 4] + fp * aw[j + 5]; 424 } 425 k += nk; 426 j = jtmp; 427 do { 428 j += 6; 429 fm = f(a + aw[j]); 430 ir += fm * aw[j + 2]; 431 fm *= aw[j + 4]; 432 *i += fm; 433 } while (fabs(fm) > errt && j < k); 434 j = jtmp; 435 do { 436 j += 6; 437 fp = f(a + aw[j + 1]); 438 ir += fp * aw[j + 3]; 439 fp *= aw[j + 5]; 440 *i += fp; 441 } while (fabs(fp) > errt && j < k); 442 } while (k < klim); 443 errd = h * (fabs(*i - 2 * iback) + fabs(ir - 2 * irback)); 444 h *= 0.5; 445 m *= 2; 446 klim = 2 * klim - NOFF; 447 } 448 *i *= h; 449 if (errd > errh) { 450 *err = -errd * m; 451 } else { 452 *err *= aw[2] * m; 453 } 454 } 455 456 457 unittest 458 { 459 real f3(real x) { return 1.0/(1 + x*x); } 460 real f4(real x) { return exp(-x)/sqrt(x); } 461 462 real[] aw = new real[8000]; 463 real i, err; 464 real tiny = 1e-307; 465 real eps = 1e-15; 466 467 intdeiini(aw.length, tiny, eps, aw.ptr); 468 469 intdei(&f3, 0.0L, aw.ptr, &i, &err); 470 assert (approxEqual(i, PI_2, eps)); 471 472 intdei(&f4, 0.0L, aw.ptr, &i, &err); 473 assert (approxEqual(i, sqrt(PI), eps)); 474 } 475 476 477 478 479 /** I = integral of f(x) over (a,infinity), 480 f(x) has oscillatory factor : 481 --- 482 f(x) = g(x) * sin(omega * x + theta) as x -> infinity 483 --- 484 485 Examples: 486 --- 487 auto aw = new real[8000]; 488 intdeoini(aw.length, tiny, eps, aw.ptr); // initialization of aw 489 ... 490 intdeo(f, a, omega, aw.ptr, &i, &err); 491 --- 492 493 Params: 494 lenaw = length of aw 495 tiny = minimum value that 1/tiny does not 496 overflow 497 eps = relative error requested 498 aw = points and weights of the quadrature 499 formula, aw[0...lenaw-1] 500 f = integrand f(x) 501 a = lower limit of integration 502 omega = frequency of oscillation 503 i = approximation to the integral 504 err = estimate of the absolute error 505 506 Remarks: 507 <pre> 508 initial parameters 509 lenaw > 1000, 510 IEEE double : 511 lenaw = 8000; 512 tiny = 1.0e-307; 513 function 514 f(x) needs to be analytic over (a,infinity). 515 relative error 516 eps is relative error requested excluding 517 cancellation of significant digits. 518 i.e. eps means : (absolute error) / 519 (integral_a^R |f(x)| dx). 520 eps does not mean : (absolute error) / I. 521 error message 522 err >= 0 : normal termination. 523 err < 0 : abnormal termination. 524 i.e. convergent error is detected : 525 1. f(x) or (d/dx)^n f(x) has 526 discontinuous points or sharp 527 peaks over (a,infinity). 528 you must divide the interval 529 (a,infinity) at this points. 530 2. relative error of f(x) is 531 greater than eps. 532 </pre> 533 */ 534 535 void intdeoini(Real)(size_t lenaw, Real tiny, Real eps, Real *aw) 536 { 537 assert (lenaw > 1000); 538 /* ---- adjustable parameter ---- */ 539 enum LMAX = 5; 540 enum : Real { EFS = 0.1, ENOFF = 0.40, PQOFF = 2.9, PPOFF = -0.72 } 541 /* ------------------------------ */ 542 int noff0, nk0, noff, k, nk, j; 543 Real tinyln, epsln, frq4, per2, pp, pq, ehp, ehm, h, t, 544 ep, em, tk, xw, wg, xa; 545 546 tinyln = -log(tiny); 547 epsln = 1 - log(EFS * eps); 548 frq4 = M_2_PI; 549 per2 = PI; 550 pq = PQOFF / epsln; 551 pp = PPOFF - log(pq * pq * frq4); 552 ehp = exp(2 * pq); 553 ehm = 1 / ehp; 554 aw[3] = LMAX; 555 aw[4] = eps; 556 aw[5] = sqrt(EFS * eps); 557 noff0 = 6; 558 nk0 = 1 + cast(int) (ENOFF * epsln); 559 aw[1] = nk0; 560 noff = 2 * nk0 + noff0; 561 wg = 0; 562 xw = 1; 563 for (k = 1; k <= nk0; k++) { 564 wg += xw; 565 aw[noff - 2 * k] = wg; 566 aw[noff - 2 * k + 1] = xw; 567 xw = xw * (nk0 - k) / k; 568 } 569 wg = per2 / wg; 570 for (k = noff0; k <= noff - 2; k += 2) { 571 aw[k] *= wg; 572 aw[k + 1] *= wg; 573 } 574 xw = exp(pp - PI_2); 575 aw[noff] = sqrt(xw * (per2 * 0.5)); 576 aw[noff + 1] = xw * pq; 577 aw[noff + 2] = per2 * 0.5; 578 h = 2; 579 nk = 0; 580 k = noff + 3; 581 do { 582 t = h * 0.5; 583 do { 584 em = exp(2 * pq * t); 585 ep = PI_4 * em; 586 em = PI_4 / em; 587 tk = t; 588 j = k; 589 do { 590 xw = exp(pp - ep - em); 591 wg = sqrt(frq4 * xw + tk * tk); 592 xa = xw / (tk + wg); 593 wg = (pq * xw * (ep - em) + xa) / wg; 594 aw[j] = xa; 595 aw[j + 1] = xw * pq; 596 aw[j + 2] = wg; 597 ep *= ehp; 598 em *= ehm; 599 tk += 1; 600 j += 3; 601 } while (ep < tinyln && j <= lenaw - 3); 602 t += h; 603 k += nk; 604 } while (t < 1); 605 h *= 0.5; 606 if (nk == 0) { 607 if (j > lenaw - 6) j -= 3; 608 nk = j - noff; 609 k += nk; 610 aw[2] = nk; 611 } 612 } while (2 * k - noff - 3 <= lenaw); 613 aw[0] = k - 3; 614 } 615 616 617 /// ditto 618 void intdeo(Real, Func)(Func f, Real a, Real omega, Real *aw, 619 Real *i, Real *err) 620 { 621 int lenawm, nk0, noff0, nk, noff, lmax, m, k, j, jm, l; 622 Real eps, per, perw, w02, ir, h, iback, irback, t, tk, 623 xa, fm, fp, errh, s0, s1, s2, errd; 624 625 lenawm = cast(int) (aw[0] + 0.5); 626 nk0 = cast(int) (aw[1] + 0.5); 627 noff0 = 6; 628 nk = cast(int) (aw[2] + 0.5); 629 noff = 2 * nk0 + noff0; 630 lmax = cast(int) (aw[3] + 0.5); 631 eps = aw[4]; 632 per = 1 / fabs(omega); 633 w02 = 2 * aw[noff + 2]; 634 perw = per * w02; 635 *i = f(a + aw[noff] * per); 636 ir = *i * aw[noff + 1]; 637 *i *= aw[noff + 2]; 638 *err = fabs(*i); 639 h = 2; 640 m = 1; 641 k = noff; 642 do { 643 iback = *i; 644 irback = ir; 645 t = h * 0.5; 646 do { 647 if (k == noff) { 648 tk = 1; 649 k += nk; 650 j = noff; 651 do { 652 j += 3; 653 xa = per * aw[j]; 654 fm = f(a + xa); 655 fp = f(a + xa + perw * tk); 656 ir += (fm + fp) * aw[j + 1]; 657 fm *= aw[j + 2]; 658 fp *= w02 - aw[j + 2]; 659 *i += fm + fp; 660 *err += fabs(fm) + fabs(fp); 661 tk += 1; 662 } while (aw[j] > eps && j < k); 663 errh = *err * aw[5]; 664 *err *= eps; 665 jm = j - noff; 666 } else { 667 tk = t; 668 for (j = k + 3; j <= k + jm; j += 3) { 669 xa = per * aw[j]; 670 fm = f(a + xa); 671 fp = f(a + xa + perw * tk); 672 ir += (fm + fp) * aw[j + 1]; 673 fm *= aw[j + 2]; 674 fp *= w02 - aw[j + 2]; 675 *i += fm + fp; 676 tk += 1; 677 } 678 j = k + jm; 679 k += nk; 680 } 681 while (fabs(fm) > *err && j < k) { 682 j += 3; 683 fm = f(a + per * aw[j]); 684 ir += fm * aw[j + 1]; 685 fm *= aw[j + 2]; 686 *i += fm; 687 } 688 fm = f(a + perw * tk); 689 s2 = w02 * fm; 690 *i += s2; 691 if (fabs(fp) > *err || fabs(s2) > *err) { 692 l = 0; 693 for (;;) { 694 l++; 695 s0 = 0; 696 s1 = 0; 697 s2 = fm * aw[noff0 + 1]; 698 for (j = noff0 + 2; j <= noff - 2; j += 2) { 699 tk += 1; 700 fm = f(a + perw * tk); 701 s0 += fm; 702 s1 += fm * aw[j]; 703 s2 += fm * aw[j + 1]; 704 } 705 if (s2 <= *err || l >= lmax) break; 706 *i += w02 * s0; 707 } 708 *i += s1; 709 if (s2 > *err) *err = s2; 710 } 711 t += h; 712 } while (t < 1); 713 if (m == 1) { 714 errd = 1 + 2 * errh; 715 } else { 716 errd = h * (fabs(*i - 2 * iback) + fabs(ir - 2 * irback)); 717 } 718 h *= 0.5; 719 m *= 2; 720 } while (errd > errh && 2 * k - noff <= lenawm); 721 *i *= h * per; 722 if (errd > errh) { 723 *err = -errd * per; 724 } else { 725 *err *= per * m * 0.5; 726 } 727 } 728 729 730 unittest 731 { 732 real f5(real x) { return sin(x)/x; } 733 real f6(real x) { return cos(x)/sqrt(x); } 734 735 real[] aw = new real[8000]; 736 real i, err; 737 real tiny = 1e-307; 738 real eps = 1e-15; 739 740 intdeoini(aw.length, tiny, eps, aw.ptr); 741 742 intdeo(&f5, 0.0L, 1.0L, aw.ptr, &i, &err); 743 assert (approxEqual(i, PI_2, eps)); 744 745 intdeo(&f6, 0.0L, 1.0L, aw.ptr, &i, &err); 746 assert (approxEqual(i, sqrt(PI_2), eps)); 747 }