1/* 2 * specialfunctions - special functions (e.g.: gamma, zeta, psi) 3 * 4 * Copyright (C) 2013,2021 Christoph Zurnieden 5 * 6 * Calc is open software; you can redistribute it and/or modify 7 * it under the terms of the version 2.1 of the GNU Lesser General Public 8 * License as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser 13 * General Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 * 20 * Under source code control: 2013/08/11 01:31:28 21 * File existed as early as: 2013 22 */ 23 24 25/* 26 * hide internal function from resource debugging 27 */ 28static resource_debug_level; 29resource_debug_level = config("resource_debug", 0); 30 31 32/* 33 get dependencies 34*/ 35read -once zeta2; 36 37 38/* 39 zeta2(x,s) is in the extra file "zeta2.cal" because of a different license: 40 GPL instead of LGPL. 41*/ 42define zeta(s) 43{ 44 /* Not the best way, I'm afraid, but a way. */ 45 return hurwitzzeta(s, 1); 46} 47 48define psi0(z) 49{ 50 local i k m x y eps_digits eps ret; 51 52 /* 53 * One can use the Stirling series, too, which might be faster for some 54 * values. The series used here converges very fast but needs a lot of 55 * bernoulli numbers which are quite expensive to compute. 56 */ 57 58 eps = epsilon(); 59 epsilon(eps * 1e-3); 60 if (isint(z) && z <= 0) { 61 return newerror("psi(z); z is a negative integer or 0"); 62 } 63 if (re(z) < 0) { 64 return (pi() * cot(pi() * (0 - z))) + psi0(1 - z); 65 } 66 eps_digits = digits(1 / epsilon()); 67 /* 68 * R.W. Gosper has r = .41 as the relation, empirical tests showed that 69 * for d<100 r = .375 is sufficient and r = .396 for d<200. 70 * It does not save much time but for a long series even a little 71 * can grow large. 72 */ 73 if (eps_digits < 100) { 74 k = 2 * ceil(.375 * eps_digits); 75 } else if (eps_digits < 200) { 76 k = 2 * ceil(.395 * eps_digits); 77 } else { 78 k = 2 * ceil(11 / 27 * eps_digits); 79 } 80 m = 0; 81 y = (z + k) ^ 2; 82 x = 0.0; 83 /* 84 * There is a chance to speed up the first partial sum with binary 85 * splitting. The second partial sum is dominated by the calculation 86 * of the Bernoulli numbers but can profit from binary splitting when 87 * the Bernoulli numbers are already cached. 88 */ 89 for (i = 1; i <= (k / 2); i++) { 90 m = 1 / (z + 2 * i - 1) + 1 / (z + 2 * i - 2) + m; 91 x = (x + bernoulli(k - 2 * i + 2) / (k - 2 * i + 2)) / y; 92 } 93 ret = ln(z + k) - 1 / (2 * (z + k)) - x - m; 94 epsilon(eps); 95 return ret; 96} 97 98define psi(z) 99{ 100 return psi0(z); 101} 102 103define polygamma(m, z) 104{ 105 /* 106 * TODO: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0002/ 107 */ 108 if (!isint(m)) { 109 return newerror("polygamma(m,z); m not an integer"); 110 } 111 if (m < 0) { 112 return newerror("polygamma(m,z); m is < 0"); 113 } 114 /* 115 * Reflection formula not implemented yet, needs cot-differentiation 116 * http://functions.wolfram.com/ElementaryFunctions/Cot/20/02/0003/ 117 * which is not implemented yet. 118 */ 119 if (m == 0) { 120 return psi0(z); 121 } 122 /* 123 * use factorial for m a (small) integer? 124 * use lngamma for m large? 125 */ 126 if (isodd(m + 1)) { 127 return (-1) * gamma(m + 1) * hurwitzzeta(m + 1, z) 128 } 129 return gamma(m + 1) * hurwitzzeta(m + 1, z); 130} 131 132/* 133 Cache for the variable independent coefficients in the sum for the 134 Gamma-function. 135*/ 136static __CZ__Ck; 137/* 138 Log-Gamma function for Re(z) > 0. 139*/ 140define __CZ__lngammarp(z) 141{ 142 local epsilon accuracy a factrl sum n ret holds_enough term; 143 144 epsilon = epsilon(); 145 accuracy = digits(int (1 / epsilon)) + 3; 146 147 epsilon(1e-18); 148 a = ceil(1.252850440912568095810522965 * accuracy); 149 150 epsilon(epsilon * 10 ^ (-(digits(1 / epsilon) //2))); 151 152 holds_enough = 0; 153 154 if (size(__CZ__Ck) != a) { 155 __CZ__Ck = mat[a]; 156 holds_enough = 1; 157 } 158 159 factrl = 1.0; 160 161 __CZ__Ck[0] = sqrt(2 * pi()); /* c_0 */ 162 for (n = 1; n < a; n++) { 163 if (holds_enough == 1) { 164 __CZ__Ck[n] = (a - n) ^ (n - 0.5) * exp(a - n) / factrl; 165 factrl *= -n} 166 } 167 sum = __CZ__Ck[0]; 168 for (n = 1; n < a; n++) { 169 sum += __CZ__Ck[n] / (z + n); 170 } 171 172 ret = ln(sum) + (-(z + a)) + ln(z + a) * (z + 0.5); 173 ret = ret - ln(z); 174 175 /* 176 * Will take some time for large im(z) but almost all time is spend above 177 * in that case. 178 */ 179 if (im(ret)) { 180 ret = re(ret) + ln(exp(im(ret) * 1i)); 181 } 182 183 epsilon(epsilon); 184 return ret; 185} 186 187/* Simple lngamma with low precision*/ 188define __CZ__lngamma_lanczos(z) 189{ 190 local a k g zghalf lanczos; 191 mat lanczos[15] = { 192 9.9999999999999709182e-1, 193 5.7156235665862923516e1, 194 -5.9597960355475491248e1, 195 1.4136097974741747173e1, 196 -4.9191381609762019978e-1, 197 3.3994649984811888638e-5, 198 4.6523628927048576010e-5, 199 -9.8374475304879566105e-5, 200 1.5808870322491249322e-4, 201 -2.1026444172410489480e-4, 202 2.1743961811521265523e-4, 203 -1.6431810653676390482e-4, 204 8.4418223983852751308e-5, 205 -2.6190838401581411237e-5, 206 3.6899182659531626821e-6 207 }; 208 g = 607 / 128; 209 z = z - 1; 210 a = 0; 211 for (k = 12; k >= 1; k--) { 212 a += lanczos[k] / (z + k); 213 } 214 a += lanczos[0]; 215 zghalf = z + g + .5; 216 return (ln(sqrt(2 * pi())) + ln(a) - zghalf) + (z + .5) * ln(zghalf); 217} 218 219/* Prints the Spouge coefficients actually in use. */ 220define __CZ__print__CZ__Ck() 221{ 222 local k; 223 if (size(__CZ__Ck) <= 1) { 224 __CZ__lngammarp(2.2 - 2.2i); 225 } 226 for (k = 0; k < size(__CZ__Ck); k++) { 227 print __CZ__Ck[k]; 228 } 229} 230 231/*Kronecker delta function */ 232define kroneckerdelta(i, j) 233{ 234 if (isnull(j)) { 235 if (i == 0) { 236 return 1; 237 } else { 238 return 0; 239 } 240 } 241 if (i != j) { 242 return 0; 243 } else { 244 return 1; 245 } 246} 247 248/* (Pre-)Computed terms of the Stirling series */ 249static __CZ__precomp_stirling; 250/* Number of terms in mat[0,0], precision in mat[0,1] with |z| >= mat[0,2]*/ 251static __CZ__stirling_params; 252 253define __CZ__lngstirling(z, n) 254{ 255 local k head sum z2 bernterm zz; 256 head = (z - 1 / 2) * ln(z) - z + (ln(2 * pi()) / 2); 257 sum = 0; 258 bernterm = 0; 259 zz = z; 260 z2 = z ^ 2; 261 262 if (size(__CZ__precomp_stirling) < n) { 263 __CZ__precomp_stirling = mat[n]; 264 for (k = 1; k <= n; k++) { 265 bernterm = bernoulli(2 * k) / (2 * k * (2 * k - 1)); 266 __CZ__precomp_stirling[k - 1] = bernterm; 267 sum += bernterm / zz; 268 zz *= z2; 269 } 270 } else { 271 for (k = 1; k <= n; k++) { 272 bernterm = __CZ__precomp_stirling[k - 1]; 273 sum += bernterm / zz; 274 zz *= z2; 275 } 276 } 277 return head + sum; 278} 279 280/* Compute error for Stirling series for "z" with "k" terms*/ 281define R(z, k) 282{ 283 local a b ab stirlingterm; 284 stirlingterm = bernoulli(2 * k) / (2 * k * (2 * k - 1) * z ^ (2 * k)); 285 a = abs(stirlingterm); 286 b = (cos(.5 * arg(z) ^ (2 * k))); 287 ab = a / b; 288 /* a/b might round to zero */ 289 if (abs(ab) < epsilon()) { 290 return digits(1 / epsilon()); 291 } 292 return digits(1 / (a / b)); 293} 294 295/*Low precision lngamma(z) for testing the branch correction */ 296define lngammasmall(z) 297{ 298 local ret eps flag corr; 299 flag = 0; 300 if (isint(z)) { 301 if (z <= 0) { 302 return newerror("gamma(z): z is a negative integer"); 303 } else { 304 /* may hold up accuracy a bit longer, but YMMV */ 305 if (z < 20) { 306 return ln((z - 1) !); 307 } else { 308 return __CZ__lngamma_lanczos(z); 309 } 310 } 311 } else { 312 if (re(z) < 0) { 313 if (im(z) && im(z) < 0) { 314 z = conj(z); 315 flag = 1; 316 } 317 318 ret = ln(pi() / sin(pi() * z)) - __CZ__lngamma_lanczos(1 - z); 319 320 if (!im(z)) { 321 if (abs(z) <= 1 / 2) { 322 ret = re(ret) - pi() * 1i; 323 } else if (isodd(floor(abs(re(z))))) { 324 ret = re(ret) + (ceil(abs(z)) * pi() * 1i); 325 } else if (iseven(floor(abs(re(z))))) { 326 /* < n+1/2 */ 327 if (iseven(floor(abs(z)))) { 328 ret = re(ret) + (ceil(abs(z) - 1 / 2 - 1) * pi() * 1i); 329 } else { 330 ret = re(ret) + (ceil(abs(z) - 1 / 2) * pi() * 1i); 331 } 332 } 333 } else { 334 corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4); 335 ret = ret + 2 * (corr * pi()) * 1i; 336 } 337 if (flag == 1) { 338 ret = conj(ret); 339 } 340 return ret; 341 } 342 ret = (__CZ__lngamma_lanczos(z)); 343 return ret; 344 } 345} 346 347 348/* 349 The logarithmic gamma function lngamma(z) 350 It has a different principal branch than log(gamma(z)) 351*/ 352define lngamma(z) 353{ 354 local ret eps flag increasedby Z k tmp tmp2 decdigits; 355 local corr d37 d52 termcount; 356 /* z \in \mathbf{Z} */ 357 if (isint(z)) { 358 /*The gamma-function has poles at zero and the negative integers */ 359 if (z <= 0) { 360 return newerror("lngamma(z): z is a negative integer"); 361 } else { 362 /* may hold up accuracy a bit longer, but YMMV */ 363 if (z < 20) { 364 return ln((z - 1) !); 365 } else { 366 return __CZ__lngammarp(z); 367 } 368 } 369 } else { /*if(isint(z)) */ 370 if (re(z) > 0) { 371 flag = 0; 372 eps = epsilon(); 373 epsilon(eps * 1e-3); 374 375 /* Compute values on the real line with Spouge's algorithm */ 376 if (!im(z)) { 377 ret = __CZ__lngammarp(z); 378 epsilon(eps); 379 return ret; 380 } 381 /* Do the rest with the Stirling series. */ 382 /* This code repeats down under. Booh! */ 383 /* Make it a positive im(z) */ 384 if (im(z) < 0) { 385 z = conj(z); 386 flag = 1; 387 } 388 /* Evaluate the number of terms needed */ 389 decdigits = floor(digits(1 / eps)); 390 /* 20 dec. digits is the default in calc(?) */ 391 if (decdigits <= 21) { 392 /* set 20 as the minimum */ 393 epsilon(1e-22); 394 increasedby = 0; 395 /* inflate z */ 396 Z = z; 397 while (abs(z) < 10) { 398 z++; 399 increasedby++; 400 } 401 402 ret = __CZ__lngstirling(z, 11); 403 /* deflate z */ 404 if (increasedby > 1) { 405 for (k = 0; k < increasedby; k++) { 406 ret = ret - ln(Z + (k)); 407 } 408 } 409 /* Undo conjugate if one has been applied */ 410 if (flag == 1) { 411 ret = conj(ret); 412 } 413 epsilon(eps); 414 return ret; 415 } else { 416 d37 = decdigits * .37; 417 d52 = decdigits * .52; 418 termcount = ceil(d52); 419 if (abs(z) >= d52) { 420 if (abs(im(z)) >= d37) { 421 termcount = ceil(d37); 422 } else { 423 termcount = ceil(d52); 424 } 425 } 426 427 Z = z; 428 increasedby = 0; 429 /* inflate z */ 430 if (abs(im(z)) >= d37) { 431 while (abs(z) < d52 + 1) { 432 z++; 433 increasedby++; 434 } 435 } else { 436 tmp = R(z, termcount); 437 tmp2 = tmp; 438 while (tmp2 < decdigits) { 439 z++; 440 increasedby++; 441 tmp2 = R(z, termcount); 442 if (tmp2 < tmp) { 443 return 444 newerror(strcat 445 ("lngamma(1): something happened ", 446 "that shouldn't have happened")); 447 } 448 } 449 } 450 451 corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4); 452 453 ret = __CZ__lngstirling(z, termcount); 454 455 /* deflate z */ 456 if (increasedby > 1) { 457 for (k = 0; k < increasedby; k++) { 458 ret = ret - ln(Z + (k)); 459 } 460 } 461 /* Undo conjugate if one has been applied */ 462 if (flag == 1) { 463 ret = conj(ret); 464 } 465 epsilon(eps); 466 return ret; 467 } /* if(decdigits <= 20){...} else */ 468 } else { /* re(z)<0 */ 469 eps = epsilon(); 470 epsilon(eps * 1e-3); 471 472 /* Use Spouge's approximation on the real line */ 473 if (!im(z)) { 474 /* reflection */ 475 ret = ln(pi() / sin(pi() * z)) - __CZ__lngammarp(1 - z); 476 /* it is log(gamma) and im(log(even(-x))) = k\pi, therefore: */ 477 if (abs(z) <= 1 / 2) { 478 ret = re(ret) - pi() * 1i; 479 } else if (isodd(floor(abs(re(z))))) { 480 ret = re(ret) + (ceil(abs(z)) * pi() * 1i); 481 } else if (iseven(floor(abs(re(z))))) { 482 /* < n+1/2 */ 483 if (iseven(floor(abs(z)))) { 484 ret = re(ret) + (ceil(abs(z) - 1 / 2 - 1) * pi() * 1i); 485 } else { 486 ret = re(ret) + (ceil(abs(z) - 1 / 2) * pi() * 1i); 487 } 488 } 489 epsilon(eps); 490 return ret; 491 } else { 492 /* Make it a positive im(z) */ 493 if (im(z) < 0) { 494 z = conj(z); 495 flag = 1; 496 } 497 /* Evaluate the number of terms needed */ 498 decdigits = floor(digits(1 / eps)); 499 /* 20 dec. digits is the default in calc(?) */ 500 if (decdigits <= 21) { 501 /* set 20 as the minimum */ 502 epsilon(1e-22); 503 termcount = 11; 504 increasedby = 0; 505 Z = z; 506 /* inflate z */ 507 if (im(z) >= digits(1 / epsilon()) * .37) { 508 while (abs(1 - z) < 10) { 509 z--; 510 increasedby++; 511 } 512 } else { 513 tmp = R(1 - z, termcount); 514 tmp2 = tmp; 515 while (tmp2 < 21) { 516 z--; 517 increasedby++; 518 tmp2 = R(1 - z, termcount); 519 if (tmp2 < tmp) { 520 return 521 newerror(strcat 522 ("lngamma(1): something happened ", 523 "that should not have happened")); 524 } 525 } 526 } 527 528 corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4); 529 530 /* reflection */ 531 ret = 532 ln(pi() / sin(pi() * z)) - __CZ__lngstirling(1 - z, 533 termcount); 534 535 /* deflate z */ 536 if (increasedby > 0) { 537 for (k = 0; k < increasedby; k++) { 538 ret = ret + ln(z + (k)); 539 } 540 } 541 /* Apply correction term */ 542 ret = ret + 2 * (corr * pi()) * 1i; 543 /* Undo conjugate if one has been applied */ 544 if (flag == 1) { 545 ret = conj(ret); 546 } 547 548 epsilon(eps); 549 return ret; 550 } else { 551 d37 = decdigits * .37; 552 d52 = decdigits * .52; 553 termcount = ceil(d52); 554 if (abs(z) >= d52) { 555 if (abs(im(z)) >= d37) { 556 termcount = ceil(d37); 557 } else { 558 termcount = ceil(d52); 559 } 560 } 561 increasedby = 0; 562 Z = z; 563 /* inflate z */ 564 if (abs(im(z)) >= d37) { 565 while (abs(1 - z) < d52 + 1) { 566 z--; 567 increasedby++; 568 } 569 } else { 570 tmp = R(1 - z, termcount); 571 tmp2 = tmp; 572 while (tmp2 < decdigits) { 573 z--; 574 increasedby++; 575 tmp2 = R(1 - z, termcount); 576 if (tmp2 < tmp) { 577 return 578 newerror(strcat 579 ("lngamma(1): something happened ", 580 "that should not have happened")); 581 } 582 } 583 } 584 corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4); 585 /* reflection */ 586 ret = 587 ln(pi() / sin(pi() * (z))) - __CZ__lngstirling(1 - z, 588 termcount); 589 /* deflate z */ 590 if (increasedby > 0) { 591 for (k = 0; k < increasedby; k++) { 592 ret = ret + ln(z + (k)); 593 } 594 } 595 /* Apply correction term */ 596 ret = ret + 2 * ((corr) * pi()) * 1i; 597 /* Undo conjugate if one has been applied */ 598 if (flag == 1) { 599 ret = conj(ret); 600 } 601 epsilon(eps); 602 return ret; 603 } /* if(decdigits <= 20){...} else */ 604 } /*if(!im(z)){...} else */ 605 } /* if(re(z) > 0){...} else */ 606 } /*if(isint(z)){...} else */ 607} 608 609/* Warning about large values? */ 610define gamma(z) 611{ 612 613 /* exp(log(gamma(z))) = exp(lngamma(z)), so use Spouge here? */ 614 local ret eps; 615 if (isint(z)) { 616 if (z <= 0) { 617 return newerror("gamma(z): z is a negative integer"); 618 } else { 619 /* may hold up accuracy a bit longer, but YMMV */ 620 if (z < 20) { 621 return (z - 1) ! *1.0; 622 } else { 623 eps = epsilon(epsilon() * 1e-2); 624 ret = exp(lngamma(z)); 625 epsilon(eps); 626 return ret; 627 } 628 } 629 } else { 630 eps = epsilon(epsilon() * 1e-2); 631 ret = exp(lngamma(z)); 632 epsilon(eps); 633 return ret; 634 } 635} 636 637define __CZ__harmonicF(a, b, s) 638{ 639 local c; 640 if (b == a) { 641 return s; 642 } 643 if (b - a > 1) { 644 c = (b + a) >> 1; 645 return (__CZ__harmonicF(a, c, 1 / a) + 646 __CZ__harmonicF(c + 1, b, 1 / b)); 647 } 648 return (1 / a + 1 / b); 649} 650 651define harmonic(limit) 652{ 653 if (!isint(limit)) { 654 return newerror("harmonic(limit): limit is not an integer"); 655 } 656 if (limit <= 0) { 657 return newerror("harmonic(limit): limit is <=0"); 658 } 659 /* The binary splitting algorithm returns 0 for f(1,1,0) */ 660 if (limit == 1) { 661 return 1; 662 } 663 return __CZ__harmonicF(1, limit, 0); 664} 665 666/* lower incomplete gamma function */ 667 668/* lower 669 for z <= 1.1 670 */ 671define __CZ__gammainc_series_lower(a, z) 672{ 673 local k ret tmp eps fact; 674 ret = 0; 675 k = 0; 676 tmp = 1; 677 fact = 1; 678 eps = epsilon(); 679 while (abs(tmp - ret) > eps) { 680 tmp = ret; 681 ret += (z ^ (k + a)) / ((a + k) * fact); 682 k++; 683 fact *= -k; 684 } 685 return gamma(a) - ret; 686} 687 688/* lower 689 for z > 1.1 690 */ 691define __CZ__gammainc_cf(a, z, n) 692{ 693 local ret k num1 denom1 num2 denom2; 694 ret = 0; 695 for (k = n + 1; k > 1; k--) { 696 ret = ((1 - k) * (k - 1 - a)) / (2 * k - 1 + z - a + ret); 697 } 698 return ((z ^ a * exp(-z)) / (1 + z - a + ret)); 699} 700 701/* G(n,z) lower*/ 702define __CZ__gammainc_integer_a(a, z) 703{ 704 local k sum fact zz; 705 for (k = 0; k < a; k++) { 706 sum += (z ^ k) / (k !); 707 } 708 return (a - 1) ! *exp(-z) * sum; 709} 710 711/* 712 P = precision in dec digits 713 n = 1,2,3... 714 a,z => G(a,z) 715*/ 716define __CZ__endcf(n, a, z, P) 717{ 718 local ret; 719 720 ret = 721 P * ln(10) + ln(4 * pi() * sqrt(n)) + re(z + (3 / 2 - a) * ln(z) - 722 lngamma(1 - a)); 723 ret = ret / (sqrt(8 * (abs(z) + re(z)))); 724 return ret ^ 2; 725 726} 727 728/* lower incomplete gamma function */ 729define gammainc(a, z) 730{ 731 local ret nterms eps epsilon tmp_before tmp_after n A B C sum k; 732 if (z == 0) { 733 return 1; 734 } 735 if (isint(a)) { 736 if (a > 0) { 737 if (a == 1) { 738 return exp(-z); 739 } 740 return __CZ__gammainc_integer_a(a, z); 741 } else { 742 if (a == 0) { 743 return -expoint(-z) + 1 / 2 * (ln(-z) - ln(-1 / z)) - ln(z); 744 } else if (a == -1) { 745 return expoint(-z) + 1 / 2 * (ln(-1 / z) - ln(-z)) + ln(z) + 746 (exp(-z) / z); 747 } else { 748 A = (-1) ^ ((-a) - 1) / ((-a) !); 749 B = expoint(-z) - 1 / 2 * (ln(-z) - ln(1 / -z)) + ln(z); 750 C = exp(-z); 751 sum = 0; 752 for (k = 1; k < -a; k++) { 753 sum += (z ^ (k - a - 1)) / (fallingfactorial(-a, k)); 754 } 755 return A * B - C * sum; 756 } 757 } 758 } 759 if (re(z) <= 1.1 || re(z) < a + 1) { 760 eps = epsilon(epsilon() * 1e-10); 761 ret = __CZ__gammainc_series_lower(a, z); 762 epsilon(eps); 763 return ret; 764 } else { 765 eps = epsilon(epsilon() * 1e-10); 766 if (abs(exp(-z)) <= eps) { 767 return 0; 768 } 769 tmp_before = 0; 770 tmp_after = 1; 771 n = 1; 772 while (ceil(tmp_before) != ceil(tmp_after)) { 773 tmp_before = tmp_after; 774 tmp_after = __CZ__endcf(n++, a, z, digits(1 / eps)); 775 /* a still quite arbitrary limit */ 776 if (n > 10) { 777 return 778 newerror(strcat 779 ("gammainc: evaluating limit for continued ", 780 "fraction does not converge")); 781 } 782 } 783 ret = __CZ__gammainc_cf(a, z, ceil(tmp_after)); 784 epsilon(eps); 785 return ret; 786 } 787} 788 789 790define heavisidestep(x) 791{ 792 return (1 + sign(x)) / 2; 793} 794 795define NUMBER_POSITIVE_INFINITY() 796{ 797 return 1 / epsilon(); 798} 799 800define NUMBER_NEGATIVE_INFINITY() 801{ 802 return -(1 / epsilon()); 803} 804 805static TRUE = 1; 806static FALSE = 0; 807 808define g(prec) 809{ 810 local eps ret; 811 if (!isnull(prec)) { 812 eps = epsilon(prec); 813 ret = -psi(1); 814 epsilon(eps); 815 return ret; 816 } 817 return -psi(1); 818} 819 820define __CZ__series_converged(new, old, max) 821{ 822 local eps; 823 if (isnull(max)) { 824 eps = epsilon(); 825 } else { 826 eps = max; 827 } 828 if (abs(re(new - old)) <= eps * abs(re(new)) 829 && abs(im(new - old)) <= eps * abs(im(new))) { 830 return TRUE; 831 } 832 return FALSE; 833} 834 835define __CZ__ei_power(z) 836{ 837 local tmp ei k old; 838 ei = g() + ln(abs(z)) + sgn(im(z)) * 1i * abs(arg(z));; 839 tmp = 1; 840 k = 1; 841 while (k) { 842 tmp *= z / k; 843 old = ei; 844 ei += tmp / k; 845 if (__CZ__series_converged(ei, old)) { 846 break; 847 } 848 k++; 849 } 850 return ei; 851} 852 853define __CZ__ei_asymp(z) 854{ 855 local ei old tmp k; 856 ei = sgn(im(z)) * 1i * pi(); 857 tmp = exp(z) / z; 858 for (k = 1; k <= floor(abs(z)) + 1; k++) { 859 old = ei; 860 ei += tmp; 861 if (__CZ__series_converged(ei, old)) { 862 return ei; 863 } 864 tmp *= k / z; 865 } 866 return newerror("expoint: asymptotic series does not converge"); 867} 868 869define __CZ__ei_cf(z) 870{ 871 local ei c d k old; 872 ei = sgn(im(z)) * 1i * pi(); 873 if (ei != 0) { 874 c = 1 / ei; 875 d = 0; 876 c = 1 / (1 - z - exp(z) * c); 877 d = 1 / (1 - z - exp(z) * d); 878 ei *= d / c; 879 } else { 880 c = NUMBER_POSITIVE_INFINITY(); 881 d = 0; 882 c = 0; 883 d = 1 / (1 - z - exp(z) * d); 884 ei = d * (-exp(z)); 885 } 886 k = 1; 887 while (1) { 888 c = 1 / (2 * k + 1 - z - k * k * c); 889 d = 1 / (2 * k + 1 - z - k * k * d); 890 old = ei; 891 ei *= d / c; 892 if (__CZ__series_converged(ei, old)) { 893 break; 894 } 895 k++; 896 } 897 return ei; 898} 899 900define expoint(z) 901{ 902 local ei eps ret; 903 eps = epsilon(epsilon() * 1e-5); 904 if (abs(z) >= NUMBER_POSITIVE_INFINITY()) { 905 if (config("user_debug") > 0) { 906 print "expoint: abs(z) > +inf"; 907 } 908 ret = sgn(im(z)) * 1i * pi() + exp(z) / z; 909 epsilon(eps); 910 return ret; 911 } 912 if (abs(z) > 2 - 1.035 * log(epsilon())) { 913 if (config("user_debug") > 0) { 914 print "expoint: using asymptotic series"; 915 } 916 ei = __CZ__ei_asymp(z); 917 if (!iserror(ei)) { 918 ret = ei; 919 epsilon(eps); 920 return ret; 921 } 922 } 923 if (abs(z) > 1 && (re(z) < 0 || abs(im(z)) > 1)) { 924 if (config("user_debug") > 0) { 925 print "expoint: using continued fraction"; 926 } 927 ret = __CZ__ei_cf(z); 928 epsilon(eps); 929 return ret; 930 } 931 if (abs(z) > 0) { 932 if (config("user_debug") > 0) { 933 print "expoint: using power series"; 934 } 935 ret = __CZ__ei_power(z); 936 epsilon(eps); 937 return ret; 938 } 939 if (abs(z) == 0) { 940 if (config("user_debug") > 0) { 941 print "expoint: abs(z) = zero "; 942 } 943 epsilon(eps); 944 return NUMBER_NEGATIVE_INFINITY(); 945 } 946} 947 948define erf(z) 949{ 950 return sqrt(z ^ 2) / z * (1 - 1 / sqrt(pi()) * gammainc(1 / 2, z ^ 2)); 951} 952 953define erfc(z) 954{ 955 return 1 - erf(z); 956} 957 958define erfi(z) 959{ 960 return -1i * erf(1i * z); 961} 962 963define faddeeva(z) 964{ 965 return exp(-z ^ 2) * erfc(-1i * z); 966} 967 968define gammap(a, z) 969{ 970 return gammainc(a, z) / gamma(a); 971} 972 973define gammaq(a, z) 974{ 975 return 1 - gammap(a, z); 976} 977 978define lnbeta(a, b) 979{ 980 local ret eps; 981 eps = epsilon(epsilon() * 1e-3); 982 ret = (lngamma(a) + lngamma(b)) - lngamma(a + b); 983 epsilon(eps); 984 return ret; 985} 986 987define beta(a, b) 988{ 989 return exp(lnbeta(a, b)); 990} 991 992 993define __CZ__ibetacf_a(a, b, z, n) 994{ 995 local A B m places; 996 if (n == 1) { 997 return 1; 998 } 999 m = n - 1; 1000 places = highbit(1 + int (1 / epsilon())) +1; 1001 A = bround((a + m - 1) * (a + b + m - 1) * m * (b - m) * z ^ 2, places++); 1002 B = bround((a + 2 * (m) - 1) ^ 2, places++); 1003 return A / B; 1004} 1005 1006define __CZ__ibetacf_b(a, b, z, n) 1007{ 1008 local A B m places; 1009 places = highbit(1 + int (1 / epsilon())) +1; 1010 m = n - 1; 1011 A = bround((m * (b - m) * z) / (a + 2 * m - 1), places++); 1012 B = bround(((a + m) * (a - (a + b) * z + 1 + m * (2 - z))) / (a + 2 * m + 1013 1), places++); 1014 return m + A + B; 1015} 1016 1017/* Didonato-Morris */ 1018define __CZ__ibeta_cf_var_dm(a, b, z, max) 1019{ 1020 local m f c d check del h qab qam qap eps places; 1021 1022 eps = epsilon(); 1023 1024 if (isnull(max)) { 1025 max = 100; 1026 } 1027 places = highbit(1 + int (1 / epsilon())) + 1; 1028 f = eps; 1029 c = f; 1030 d = 0; 1031 for (m = 1; m <= max; m++) { 1032 d = bround(__CZ__ibetacf_b(a, b, z, m) + 1033 __CZ__ibetacf_a(a, b, z, m) * d, places++); 1034 if (abs(d) < eps) { 1035 d = eps; 1036 } 1037 c = bround(__CZ__ibetacf_b(a, b, z, m) + 1038 __CZ__ibetacf_a(a, b, z, m) / c, places++); 1039 if (abs(c) < eps) { 1040 c = eps; 1041 } 1042 d = 1 / d; 1043 check = c * d; 1044 f = f * check; 1045 if (abs(check - 1) < eps) { 1046 break; 1047 } 1048 } 1049 if (m > max) { 1050 return newerror("ibeta: continuous fraction does not converge"); 1051 } 1052 return f; 1053} 1054 1055define betainc_complex(z, a, b) 1056{ 1057 local factor ret eps cf sum k N places tmp tmp2; 1058 1059 if (z == 0) { 1060 if (re(a) > 0) { 1061 return 0; 1062 } 1063 if (re(a) < 0) { 1064 return newerror("betainc_complex: z == 0 and re(a) < 0"); 1065 } 1066 } 1067 if (z == 1) { 1068 if (re(b) > 0) { 1069 return 1; 1070 } else { 1071 return newerror("betainc_complex: z == 1 and re(b) < 0"); 1072 } 1073 } 1074 if (b <= 0) { 1075 if (isint(b)) { 1076 return 0; 1077 } else { 1078 return newerror("betainc_complex: b <= 0"); 1079 } 1080 } 1081 if (z == 1 / 2 && (a == b)) { 1082 return 1 / 2; 1083 } 1084 if (isint(a) && isint(b)) { 1085 eps = epsilon(epsilon() * 1e-10); 1086 N = a + b - 1; 1087 sum = 0; 1088 for (k = a; k <= N; k++) { 1089 tmp = ln(z) * k + ln(1 - z) * (N - k); 1090 tmp2 = exp(ln(comb(N, k)) + tmp); 1091 sum += tmp2; 1092 } 1093 epsilon(eps); 1094 return sum; 1095 } else if (re(z) <= re((a + 1) / (a + b + 2))) { 1096 eps = epsilon(epsilon() * 1e-10); 1097 places = highbit(1 + int (1 / epsilon())) + 1; 1098 factor = bround((ln(z ^ a * (1 - z) ^ b) - lnbeta(a, b)), places); 1099 cf = bround(__CZ__ibeta_cf_var_dm(a, b, z), places); 1100 ret = factor + ln(cf); 1101 if (abs(ret//ln(2)) >= places) { 1102 ret = 0; 1103 } else { 1104 ret = bround(exp(factor + ln(cf)), places); 1105 } 1106 epsilon(eps); 1107 return ret; 1108 } else if (re(z) > re((a + 1) / (a + b + 2)) 1109 || re(1 - z) < re((b + 1) / (a + b + 2))) { 1110 ret = 1 - betainc_complex(1 - z, b, a); 1111 } 1112 return ret; 1113} 1114 1115 1116 1117 1118 1119/******************************************************************************/ 1120/* 1121 Purpose: 1122 1123 __CZ__ibetaas63 computes the incomplete Beta function ratio. 1124 1125 Licensing: 1126 1127 This code is distributed under the GNU LGPL license. 1128 1129 Modified: 1130 1131 2013-08-03 20:52:05 +0000 1132 1133 Author: 1134 1135 Original FORTRAN77 version by KL Majumder, GP Bhattacharjee. 1136 C version by John Burkardt 1137 Calc version by Christoph Zurnieden 1138 1139 Reference: 1140 1141 KL Majumder, GP Bhattacharjee, 1142 Algorithm AS 63: 1143 The incomplete Beta Integral, 1144 Applied Statistics, 1145 Volume 22, Number 3, 1973, pages 409-411. 1146 1147 Parameters: 1148 1149 Input, x, the argument, between 0 and 1. 1150 1151 Input, a, b, the parameters, which 1152 must be positive. 1153 1154 1155 Output, the value of the incomplete 1156 Beta function ratio. 1157*/ 1158define __CZ__ibetaas63(x, a, b, beta) 1159{ 1160 local ai betain cx indx ns aa asb bb rx temp term value xx acu places; 1161 acu = epsilon(); 1162 1163 value = x; 1164 /* inverse incbeta calculates it already */ 1165 if (isnull(beta)) 1166 beta = lnbeta(a, b); 1167 1168 if (a <= 0.0 || b <= 0.0) { 1169 return newerror("betainc: domain error: a < 0 and/or b < 0"); 1170 } 1171 if (x < 0.0 || 1.0 < x) { 1172 return newerror("betainc: domain error: x<0 or x>1"); 1173 } 1174 if (x == 0.0 || x == 1.0) { 1175 return value; 1176 } 1177 asb = a + b; 1178 cx = 1.0 - x; 1179 1180 if (a < asb * x) { 1181 xx = cx; 1182 cx = x; 1183 aa = b; 1184 bb = a; 1185 indx = 1; 1186 } else { 1187 xx = x; 1188 aa = a; 1189 bb = b; 1190 indx = 0; 1191 } 1192 1193 term = 1.0; 1194 ai = 1.0; 1195 value = 1.0; 1196 ns = floor(bb + cx * asb); 1197 1198 rx = xx / cx; 1199 temp = bb - ai; 1200 if (ns == 0) { 1201 rx = xx; 1202 } 1203 places = highbit(1 + int (1 / acu)) + 1; 1204 while (1) { 1205 term = bround(term * temp * rx / (aa + ai), places++); 1206 value = value + term;; 1207 temp = abs(term); 1208 1209 if (temp <= acu && temp <= abs(acu * value)) { 1210 value = value * exp(aa * ln(xx) 1211 + (bb - 1.0) * ln(cx) - beta) / aa; 1212 1213 if (indx) { 1214 value = 1.0 - value; 1215 } 1216 break; 1217 } 1218 1219 ai = ai + 1.0; 1220 ns = ns - 1; 1221 1222 if (0 <= ns) { 1223 temp = bb - ai; 1224 if (ns == 0) { 1225 rx = xx; 1226 } 1227 } else { 1228 temp = asb; 1229 asb = asb + 1.0; 1230 } 1231 } 1232 epsilon(acu); 1233 return value; 1234} 1235 1236/* 1237 z 1238 / 1239 [ b - 1 a - 1 1240 1/beta(a,b) * I (1 - t) t dt 1241 ] 1242 / 1243 0 1244 1245*/ 1246 1247define betainc(z, a, b) 1248{ 1249 local factor ret eps cf sum k N places tmp tmp2; 1250 1251 if (im(z) || im(a) || im(b)) 1252 return betainc_complex(z, a, b); 1253 1254 if (z == 0) { 1255 if (re(a) > 0) { 1256 return 0; 1257 } 1258 if (re(a) < 0) { 1259 return newerror("betainc: z == 0 and re(a) < 0"); 1260 } 1261 } 1262 if (z == 1) { 1263 if (re(b) > 0) { 1264 return 1; 1265 } else { 1266 return newerror("betainc: z == 1 and re(b) < 0"); 1267 } 1268 } 1269 if (b <= 0) { 1270 if (isint(b)) { 1271 return 0; 1272 } else { 1273 return newerror("betainc: b <= 0"); 1274 } 1275 } 1276 if (z == 1 / 2 && a == b) { 1277 return 1 / 2; 1278 } 1279 return __CZ__ibetaas63(z, a, b); 1280 1281} 1282 1283define __CZ__erfinvapprox(x) 1284{ 1285 local a; 1286 a = 0.147; 1287 return sgn(x) * 1288 sqrt(sqrt 1289 ((2 / (pi() * a) + (ln(1 - x ^ 2)) / 2) ^ 2 - (ln(1 - x ^ 2)) / a) 1290 - (2 / (pi() * a) + (ln(1 - x ^ 2)) / 2)); 1291} 1292 1293/* complementary inverse error function, faster at about x < 1-.91 1294 Henry E. Fettis. "A stable algorithm for computing the inverse error function 1295 in the 'tail-end' region" Math. Comp., 28:585-587, 1974. 1296*/ 1297define __CZ__inverffettis(x, n) 1298{ 1299 local y sqrtpi oldy k places; 1300 if (isnull(n)) 1301 n = 205; 1302 y = erfinvapprox(1 - x); 1303 places = highbit(1 + int (1 / epsilon())) + 1; 1304 sqrtpi = sqrt(pi()); 1305 do { 1306 oldy = y; 1307 k++; 1308 y = bround((ln(__CZ__fettiscf(y, n) / (sqrtpi * x))) ^ (1 / 2), places); 1309 } while (abs(y - oldy) / y > epsilon()); 1310 return y; 1311} 1312 1313/* cf for erfc() */ 1314define __CZ__fettiscf(y, n) 1315{ 1316 local k t tt r a b; 1317 t = 1 / y; 1318 tt = t ^ 2 / 2; 1319 for (k = n; k > 0; k--) { 1320 a = 1; 1321 b = k * tt; 1322 r = b / (a + r); 1323 } 1324 return t / (1 + r); 1325} 1326 1327/* inverse error function, faster at about x<=.91*/ 1328define __CZ__inverfbin(x) 1329{ 1330 local places approx flow fhigh eps high low mid fmid epsilon; 1331 approx = erfinvapprox(x); 1332 epsilon = epsilon(); 1333 high = approx + 1e-4; 1334 low = -1; 1335 places = highbit(1 + int (1 / epsilon)) +1; 1336 fhigh = x - erf(high); 1337 flow = x - erf(low); 1338 while (1) { 1339 mid = bround(high - fhigh * (high - low) / (fhigh - flow), places); 1340 if ((mid == low) || (mid == high)) { 1341 places++; 1342 } 1343 fmid = x - erf(mid); 1344 if (abs(fmid) < epsilon) { 1345 return mid; 1346 } 1347 if (sgn(fmid) == sgn(flow)) { 1348 low = mid; 1349 flow = fmid; 1350 } else { 1351 high = mid; 1352 fhigh = fmid; 1353 } 1354 } 1355} 1356 1357define erfinv(x) 1358{ 1359 local ret approx a eps y old places errfunc sqrtpihalf flag k; 1360 if (x < -1 || x > 1) 1361 return newerror("erfinv: input out of domain (-1<=x<=1)"); 1362 if (x == 0) 1363 return 0; 1364 if (x == -1) 1365 return NUMBER_NEGATIVE_INFINITY(); 1366 if (x == +1) 1367 return NUMBER_POSITIVE_INFINITY(); 1368 1369 if (x < 0) { 1370 x = -x; 1371 flag = 1; 1372 } 1373 /* No need for full precision */ 1374 eps = epsilon(1e-20); 1375 if (eps >= 1e-40) { 1376 /* Winitzki, Sergei (6 February 2008). "A handy approximation for the 1377 * error function and its inverse" */ 1378 a = 0.147; 1379 y = sgn(x) * sqrt(sqrt((2 / (pi() * a) 1380 + (ln(1 - x ^ 2)) / 2) ^ 2 1381 - (ln(1 - x ^ 2)) / a) 1382 - (2 / (pi() * a) + (ln(1 - x ^ 2)) / 2)); 1383 1384 } else { 1385 /* 20 digits instead of 5 */ 1386 if (x <= .91) { 1387 y = __CZ__inverfbin(x); 1388 } else { 1389 y = __CZ__inverffettis(1 - x); 1390 } 1391 1392 if (eps <= 1e-20) { 1393 epsilon(eps); 1394 return y; 1395 } 1396 } 1397 epsilon(eps); 1398 /* binary digits in number (here: number = epsilon()) */ 1399 places = highbit(1 + int (1 / eps)) +1; 1400 sqrtpihalf = 2 / sqrt(pi()); 1401 k = 0; 1402 /* 1403 * Do some Newton-Raphson steps to reach final accuracy. 1404 * Only a couple of steps are necessary but calculating the error function 1405 * at higher precision is quite costly; 1406 */ 1407 do { 1408 old = y; 1409 errfunc = bround(erf(y), places); 1410 if (abs(errfunc - x) <= eps) { 1411 break; 1412 } 1413 y = bround(y - (errfunc - x) / (sqrtpihalf * exp(-y ^ 2)), places); 1414 k++; 1415 } while (1); 1416 /* 1417 * This is not really necessary but e.g: 1418 * ; epsilon(1e-50) 1419 * 0.00000000000000000000000000000000000000000000000001 1420 * ; erfinv(.9999999999999999999999999999999) 1421 * 8.28769266865549025938 1422 * ; erfinv(.999999999999999999999999999999) 1423 * 8.14861622316986460738453487549552168842204512959346 1424 * ; erf(8.28769266865549025938) 1425 * 0.99999999999999999999999999999990000000000000000000 1426 * ; erf(8.14861622316986460738453487549552168842204512959346) 1427 * 0.99999999999999999999999999999900000000000000000000 1428 * The precision "looks too short". 1429 */ 1430 if (k == 0) { 1431 y = bround(y - (errfunc - x) / (sqrtpihalf * exp(-y ^ 2)), places); 1432 } 1433 if (flag == 1) { 1434 y = -y; 1435 } 1436 return y; 1437} 1438 1439 1440/* 1441 * restore internal function from resource debugging 1442 */ 1443config("resource_debug", resource_debug_level),; 1444if (config("resource_debug") & 3) { 1445 print "zeta(z)"; 1446 print "psi(z)"; 1447 print "polygamma(m,z)"; 1448 print "lngamma(z)"; 1449 print "gamma(z)"; 1450 print "harmonic(limit)"; 1451 print "gammainc(a,z)"; 1452 print "heavisidestep(x)"; 1453 print "expoint(z)"; 1454 print "erf(z)"; 1455 print "erfinv(x)"; 1456 print "erfc(z)"; 1457 print "erfi(z)"; 1458 print "erfinv(x)"; 1459 print "faddeeva(z)"; 1460 print "gammap(a,z)"; 1461 print "gammaq(a,z)"; 1462 print "beta(a,b)"; 1463 print "lnbeta(a,b)"; 1464 print "betainc(z,a,b)"; 1465} 1466