tests/rtest_gamma.mac

   1 /******************************************************************************
   2   rtest_gamma.mac
   3   Test for Factorial, Gamma function and related functions ...
   4 ******************************************************************************/
   5
   6 kill(all);
   7 done;
   8
   9 (oldfpprec:fpprec, fpprec:16, done);
  10 done;
  11
  12 /* Two definitions for numerical test functions
  13    For big results relerror is used */
  14
  15 (closeto(value,compare,tol):=
  16   block(
  17     [abse],
  18     abse:abs(value-compare),if(abse<tol) then true else abse),
  19     done);
  20 done;
  21
  22 (relerror(value,compare,tol):=
  23   block(
  24     [abse],
  25     abse:abs((value-compare)/compare),
  26     if(abse<tol) then true else abse),
  27     done);
  28 done;
  29
  30 /******************************************************************************
  31   Factorial
  32 ******************************************************************************/
  33
  34 /* Factorial has mirror symmetry */
  35
  36 declare(z,complex);
  37 done;
  38
  39 conjugate(factorial(z));
  40 factorial(conjugate(z));
  41
  42 conjugate(factorial(x+%i*y));
  43 factorial(x-%i*y);
  44
  45 /* some small positive integers or the real representation */
  46
  47 fpprec:16;
  48 16;
  49
  50 map(factorial, [0,1,2,3,4]);
  51 [1,1,2,6,24];
  52
  53 closeto(factorial(0.0),1.0,1e-13);
  54 true;
  55 closeto(factorial(1.0),1.0,1e-13);
  56 true;
  57 closeto(factorial(2.0),2.0,1e-13);
  58 true;
  59 closeto(factorial(3.0),6.0,1e-13);
  60 true;
  61 closeto(factorial(4.0),24.0,1e-13);
  62 true;
  63
  64 closeto(factorial(0.0b0),1.0b0,1e-13);
  65 true;
  66 closeto(factorial(1.0b0),1.0b0,1e-13);
  67 true;
  68 closeto(factorial(2.0b0),2.0b0,1e-13);
  69 true;
  70 closeto(factorial(3.0b0),6.0b0,1e-13);
  71 true;
  72 closeto(factorial(4.0b0),24.0b0,1e-13);
  73 true;
  74
  75 /* negative integers or there real representation */
  76
  77 errcatch(factorial(-1));
  78 [];
  79
  80 errcatch(factorial(-1.0));
  81 [];
  82
  83 errcatch(factorial(-1.0b0));
  84 [];
  85
  86 errcatch(factorial(-10));
  87 [];
  88
  89 errcatch(factorial(-10.0));
  90 [];
  91
  92 errcatch(factorial(-10.0b0));
  93 [];
  94
  95 /* half integral values */
  96
  97 factorial(1/2);
  98 sqrt(%pi)/2;
  99 factorial(-1/2);
 100 sqrt(%pi);
 101 factorial(3/2);
 102 3*sqrt(%pi)/4;
 103 factorial(-3/2);
 104 -2*sqrt(%pi);
 105 factorial(5/2);
 106 15*sqrt(%pi)/8;
 107 factorial(-5/2);
 108 4*sqrt(%pi)/3;
 109
 110 /* Expansion for factorial(z+n) and integer n */
 111
 112 factorial_expand:true;
 113 true;
 114
 115 factorial(z+1);
 116 (z+1)*factorial(z);
 117
 118 factorial(z+2);
 119 (z+1)*(z+2)*factorial(z);
 120
 121 factorial(z+3);
 122 (z+1)*(z+2)*(z+3)*factorial(z);
 123
 124 factorial(z-1);
 125 factorial(z)/z;
 126
 127 factorial(z-1);
 128 factorial(z)/z;
 129
 130 factorial(z-2);
 131 factorial(z)/(z*(z-1));
 132
 133 factorial(z-3);
 134 factorial(z)/(z*(z-1)*(z-2));
 135
 136 /* Nested factorials simplifies too, see SF[1486452] */
 137
 138 factorial(factorial(n)/factorial(n-1));
 139 factorial(n);
 140
 141 factorial(sin(factorial(n)/factorial(n-1)));
 142 factorial(sin(n));
 143
 144 factorial_expand:false;
 145 false;
 146
 147 /* minfactorial does not do this job */
 148
 149 minfactorial(factorial(factorial(n)/factorial(n-1)));
 150 factorial(factorial(n)/factorial(n-1));
 151
 152 /* factcomb is the inverse operation to minfactorial
 153    factorial_expand has to be false
 154 */
 155
 156 factcomb((n+1)*(n+2)*(n+3)*n!);
 157 (n+3)!;
 158
 159 factcomb(n!/(n*(n-1)*(n-2)));
 160 (n-3)!;
 161
 162 /* No simplifcation for infinities and undeterminates
 163    with the exception of inf! -> inf */
 164
 165 map(factorial, [inf,minf,infinity,und,ind]);
 166 [inf,factorial(minf),factorial(infinity),factorial(und),factorial(ind)];
 167
 168 /* factlim is set to the value 100,000. This should work.  */
 169
 170 fpprec:64;
 171 64;
 172
 173 relerror(
 174   factorial(factlim),
 175   factorial(bfloat(factlim)),
 176   1b-58);                      /* We loose a lost of digits in relerror */
 177 true;
 178
 179 factorial(factlim+1);
 180 factorial(100000+1);
 181
 182 /* Some real values in double float and bigfloat precision */
 183
 184 fpprec:64;
 185 64;
 186
 187 closeto(
 188   factorial(1.3),
 189   1.166711905198160345041881441202917938533994349719468893970206664b0,
 190   1e-14);
 191 true;
 192
 193 closeto(
 194   factorial(2.3),
 195   2.683437381955768793596327314766711258628187004354778456131475327b0,
 196   1e-14);
 197 true;
 198
 199 closeto(
 200   factorial(3.3),
 201   8.855343360454037018867880138730147153473017114370768905233868579b0,
 202   1e-14);
 203 true;
 204
 205 closeto(
 206   factorial(1.3b0),
 207   1.166711905198160345041881441202917938533994349719468893970206664b0,
 208   1e-62);
 209 true;
 210
 211 closeto(
 212   factorial(2.3b0),
 213   2.683437381955768793596327314766711258628187004354778456131475327b0,
 214   1e-62);
 215 true;
 216
 217 closeto(
 218   factorial(3.3b0),
 219   8.855343360454037018867880138730147153473017114370768905233868579b0,
 220   1e-61);
 221 true;
 222
 223 /* some complex values in double float and bigfloat precision */
 224
 225 closeto(
 226   factorial(1.3+%i),
 227   (0.7191409365372817791473038599462048083254863806205029128993808432b0
 228   +0.5406144679098492753783510221774150545811250310680842509749769021b0*%i),
 229   1e-14);
 230 true;
 231
 232 closeto(
 233   factorial(2.3+%i),
 234   (1.113409686125898816660447855698856004567493644359072448693599037b0
 235   +1.962554212729935112517511210954259433862073952077096690141827718b0*%i),
 236   1e-14);
 237 true;
 238
 239 closeto(
 240   factorial(3.3+%i),
 241   (1.711697751485530982461966712851965381210655074307842390547049105b0
 242   +7.589838588134684687968234851847912136312337686213491526161630507b0*%i),
 243   5.0e-14);
 244 true;
 245
 246 closeto(
 247   factorial(1.3b0+%i),
 248   (0.7191409365372817791473038599462048083254863806205029128993808432b0
 249   +0.5406144679098492753783510221774150545811250310680842509749769021b0*%i),
 250   1e-62);
 251 true;
 252
 253 closeto(
 254   factorial(2.3b0+%i),
 255   (1.113409686125898816660447855698856004567493644359072448693599037b0
 256   +1.962554212729935112517511210954259433862073952077096690141827718b0*%i),
 257   1e-62);
 258 true;
 259
 260 closeto(
 261   factorial(3.3b0+%i),
 262   (1.711697751485530982461966712851965381210655074307842390547049105b0
 263   +7.589838588134684687968234851847912136312337686213491526161630507b0*%i),
 264   1e-61);
 265 true;
 266
 267 /******************************************************************************
 268   General factorial: Tests for genfact(x,y,z)
 269 ******************************************************************************/
 270
 271 genfact(0,0,1);
 272 1;
 273 genfact(1,1,1);
 274 1;
 275 genfact(2,2,1);
 276 2;
 277 genfact(3,3,1);
 278 6;
 279 genfact(4,4,1);
 280 24;
 281
 282 genfact(0,0/2,2);
 283 1;
 284 genfact(1,1/2,2);
 285 1;
 286 genfact(2,2/2,2);
 287 2;
 288 genfact(3,3/2,2);
 289 3;
 290 genfact(4,4/2,2);
 291 8;
 292
 293 genfact(10,10,1);
 294 3628800;
 295 genfact(10,9,1);
 296 3628800;
 297 genfact(10,8,1);
 298 1814400;
 299 genfact(10,7,1);
 300 604800;
 301 genfact(10,6,1);
 302 151200;
 303 genfact(10,5,1);
 304 30240;
 305 genfact(10,4,1);
 306 5040;
 307 genfact(10,3,1);
 308 720;
 309 genfact(10,2,1);
 310 90;
 311 genfact(10,1,1);
 312 10;
 313 genfact(10,0,1);
 314 1;
 315
 316 genfact(10,5,2);
 317 3840;
 318 genfact(10,4,2);
 319 1920;
 320 genfact(10,3,2);
 321 480;
 322 genfact(10,2,2);
 323 80;
 324 genfact(10,1,2);
 325 10;
 326 genfact(10,0,2);
 327 1;
 328
 329 /* for non valid integers we get an error */
 330 errcatch(genfact(-2,-2,1));
 331 [];
 332 errcatch(genfact(2,5,2));
 333 [];
 334
 335 /* for all other numbers we get a noun form */
 336 (3/2)!!;
 337 genfact(3/2,0,2);
 338 (2.5)!!;
 339 genfact(2.5,1,2);
 340 (2.5b0)!!;
 341 genfact(2.5b0,1,2);
 342
 343 /******************************************************************************
 344   Double factorial
 345 ******************************************************************************/
 346
 347 /* Double factorial has mirror symmetry */
 348
 349 declare(z,complex);
 350 done;
 351
 352 conjugate(double_factorial(z));
 353 double_factorial(conjugate(z));
 354
 355 conjugate(double_factorial(x+%i*y));
 356 double_factorial(x-%i*y);
 357
 358 /* No simplifcation for infinities and undeterminates
 359    with the exception of inf: inf! -> inf */
 360
 361 map(factorial, [inf,minf,infinity,und,ind]);
 362 [inf,factorial(minf),factorial(infinity),factorial(und),factorial(ind)];
 363
 364 /* Test the expansion of Double factorial */
 365
 366 double_factorial(n+1);
 367 double_factorial(n+1);
 368
 369 double_factorial(n+2),factorial_expand:true;
 370 (n+2)*double_factorial(n);
 371
 372 double_factorial(n+3),factorial_expand:true;
 373 double_factorial(n+3);
 374
 375 double_factorial(n+4),factorial_expand:true;
 376 (n+2)*(n+4)*double_factorial(n);
 377
 378 double_factorial(n-1),factorial_expand:true;
 379 factorial(n)/double_factorial(n);
 380
 381 double_factorial(n-2),factorial_expand:true;
 382 double_factorial(n)/n;
 383
 384 double_factorial(n-3),factorial_expand:true;
 385 double_factorial(n-3);
 386
 387 double_factorial(n-4),factorial_expand:true;
 388 double_factorial(n)/(n*(n-2));
 389
 390 makegamma(double_factorial(n));
 391 (gamma(n/2+1)*2^((1-cos(%pi*n))/4+n/2))/%pi^((1-cos(%pi*n))/4);
 392
 393 /* Some small numbers */
 394
 395 double_factorial(-3);
 396 -1;
 397 errcatch(double_factorial(-2));
 398 [];
 399 double_factorial(-1);
 400 1;
 401 double_factorial(0);
 402 1;
 403 double_factorial(1);
 404 1;
 405 double_factorial(2);
 406 2;
 407 double_factorial(3);
 408 3;
 409 double_factorial(4);
 410 8;
 411 double_factorial(5);
 412 15;
 413 double_factorial(6);
 414 48;
 415 double_factorial(7);
 416 105;
 417 double_factorial(8);
 418 384;
 419 double_factorial(9);
 420 945;
 421 double_factorial(10);
 422 3840;
 423
 424 /* The same for double float */
 425
 426 closeto(
 427   double_factorial(-3.0),
 428   -1.0,
 429   1e-13);
 430 true;
 431
 432 errcatch(double_factorial(-2.0));
 433 [];
 434
 435 closeto(
 436   double_factorial(-1.0),
 437   1.0,
 438   1e-13);
 439 true;
 440
 441 closeto(
 442   double_factorial(0.0),
 443   1.0,
 444   1e-13);
 445 true;
 446
 447 closeto(
 448   double_factorial(1.0),
 449   1.0,
 450   1e-13);
 451 true;
 452
 453 closeto(
 454   double_factorial(2.0),
 455   2.0,
 456   1e-13);
 457 true;
 458
 459 closeto(
 460   double_factorial(3.0),
 461   3.0,
 462   1e-13);
 463 true;
 464
 465 closeto(
 466   double_factorial(4.0),
 467   8.0,
 468   1e-13);
 469 true;
 470
 471 closeto(
 472   double_factorial(5.0),
 473   15.0,
 474   1e-13);
 475 true;
 476
 477 closeto(
 478   double_factorial(6.0),
 479   48.0,
 480   1e-13);
 481 true;
 482
 483 closeto(
 484   double_factorial(7.0),
 485   105.0,
 486   1.279e-13);
 487 true;
 488
 489 closeto(
 490   double_factorial(8.0),
 491   384.0,
 492   1e-12);
 493 true;
 494
 495 closeto(
 496   double_factorial(9.0),
 497   945.0,
 498   1e-11);
 499 true;
 500
 501 closeto(
 502   double_factorial(10.0),
 503   3840.0,
 504   1.5E-12);
 505 true;
 506
 507 /* The same with bigfloat */
 508 fpprec:64;
 509 64;
 510
 511 closeto(
 512   double_factorial(-3.0b0),
 513   -1.0b0,
 514   1e-13);
 515 true;
 516
 517 errcatch(double_factorial(-2.0b0));
 518 [];
 519
 520 closeto(
 521   double_factorial(-1.0b0),
 522   1.0b0,
 523   1e-13);
 524 true;
 525
 526 closeto(
 527   double_factorial(0.0b0),
 528   1.0b0,
 529   1e-13);
 530 true;
 531
 532 closeto(
 533   double_factorial(1.0b0),
 534   1.0b0,
 535   1e-13);
 536 true;
 537
 538 closeto(
 539   double_factorial(2.0b0),
 540   2.0b0,
 541   1e-13);
 542 true;
 543
 544 closeto(
 545   double_factorial(3.0b0),
 546   3.0b0,
 547   1e-13);
 548 true;
 549
 550 closeto(
 551   double_factorial(4.0b0),
 552   8.0b0,
 553   1e-13);
 554 true;
 555
 556 closeto(
 557   double_factorial(5.0b0),
 558   15.0b0,
 559   1e-13);
 560 true;
 561
 562 closeto(
 563   double_factorial(6.0b0),
 564   48.0b0,
 565   1e-13);
 566 true;
 567
 568 closeto(
 569   double_factorial(7.0b0),
 570   105.0b0,
 571   1e-13);
 572 true;
 573
 574 closeto(
 575   double_factorial(8.0b0),
 576   384.0b0,
 577   1e-13);
 578 true;
 579
 580 closeto(
 581   double_factorial(9.0b0),
 582   945.0b0,
 583   1e-13);
 584 true;
 585
 586 closeto(
 587   double_factorial(10.0b0),
 588   3840.0b0,
 589   1e-13);
 590 true;
 591
 592 /* Some real and complex values */
 593
 594 closeto(
 595   double_factorial(-3.5),
 596   -1.283770376595223397225456287264697304361344685971440894669095353b0,
 597   1e-13);
 598 true;
 599
 600 closeto(
 601   double_factorial(-3.5b0),
 602   -1.283770376595223397225456287264697304361344685971440894669095353b0,
 603   1e-60);
 604 true;
 605
 606 closeto(
 607   double_factorial(-3.5+%i),
 608   (-0.0026442534512730229977827874410755514695008373007370518369259413b0
 609    +0.4140148090845355309500755922424659939330568167751526009311942842b0*%i),
 610   1e-13);
 611 true;
 612
 613 closeto(
 614   double_factorial(3.5),
 615   4.832319386136852665658314936437452651454869331098044546829825309b0,
 616   1e-13);
 617 true;
 618
 619 closeto(
 620   double_factorial(3.5b0),
 621   4.832319386136852665658314936437452651454869331098044546829825309b0,
 622   1e-60);
 623 true;
 624
 625 closeto(
 626   double_factorial(3.5+%i),
 627   (-2.165793510810110416038389252512222520262890874310470919228355939b0
 628    +4.032141259508464573377851775093179996368679285808989461893416849b0*%i),
 629   1e-13);
 630 true;
 631
 632 closeto(
 633   double_factorial(-3.5b0+%i),
 634   (-0.0026442534512730229977827874410755514695008373007370518369259413b0
 635    +0.4140148090845355309500755922424659939330568167751526009311942842b0*%i),
 636   1e-60);
 637 true;
 638
 639 closeto(
 640   double_factorial(3.5b0+%i),
 641   (-2.165793510810110416038389252512222520262890874310470919228355939b0
 642    +4.032141259508464573377851775093179996368679285808989461893416849b0*%i),
 643   1b-60);
 644 true;
 645
 646 closeto(
 647   double_factorial(3.3b0+%i),
 648   (-0.401169963963553982868990904015984192029807700247132080411340721b0
 649   + 1.778201955902329072861901606357849890890501421219437116360540910b0*%i),
 650   1b-60);
 651 true;
 652
 653 /******************************************************************************
 654
 655   Test the Gamma function
 656
 657   Numerical values are taken from functions/wolfram.com.
 658 ******************************************************************************/
 659
 660 /* The Gamma function has mirror symmetry */
 661
 662 declare(z,complex);
 663 done;
 664
 665 conjugate(gamma(z));
 666 gamma(conjugate(z));
 667
 668 conjugate(gamma(x+%i*y));
 669 gamma(x-%i*y);
 670
 671 /* Check some simple values for integer, float and bigfloat */
 672
 673 map('gamma,[1,2,3,4,5]);
 674 [1,1,2,6,24];
 675
 676 closeto(gamma(1.0),1.0,1e-13);
 677 true;
 678
 679 closeto(gamma(2.0),1.0,1e-13);
 680 true;
 681
 682 closeto(gamma(3.0),2.0,1e-13);
 683 true;
 684
 685 closeto(gamma(4.0),6.0,1e-13);
 686 true;
 687
 688 closeto(gamma(5.0),24.0,5e-13);
 689 true;
 690
 691 closeto(gamma(1.0b0),1.0b0,1e-13);
 692 true;
 693
 694 closeto(gamma(2.0b0),1.0b0,1e-13);
 695 true;
 696
 697 closeto(gamma(3.0b0),2.0b0,1e-13);
 698 true;
 699
 700 closeto(gamma(4.0b0),6.0b0,1e-13);
 701 true;
 702
 703 closeto(gamma(5.0b0),24.0b0,1e-13);
 704 true;
 705
 706 /* Check for a zero argument */
 707
 708 errcatch(gamma(0));
 709 [];
 710 errcatch(gamma(0.0));
 711 [];
 712 errcatch(gamma(0.0b0));
 713 [];
 714
 715 /* Check test for negative integer or a representation of a negative integer */
 716
 717 errcatch(gamma(-2));
 718 [];
 719 errcatch(gamma(-2.0));
 720 [];
 721 errcatch(gamma(-2.b0));
 722 [];
 723
 724 /* Check the correct handling of the $numer flag */
 725
 726 gamma(%e),numer;
 727 1.5674682557740529;
 728 gamma(%e),numer,%enumer;
 729 1.5674682557740529;
 730 gamma(a+b),numer;
 731 gamma(a+b);
 732 gamma(%i);
 733 gamma(%i);
 734 gamma(%i),numer;
 735 gamma(1.0*%i);   /* Evaluates to a complex number */
 736
 737 /* Check half integral integers as values */
 738
 739 gamma(1/2);
 740 sqrt(%pi);
 741 gamma(-1/2);
 742 -2*sqrt(%pi);
 743 gamma(3/2);
 744 sqrt(%pi)/2;
 745 gamma(-3/2);
 746 4*sqrt(%pi)/3;
 747 gamma(5/2);
 748 3*sqrt(%pi)/4;
 749 gamma(-5/2);
 750 -8*sqrt(%pi)/15;
 751
 752 /* Check expansion of the Gamma function */
 753
 754 gamma_expand:true;
 755 true;
 756
 757 gamma(z+1)/gamma(z);
 758 z;
 759
 760 gamma(gamma(z+1)/gamma(z));
 761 gamma(z);
 762
 763 gamma(z+1)/gamma(z-1);
 764 (z-1)*z;
 765
 766 gamma(z+2)/gamma(z-2);
 767 (z-2)*(z-1)*z*(z+1);
 768
 769 gamma_expand:false;
 770 false;
 771
 772 /* We check that the default values for $factlim:100000 and
 773    $gammalim:10000 work. */
 774
 775 fpprec:16;
 776 16;
 777
 778 relerror(
 779   bfloat(gamma(factlim)),
 780   gamma(bfloat(factlim)),
 781   1b-14);
 782 true;
 783
 784 relerror(
 785   bfloat(gamma((gammalim-1)+1/2)),
 786   gamma(bfloat((gammalim-1)+1/2)),
 787   1b-15);
 788 true;
 789
 790 relerror(
 791   bfloat(gamma(-gammalim+1/2)),
 792   gamma(bfloat(-gammalim+1/2)),
 793   1b-15);
 794 true;
 795
 796 /* Check test for overflow in flonum routine gamma-lanczos */
 797
 798 relerror(
 799   gamma(170.0),          /* should not overflow. For GCL 2.6.8 and */
 800   float(gamma(170)),     /* and CLISP 2.44 the limit is ~171.6243  */
 801   1e-12);
 802 true;
 803
 804 errcatch(gamma(175.0));  /* should overflow */
 805 [];
 806
 807 errcatch(gamma(250.0));  /* should overflow */
 808 [];
 809
 810 /* Simplifcation for infinities and undeterminates only for inf */
 811
 812 map('gamma, [inf,minf,infinity,und,ind]);
 813 [inf,gamma(minf),gamma(infinity),gamma(und),gamma(ind)];
 814
 815 /* Check real and complex arguments in double float precision.
 816    This is a check for the numerical routine gamma-lanczos */
 817
 818 closeto(
 819   gamma(1.5),
 820   0.8862269254527580136490837416705725913987747280611935641069038949b0,
 821   1e-14);
 822 true;
 823
 824 closeto(
 825   gamma(2.5),
 826   1.329340388179137020473625612505858887098162092091790346160355842b0,
 827   1e-14);
 828 true;
 829
 830 closeto(
 831   gamma(3.5),
 832   3.323350970447842551184064031264647217745405230229475865400889606b0,
 833   1e-14);
 834 true;
 835
 836 relerror(
 837   gamma(75.5),
 838   2.859942315653572214189951793671955438617013849084406338093590075b108,
 839   5e-14);
 840 true;
 841
 842 closeto(
 843   gamma(0.5+%i),
 844   (0.3006946172606558162173894638352104402306759641691949986162475934b0
 845   -0.4249678794331238126098496402574059704734842223340586518754297249b0*%i),
 846   1e-14);
 847 true;
 848
 849 closeto(
 850   gamma(1.5+%i),
 851   (0.5753151880634517207185443721750111905888222044186561511835535216b0
 852   +0.0882106775440939099124646437065074549939338530021656726785327309b0*%i),
 853   1e-14);
 854 true;
 855
 856 closeto(
 857   gamma(2.5+%i),
 858   (0.7747621045510836711653519145560093308892994536258185540967975514b0
 859   +0.7076312043795925855872413377347723730797229839219046602013526179b0*%i),
 860   1e-14);
 861 true;
 862
 863 closeto(
 864   gamma(3.5+%i),
 865   (1.229274056998116592326138448655250954143525650142641725040641261b0
 866   +2.543840115500065135133455258892940263588606913430580204600179096b0*%i),
 867   1e-14);
 868 true;
 869
 870 relerror(
 871   gamma(75.5+%i),
 872   (-1.092860022497734443706055997676557155572470037327121860702819811b108
 873    -2.622326961675321010452173874453854546607804545768376326095021243b108*%i),
 874   5.339e-14);
 875 true;
 876
 877 /* Check negative real arguments in double float precision.
 878    This is a check for the reflection formula of gamma-lanzos */
 879
 880 closeto(
 881   gamma(-0.5),
 882   -3.544907701811032054596334966682290365595098912244774256427615580b0,
 883   1e-14);
 884 true;
 885
 886 closeto(
 887   gamma(-1.5),
 888   2.363271801207354703064223311121526910396732608163182837618410386b0,
 889   1e-14);
 890 true;
 891
 892 closeto(
 893   gamma(-2.5),
 894   -0.9453087204829418812256893244486107641586930432652731350473641546b0,
 895   1e-14);
 896 true;
 897
 898 /* Check real arguments up to 64 digits.
 899    This is a check for the numerical routine bffac */
 900
 901 fpprec:64; /* we have saved the actual value at the beginning of the file */
 902 64;
 903
 904 closeto(
 905   gamma(0.5b0),
 906   1.772453850905516027298167483341145182797549456122387128213807790b0,
 907   5.0b-64);
 908 true;
 909
 910 closeto(
 911   gamma(1.5b0),
 912   0.8862269254527580136490837416705725913987747280611935641069038949b0,
 913   5.0b-64);
 914 true;
 915
 916 closeto(
 917   gamma(2.5b0),
 918   1.329340388179137020473625612505858887098162092091790346160355842b0,
 919   5.0b-64);
 920 true;
 921
 922 closeto(
 923   gamma(3.5b0),
 924   3.323350970447842551184064031264647217745405230229475865400889606b0,
 925   5.0b-64);
 926 true;
 927
 928 relerror(
 929   gamma(75.5b0),
 930   2.859942315653572214189951793671955438617013849084406338093590075b108,
 931   5.0b-64);
 932 true;
 933
 934 /* Check negative real arguments up to 64 digits.
 935    This is a check for the reflection formula of bffac */
 936
 937 closeto(
 938   gamma(-0.5b0),
 939   -3.544907701811032054596334966682290365595098912244774256427615580b0,
 940   5.0b-64);
 941 true;
 942
 943 closeto(
 944   gamma(-1.5b0),
 945   2.363271801207354703064223311121526910396732608163182837618410386b0,
 946   8.0b-64);
 947 true;
 948
 949 closeto(
 950   gamma(-2.5b0),
 951   -0.9453087204829418812256893244486107641586930432652731350473641546b0,
 952   5.0b-64);
 953 true;
 954
 955 /* Check complex arguments up to 64 digits.
 956    This is a check for the numerical routine cbffac */
 957
 958 closeto(
 959   gamma(0.5b0+%i),
 960   (0.3006946172606558162173894638352104402306759641691949986162475934b0
 961   -0.4249678794331238126098496402574059704734842223340586518754297249b0*%i),
 962   5.0b-64);
 963 true;
 964
 965 closeto(
 966   gamma(1.5b0+%i),
 967   (0.5753151880634517207185443721750111905888222044186561511835535216b0
 968   +0.0882106775440939099124646437065074549939338530021656726785327309b0*%i),
 969   5.0b-64);
 970 true;
 971
 972 closeto(
 973   gamma(2.5b0+%i),
 974   (0.7747621045510836711653519145560093308892994536258185540967975514b0
 975   +0.7076312043795925855872413377347723730797229839219046602013526179b0*%i),
 976   5.0b-64);
 977 true;
 978
 979 closeto(
 980   gamma(3.5b0+%i),
 981   (1.229274056998116592326138448655250954143525650142641725040641261b0
 982   +2.543840115500065135133455258892940263588606913430580204600179096b0*%i),
 983   5.0b-64);
 984 true;
 985
 986 relerror(
 987   gamma(75.5b0+%i),
 988   (-1.092860022497734443706055997676557155572470037327121860702819811b108
 989    -2.622326961675321010452173874453854546607804545768376326095021243b108*%i),
 990   5.0b-64);
 991 true;
 992
 993 (fpprec:oldfpprec,done); /* Reset the value of fpprec */
 994 done;
 995
 996 /******************************************************************************
 997
 998   Test the Incomplete Gamma function
 999
1000 ******************************************************************************/
1001
1002 (oldfpprec : fpprec,done);
1003 done;
1004
1005 /* Some special values */
1006
1007 gamma_incomplete(a,0);
1008 gamma_incomplete(a,0);
1009
1010 (assume(am < 0, ap > 0),done);
1011 done;
1012
1013 errcatch(gamma_incomplete(-1,0));
1014 [];
1015 errcatch(gamma_incomplete(-2,0));
1016 [];
1017 errcatch(gamma_incomplete(am,0));
1018 [];
1019 gamma_incomplete(1,0);
1020 1;
1021 gamma_incomplete(2,0);
1022 1;
1023 gamma_incomplete(3,0);
1024 2;
1025 gamma_incomplete(ap,0);
1026 gamma(ap);
1027
1028 gamma_incomplete(a,inf);
1029 0;
1030
1031 /* Expansion of the Incomplete Gamma function */
1032
1033 gamma_expand:true;
1034 true;
1035
1036 gamma_incomplete(0,z);
1037 -expintegral_ei(-z)+1/2*(log(-z)-log(-1/z))-log(z);
1038
1039 gamma_incomplete(1/2,z);
1040 sqrt(%pi)*erfc(sqrt(z));
1041
1042 gamma_incomplete(-1/2,z);
1043 2*%e^(-z)/sqrt(z)-2*sqrt(%pi)*erfc(sqrt(z));
1044
1045 gamma_incomplete(1,z);
1046 %e^-z;
1047
1048 gamma_incomplete(-1,z);
1049 log(z)-(log(-z)-log(-1/z))/2+expintegral_ei(-z)+%e^-z/z;
1050
1051 gamma_incomplete(3/2,z);
1052 sqrt(z)*%e^-z+sqrt(%pi)*erfc(sqrt(z))/2;
1053
1054 gamma_incomplete(-3/2,z);
1055 4*sqrt(%pi)*erfc(sqrt(z))/3-(4*z/3-2/3)*%e^-z/z^(3/2);
1056
1057 gamma_incomplete(2,z);
1058 (z+1)*%e^-z;
1059
1060 gamma_incomplete(a+1,z);
1061 z^a*%e^-z+a*gamma_incomplete(a,z);
1062
1063 gamma_incomplete(a-1,z);
1064 -z^(a-1)*%e^-z/(a-1)-gamma_incomplete(a,z)/(1-a);
1065
1066 gamma_incomplete(a+2,z);
1067 z^a*(z+a+1)*%e^-z+a*(a+1)*gamma_incomplete(a,z);
1068
1069 gamma_incomplete(a-2,z);
1070 gamma_incomplete(a,z)/((1-a)*(2-a))-z^(a-2)*(z/((a-2)*(a-1))+1/(a-2))*%e^-z;
1071
1072 gamma_expand:false;
1073 false;
1074
1075 /* The Incomplete Gamma function has not mirror symmetry on the negative
1076    real axis. We have supported a conjugate-gamma-incomplete function */
1077
1078 declare(ac, complex, zc,complex);
1079 done;
1080
1081 conjugate(gamma_incomplete(ac,zc));
1082 conjugate(gamma_incomplete(ac,zc));
1083
1084 conjugate(gamma_incomplete(a+b*%i,x+y*%i));
1085 conjugate(gamma_incomplete(a+b*%i,x+y*%i));
1086
1087 conjugate(gamma_incomplete(a+b*%i,10));
1088 gamma_incomplete(a-b*%i,10);
1089
1090 conjugate(gamma_incomplete(a+b*%i,-10));
1091 conjugate(gamma_incomplete(a+b*%i,-10));
1092
1093 /* Derivatives of the Incomplete Gamma function */
1094
1095 diff(gamma_incomplete(a,x),x);
1096 -x^(a-1)*%e^(-x);
1097
1098 /* Noun form if sign of the parameter a is not known */
1099 diff(gamma_incomplete(a,x),a);
1100 'diff(gamma_incomplete(a,x),a);
1101
1102 /* The parameter a has a positive sign */
1103 (assume(a>0),done);
1104 done;
1105
1106 expand(diff(gamma_incomplete(a,x),a)-
1107 (-(gamma(a)-gamma_incomplete(a,x))*log(x)
1108        +gamma(a)^2*hypergeometric_regularized([a,a],[a+1,a+1],-x)*x^a
1109        +psi[0](a)*gamma(a)));
1110 0;
1111
1112 (forget(a>0),done);
1113 done;
1114
1115 /* Numerical tests for the Incomplete Gamma function */
1116
1117 /* At first tests for a=0 and negative integers for a
1118    For this values the numerical algorithm of gamma-incomplete do not
1119    work. The routines of the Exponential Integral E are used. */
1120
1121 closeto(
1122   gamma_incomplete(0,0.5),
1123   0.55977359477616081174679593931509b0,
1124   5e-15);
1125 true;
1126
1127 closeto(
1128   gamma_incomplete(0,-0.5),
1129   -0.4542199048631735799205238126628b0-3.1415926535897932384626433832795b0*%i,
1130   5e-15);
1131 true;
1132
1133 closeto(
1134   gamma_incomplete(0,0.5+%i),
1135   -0.07139471104245272355588497993684b0-0.35749377365216265125485869345732b0*%i,
1136   5e-15);
1137 true;
1138
1139 closeto(
1140   gamma_incomplete(0,-0.5+%i),
1141   -0.92289919055678882179364241497461b0-0.81571273343182452677967141350955b0*%i,
1142   5e-15);
1143 true;
1144
1145 /* Now for bigfloat precision */
1146
1147 fpprec:32;
1148 32;
1149
1150 closeto(
1151   gamma_incomplete(0,0.5b0),
1152   0.55977359477616081174679593931509b0,
1153   1b-30);
1154 true;
1155
1156 closeto(
1157   gamma_incomplete(0,-0.5b0),
1158   -0.4542199048631735799205238126628b0-3.1415926535897932384626433832795b0*%i,
1159   1b-30);
1160 true;
1161
1162 closeto(
1163   gamma_incomplete(0,0.5b0+%i),
1164   -0.07139471104245272355588497993684b0-0.35749377365216265125485869345732b0*%i,
1165   1b-30);
1166 true;
1167
1168 closeto(
1169   gamma_incomplete(0,-0.5b0+%i),
1170   -0.92289919055678882179364241497461b0-0.81571273343182452677967141350955b0*%i,
1171   1b-30);
1172 true;
1173
1174 /* Tests with negative integers */
1175
1176 closeto(
1177   gamma_incomplete(-1,-0.5),
1178   -2.8432226365370827137767777629655b0+3.1415926535897932384626433832795b0*%i,
1179   5e-15);
1180 true;
1181
1182 closeto(
1183   gamma_incomplete(-5,-0.5),
1184   -12.1692735494620863710665339515891b0+0.0261799387799149436538553615273b0*%i,
1185   8e-15);
1186 true;
1187
1188 closeto(
1189   gamma_incomplete(-1,-0.5+%i),
1190   -0.54330486022427331058655939975980b0+0.65800685452922697341354545991819b0*%i,
1191   5e-15);
1192 true;
1193
1194 closeto(
1195   gamma_incomplete(-5,-0.5+%i),
1196   0.08881923337904654547301590363812b0+0.18155943299796684573147338017934b0*%i,
1197   5e-15);
1198 true;
1199
1200 /* Again for bigfloat precision */
1201
1202 closeto(
1203   gamma_incomplete(-1,-0.5b0),
1204   -2.8432226365370827137767777629655b0+3.1415926535897932384626433832795b0*%i,
1205   1b-30);
1206 true;
1207
1208 closeto(
1209   gamma_incomplete(-5,-0.5b0),
1210   -12.1692735494620863710665339515891b0+0.0261799387799149436538553615273b0*%i,
1211   1b-30);
1212 true;
1213
1214 closeto(
1215   gamma_incomplete(-1,-0.5b0+%i),
1216   -0.54330486022427331058655939975980b0+0.65800685452922697341354545991819b0*%i,
1217   1b-30);
1218 true;
1219
1220 closeto(
1221   gamma_incomplete(-5,-0.5b0+%i),
1222   0.08881923337904654547301590363812b0+0.18155943299796684573147338017934b0*%i,
1223   1b-30);
1224 true;
1225
1226 /* Test gamma_incomplete(0.25,2.5) for Float and Bigfloat */
1227
1228 closeto(
1229   gamma_incomplete(0.25,2.5),
1230   0.03340545777928488523612480546612030546638337899458717728445920914b0,
1231   1e-15);
1232 true;
1233
1234 fpprec:32;
1235 32;
1236
1237 closeto(
1238   gamma_incomplete(0.25b0,2.5b0),
1239   0.03340545777928488523612480546612030546638337899458717728445920914b0,
1240   1b-31);
1241 true;
1242
1243 fpprec:64;
1244 64;
1245
1246 closeto(
1247   gamma_incomplete(0.25b0,2.5b0),
1248   0.03340545777928488523612480546612030546638337899458717728445920914b0,
1249   1b-63);
1250 true;
1251
1252 /* Test gamma_incomplete(0.25,0.25) for Float and Bigfloat */
1253
1254 closeto(
1255   gamma_incomplete(0.25,0.25),
1256   0.9293237832774184425973508042578251762794944752213875213176435274b0,
1257   5.0e-15);
1258 true;
1259
1260 fpprec:32;
1261 32;
1262
1263 closeto(
1264   gamma_incomplete(0.25b0,0.25b0),
1265   0.9293237832774184425973508042578251762794944752213875213176435274b0,
1266   1.1b-31);
1267 true;
1268
1269 fpprec:64;
1270 64;
1271
1272 closeto(
1273   gamma_incomplete(0.25b0,0.25b0),
1274   0.9293237832774184425973508042578251762794944752213875213176435274b0,
1275   1.1b-63);
1276 true;
1277
1278 /* Test gamma_incomplete(0.25,0.50) for Float and Bigfloat */
1279
1280 closeto(
1281   gamma_incomplete(0.25,0.50),
1282   0.55658041400942713438787175086207b0,
1283   5.0e-15);
1284 true;
1285
1286 fpprec:32;
1287 32;
1288
1289 closeto(
1290   gamma_incomplete(0.25b0,0.50b0),
1291   0.55658041400942713438787175086207b0,
1292   1b-31);
1293 true;
1294
1295 fpprec:64;
1296 64;
1297
1298 closeto(
1299   gamma_incomplete(0.25b0,0.50b0),
1300   0.5565804140094271343878717508620650091658338999776480841533264361b0,
1301   2b-63);
1302 true;
1303
1304 fpprec:128;
1305 128;
1306
1307 closeto(
1308   gamma_incomplete(0.25b0,0.50b0),
1309   0.55658041400942713438787175086206500916583389997764808415332643613122015052649897833312327325822333229784708198027750127190766504b0,
1310   1.0b-127);
1311 true;
1312
1313 /* Test gamma_incomplete(0.25,1.50) for Float and Bigfloat */
1314
1315 closeto(
1316   gamma_incomplete(0.25,1.50),
1317   0.12115499104033848614860340878369b0,
1318   1e-15);
1319 true;
1320
1321 fpprec:32;
1322 32;
1323
1324 closeto(
1325   gamma_incomplete(0.25b0,1.50b0),
1326   0.12115499104033848614860340878369b0,
1327   1b-31);
1328 true;
1329
1330 fpprec:64;
1331 64;
1332
1333 closeto(
1334   gamma_incomplete(0.25b0,1.50b0),
1335   0.1211549910403384861486034087836891246955052387140720625064500006b0,
1336   1.5b-63);
1337 true;
1338
1339 fpprec:128;
1340 128;
1341
1342 closeto(
1343   gamma_incomplete(0.25b0,1.50b0),
1344   0.12115499104033848614860340878368912469550523871407206250645000059332022509505923467877887847273887882437030555876962014143410940b0,
1345   4b-127);
1346 true;
1347
1348 closeto(
1349   gamma_incomplete(0.25b0,1.50b0),
1350   1.5b0^0.25b0*expintegral_e(1.0b0-0.25b0,1.50b0),
1351   5b-128);
1352 true;
1353
1354 fpprec:34; /* Two extra digits to get 32 digits in the following tests */
1355 34;
1356
1357 relerror(
1358   gamma_incomplete(1000b0,1000b0),
1359   1.995014933549148239529838438260433407652488769526598301696165147b2564,
1360   1b-32);
1361 true;
1362
1363 relerror(
1364   gamma_incomplete(1000b0,100b0),
1365   4.023872600770937735437024339230039857193748642107146325437999104b2564,
1366   1b-32);
1367 true;
1368
1369 relerror(
1370   gamma_incomplete(1000b0,10b0),
1371   4.023872600770937735437024339230039857193748642107146325437999104b2564,
1372   1b-32);
1373 true;
1374
1375 relerror(
1376   gamma_incomplete(1000b0,1b0),
1377   4.023872600770937735437024339230039857193748642107146325437999104b2564,
1378   1b-32);
1379 true;
1380
1381 relerror(
1382   gamma_incomplete(1000b0,-1b0),
1383   4.023872600770937735437024339230039857193748642107146325437999104b2564,
1384   1b-32);
1385 true;
1386
1387 relerror(
1388   gamma_incomplete(1000b0,-10b0),
1389   4.023872600770937735437024339230039857193748642107146325437999104b2564,
1390   1b-32);
1391 true;
1392
1393 relerror(
1394   gamma_incomplete(1000b0,-100b0),
1395   4.023872600770937735437024339230039857193748642107146325437999104b2564,
1396   1b-32);
1397 true;
1398
1399 relerror(
1400   gamma_incomplete(1000b0,-1000b0),
1401   -9.852818774470566937423668137175694874333788729537950495924821627b3430,
1402   1b-32);
1403 true;
1404
1405 relerror(
1406   gamma_incomplete(100b0,100b0),
1407   4.542198120862669429369147083086235039624517049342017449058357596b155,
1408   1b-32);
1409 true;
1410
1411 relerror(
1412   gamma_incomplete(100b0,10b0),
1413   9.332621544394415268169923885626670049071596826438162146859296339b155,
1414   1b-32);
1415 true;
1416
1417 relerror(
1418   gamma_incomplete(100b0,1b0),
1419   9.332621544394415268169923885626670049071596826438162146859296390b155,
1420   1b-32);
1421 true;
1422
1423 relerror(
1424   gamma_incomplete(100b0,-1b0),
1425   9.332621544394415268169923885626670049071596826438162146859296390b155,
1426   1b-32);
1427 true;
1428
1429 relerror(
1430   gamma_incomplete(100b0,-10b0),
1431   9.332621544394415268169923885626670049071596826438162126818796880b155,
1432   1b-32);
1433 true;
1434
1435 relerror(
1436   gamma_incomplete(100b0,-100b0),
1437   -1.3474270960118181325667224386845432493096383414519386259680854024b241,
1438   1b-32);
1439 true;
1440
1441 relerror(
1442   gamma_incomplete(10.0+10.0*%i,10.0+10.0*%i),
1443   (712.747910954771249931938579893612285083502899995529160358791610b0
1444   -1614.519712336984904341104157868496978481416095290952330318983747b0*%i),
1445   7.5e-15);
1446 true;
1447
1448 relerror(
1449   gamma_incomplete(10b0+10b0*%i,10b0+10b0*%i),
1450   (712.747910954771249931938579893612285083502899995529160358791610b0
1451   -1614.519712336984904341104157868496978481416095290952330318983747b0*%i),
1452   4b-32);
1453 true;
1454
1455 relerror(
1456   gamma_incomplete(10.0+10*%i,10.0+5*%i),
1457   (3795.479456353067145208395441052660229834399956460948716792241863b0
1458   -1859.399046776284485239753978633491801182777480033526406270435152b0*%i),
1459   6.0e-15);
1460 true;
1461
1462 relerror(
1463   gamma_incomplete(10.0b0+10*%i,10.0b0+5*%i),
1464   (3795.479456353067145208395441052660229834399956460948716792241863b0
1465   -1859.399046776284485239753978633491801182777480033526406270435152b0*%i),
1466   1b-32);
1467 true;
1468
1469 relerror(
1470   gamma_incomplete(10.0+5*%i,10.0+5*%i),
1471   (22616.57428441264599471916533645601396385068769401974320192387776b0
1472   -41760.26634389514228374497096679850877647381173070930602051580693b0*%i),
1473   5.0e-15);
1474 true;
1475
1476 relerror(
1477   gamma_incomplete(10.0b0+5*%i,10.0b0+5*%i),
1478   (22616.57428441264599471916533645601396385068769401974320192387776b0
1479   -41760.26634389514228374497096679850877647381173070930602051580693b0*%i),
1480   1b-32);
1481 true;
1482
1483 relerror(
1484   gamma_incomplete(10+5*%i,10+2.5*%i),
1485   (55884.99767768350551452192526458363894624371018195106017631282130b0
1486   -30587.35558698211815103119732529095917842073159139555085583572089b0*%i),
1487   8.0e-15);
1488 true;
1489
1490 relerror(
1491   gamma_incomplete(10b0+5*%i,10b0+2.5*%i),
1492   (55884.99767768350551452192526458363894624371018195106017631282130b0
1493   -30587.35558698211815103119732529095917842073159139555085583572089b0*%i),
1494   1b-32);
1495 true;
1496
1497 relerror(
1498   gamma_incomplete(10.0+2.5*%i,10.0+2.5*%i),
1499   (98307.31859173691954817642978681043594336734907098079356276769738b0
1500   -69378.82767710646665454093742183442049572498499915146277510648781b0*%i),
1501   7.0e-15);
1502 true;
1503
1504 relerror(
1505   gamma_incomplete(10.0b0+2.5*%i,10.0b0+2.5*%i),
1506   (98307.31859173691954817642978681043594336734907098079356276769738b0
1507   -69378.82767710646665454093742183442049572498499915146277510648781b0*%i),
1508   1b-32);
1509 true;
1510
1511 relerror(
1512   gamma_incomplete(10.0+2.5*%i,10.0+1.5*%i),
1513   (119713.97958915216843109406780063078781556428789769599762881675530b0
1514   -44021.05551717694140528840726282083152859527358436513276174120234b0*%i),
1515   9.5e-15);
1516 true;
1517
1518 relerror(
1519   gamma_incomplete(10.0b0+2.5*%i,10.0b0+1.5*%i),
1520   (119713.97958915216843109406780063078781556428789769599762881675530b0
1521   -44021.05551717694140528840726282083152859527358436513276174120234b0*%i),
1522   1b-32);
1523 true;
1524
1525 relerror(
1526   gamma_incomplete(10.0+1.5*%i,10.0+1.5*%i),
1527   (-143260.5455945276009736506823530548923946268185687353440779787855b0
1528   -36427.8601104063533811294405176711748076293661563594190640674077b0*%i),
1529   1.02e-14);
1530 true;
1531
1532 relerror(
1533   gamma_incomplete(10.0b0+1.5*%i,10.0b0+1.5*%i),
1534   (-143260.5455945276009736506823530548923946268185687353440779787855b0
1535   -36427.8601104063533811294405176711748076293661563594190640674077b0*%i),
1536   1b-32);
1537 true;
1538
1539 relerror(
1540   gamma_incomplete(10.0+1.5*%i,10.0+0.5*%i),
1541   (-134422.2837349310843015830622649231296922981730868773827217434186b0
1542   -76495.4696532860249045908863041952283847028641854696387454419024b0*%i),
1543   1.271e-14);
1544 true;
1545
1546 relerror(
1547   gamma_incomplete(10.0b0+1.5*%i,10.0b0+0.5*%i),
1548   (-134422.2837349310843015830622649231296922981730868773827217434186b0
1549   -76495.4696532860249045908863041952283847028641854696387454419024b0*%i),
1550   1b-32);
1551 true;
1552
1553 relerror(
1554   gamma_incomplete(10+0.5*%i,10+0.5*%i),
1555   (70217.4190738045440197722508789471776325390496700488861347101465b0
1556   +148228.5649085354026330085827685305718428914334429314381378279173b0*%i),
1557   6.0e-15);
1558 true;
1559
1560 relerror(
1561   gamma_incomplete(10b0+0.5*%i,10b0+0.5*%i),
1562   (70217.4190738045440197722508789471776325390496700488861347101465b0
1563   +148228.5649085354026330085827685305718428914334429314381378279173b0*%i),
1564   1b-32);
1565 true;
1566
1567 relerror(
1568   gamma_incomplete(5.0+5*%i,5.0+5*%i),
1569   (-0.4806117328699535298510981197039733622773799503543787399412087606b0
1570    +0.8919199556012029365414433316086474496955099800079200844577588174b0*%i),
1571   5.0e-15);
1572 true;
1573
1574 relerror(
1575   gamma_incomplete(5.0b0+5*%i,5.0b0+5*%i),
1576   (-0.4806117328699535298510981197039733622773799503543787399412087606b0
1577    +0.8919199556012029365414433316086474496955099800079200844577588174b0*%i),
1578   1b-32);
1579 true;
1580
1581 relerror(
1582   gamma_incomplete(5.0+5*%i,5.0+2.5*%i),
1583   (-1.564618625118515702134408016446776958254116141951281337045968096b0
1584    +0.763213115623024892186357289469889653460737525026968402880536542b0*%i),
1585   5.0e-15);
1586 true;
1587
1588 relerror(
1589   gamma_incomplete(5.0b0+5*%i,5.0b0+2.5*%i),
1590   (-1.564618625118515702134408016446776958254116141951281337045968096b0
1591    +0.763213115623024892186357289469889653460737525026968402880536542b0*%i),
1592   1b-32);
1593 true;
1594
1595 relerror(
1596   gamma_incomplete(5.0+2.5*%i,5.0+2.5*%i),
1597   (-3.966094476128530812031476427059327525923710898646502835891758741b0
1598    -3.843825405733108837026233752472379492035828310495276258481876760b0*%i),
1599   5.0e-15);
1600 true;
1601
1602 relerror(
1603   gamma_incomplete(5.0b0+2.5*%i,5.0b0+2.5*%i),
1604   (-3.966094476128530812031476427059327525923710898646502835891758741b0
1605    -3.843825405733108837026233752472379492035828310495276258481876760b0*%i),
1606   1b-32);
1607 true;
1608
1609 relerror(
1610   gamma_incomplete(5.0+2.5*%i,5.0-2.5*%i),
1611   (24.28851242625584709660092800462769890376139980405107830536077980b0
1612   -13.30717877353455881129062577273028018783505328101266063709140224b0*%i),
1613   2.0e-15);
1614 true;
1615
1616 relerror(
1617   gamma_incomplete(5.0b0+2.5*%i,5.0b0-2.5*%i),
1618   (24.28851242625584709660092800462769890376139980405107830536077980b0
1619   -13.30717877353455881129062577273028018783505328101266063709140224b0*%i),
1620   1b-32);
1621 true;
1622
1623 relerror(
1624   gamma_incomplete(5.0-2.5*%i,5.0-2.5*%i),
1625   (-3.966094476128530812031476427059327525923710898646502835891758741b0
1626    +3.843825405733108837026233752472379492035828310495276258481876760b0*%i),
1627   5.0e-15);
1628 true;
1629
1630 relerror(
1631   gamma_incomplete(5.0b0-2.5*%i,5.0b0-2.5*%i),
1632   (-3.966094476128530812031476427059327525923710898646502835891758741b0
1633    +3.843825405733108837026233752472379492035828310495276258481876760b0*%i),
1634   1b-32);
1635 true;
1636
1637 /* Further tests after modication of the code */
1638
1639 fpprec:35;
1640 35;
1641
1642 /* We start with values on the real negative axis */
1643
1644 /* Problem 6 */
1645 relerror(
1646   gamma_incomplete(0.5,-10),
1647   1.7724538509055160272981674833 - 7388.5381938108184552671664573665*%i,
1648   2.0e-12);
1649 true;
1650
1651 /* Problem 7 */
1652 relerror(
1653   gamma_incomplete(0.5,-10b0),
1654   1.7724538509055160272981674833b0 - 7388.5381938108184552671664573665b0*%i,
1655   3.0b-32);
1656 true;
1657
1658 relerror(
1659   gamma_incomplete(0.5,-5),
1660   1.772453850905516027298167483341 - 76.796224205322062453935496965541*%i,
1661   1.5e-14);
1662 true;
1663
1664 relerror(
1665   gamma_incomplete(0.5,-5b0),
1666   1.772453850905516027298167483341b0 - 76.796224205322062453935496965541b0*%i,
1667   1b-32);
1668 true;
1669
1670 relerror(
1671   gamma_incomplete(0.5,-2.5),
1672   1.7724538509055160272981674833411 - 9.8735082388772780413725529343873*%i,
1673   1e-15);
1674 true;
1675
1676 relerror(
1677   gamma_incomplete(0.5,-2.5b0),
1678   1.7724538509055160272981674833411b0 - 9.8735082388772780413725529343873b0*%i,
1679   1b-32);
1680 true;
1681
1682 relerror(
1683   gamma_incomplete(0.5,-0.25),
1684   1.7724538509055160272981674833411 - 1.0899742083672444473248402628140*%i,
1685   5.0e-15);
1686 true;
1687
1688 relerror(
1689   gamma_incomplete(0.5,-0.25b0),
1690   1.7724538509055160272981674833411b0 - 1.0899742083672444473248402628140b0*%i,
1691   5.0b-32);
1692 true;
1693
1694 /* We add a small imaginary part
1695    The continued fraction will fail, Maxima use the series expansion.
1696 */
1697
1698 relerror(
1699   gamma_incomplete(0.5,-10+0.1*%i),
1700   -693.7125259652496257225339009893 - 7355.4778982854340599489999722765*%i,
1701   2.0e-12);
1702 true;
1703
1704 relerror(
1705   gamma_incomplete(0.5,-10.0b0+0.1b0*%i),
1706   -693.7125259652496257225339009893b0 - 7355.4778982854340599489999722765b0*%i,
1707   1b-32);
1708 true;
1709
1710 relerror(
1711   gamma_incomplete(0.5,-5+0.1*%i),
1712   -4.855606925906958733850786731743 - 76.497762747671653356926921414108*%i,
1713   5.0e-15);
1714 true;
1715
1716 relerror(
1717   gamma_incomplete(0.5,-5b0+0.1b0*%i),
1718   -4.855606925906958733850786731743b0 - 76.497762747671653356926921414108b0*%i,
1719   1b-32);
1720 true;
1721
1722 relerror(
1723   gamma_incomplete(0.5,-2.5+0.1*%i),
1724   1.0028894769544495897063089678047 - 9.8427092366892270639885076130022*%i,
1725   1.21e-15);
1726 true;
1727
1728 relerror(
1729   gamma_incomplete(0.5,-2.5b0+0.1b0*%i),
1730   1.0028894769544495897063089678047b0 - 9.8427092366892270639885076130022b0*%i,
1731   1b-32);
1732 true;
1733
1734 relerror(
1735   gamma_incomplete(0.5,-0.25+0.1*%i),
1736   1.5192534335558569724359463295125 - 1.1019363602857804797837310045806*%i,
1737   2.0e-15);
1738 true;
1739
1740 relerror(
1741   gamma_incomplete(0.5,-0.25b0+0.1b0*%i),
1742   1.5192534335558569724359463295125b0 - 1.1019363602857804797837310045806b0*%i,
1743   2.0b-32);
1744 true;
1745
1746 /* We add an imaginary part above the threshold for the series expansion
1747 */
1748
1749 relerror(
1750   gamma_incomplete(0.5,-10+1.1*%i),
1751   -6334.1948696342193438649689587744 - 3741.5383594183284610691163553403*%i,
1752   4.0e-12);
1753 true;
1754
1755 relerror(
1756   gamma_incomplete(0.5,-10b0+1.1b0*%i),
1757   -6334.1948696342193438649689587744b0 - 3741.5383594183284610691163553403b0*%i,
1758   5.0b-32);
1759 true;
1760
1761 relerror(
1762   gamma_incomplete(0.5,-5b0+1.1b0*%i),
1763   -59.650578505004222281783394220768b0 - 43.683456195166750403224017552560b0*%i,
1764   5.0b-32);
1765 true;
1766
1767 relerror(
1768   gamma_incomplete(0.5,-2.5+1.1*%i),
1769   -5.5335063099201920292609023851020 - 6.4351290837917628991361018359964*%i,
1770   5.0e-15);
1771 true;
1772
1773 relerror(
1774   gamma_incomplete(0.5,-2.5b0+1.1b0*%i),
1775   -5.5335063099201920292609023851020b0 - 6.4351290837917628991361018359964b0*%i,
1776   1b-32);
1777 true;
1778
1779 /* This is the serious expansion which works */
1780
1781 relerror(
1782   gamma_incomplete(0.5,-0.25+1.1*%i),
1783   -0.13794253952035726152165921998252 - 1.06583290164642973789895384973429*%i,
1784   2.0e-15);
1785 true;
1786
1787 relerror(
1788   gamma_incomplete(0.5,-0.25b0+1.1b0*%i),
1789   -0.13794253952035726152165921998252b0 - 1.06583290164642973789895384973429b0*%i,
1790   1b-32);
1791 true;
1792
1793 /* Values on the imaginary axis */
1794
1795 relerror(
1796   gamma_incomplete(0.5,-100*%i),
1797   0.096899159215326861150517776631041 + 0.024684086404223368298902429118542*%i,
1798   1e-15);
1799 true;
1800
1801 relerror(
1802   gamma_incomplete(0.5,-100b0*%i),
1803   0.096899159215326861150517776631041b0 + 0.024684086404223368298902429118542b0*%i,
1804   1b-32);
1805 true;
1806
1807 relerror(
1808   gamma_incomplete(0.5,-50*%i),
1809   0.123398939447119035626603437814684 + 0.069012601491650179241716316049576*%i,
1810   1e-15);
1811 true;
1812
1813 relerror(
1814   gamma_incomplete(0.5,-50b0*%i),
1815   0.123398939447119035626603437814684b0 + 0.069012601491650179241716316049576b0*%i,
1816   1b-32);
1817 true;
1818
1819 relerror(
1820   gamma_incomplete(0.5,-10*%i),
1821   -0.08046978022707502678937802264904 - 0.30392674968841293510316822780670*%i,
1822   1e-15);
1823 true;
1824
1825 relerror(
1826   gamma_incomplete(0.5,-10b0*%i),
1827   -0.08046978022707502678937802264904b0 - 0.30392674968841293510316822780670b0*%i,
1828   2.0b-32);
1829 true;
1830
1831 relerror(
1832   gamma_incomplete(0.5,-5*%i),
1833   0.36441984106355895337750863822965 - 0.24368559063811288395048612967632*%i,
1834   3.0e-15);
1835 true;
1836
1837 relerror(
1838   gamma_incomplete(0.5,-5b0*%i),
1839   0.36441984106355895337750863822965b0 - 0.24368559063811288395048612967632b0*%i,
1840   1b-32);
1841 true;
1842
1843 relerror(
1844   gamma_incomplete(0.5,-2.5*%i),
1845   -0.59691417904238855062194720247331 + 0.00921495731742953647951029973386*%i,
1846   3.0788e-15);
1847 true;
1848
1849 relerror(
1850   gamma_incomplete(0.5,-2.5b0*%i),
1851   -0.59691417904238855062194720247331b0 + 0.00921495731742953647951029973386b0*%i,
1852   1b-32);
1853 true;
1854
1855 relerror(
1856   gamma_incomplete(0.5,-0.25*%i),
1857   1.01109069076165681623650036780470 + 0.64403710594044452447886365988086*%i,
1858   1e-15);
1859 true;
1860
1861 relerror(
1862   gamma_incomplete(0.5,-0.25b0*%i),
1863   1.01109069076165681623650036780470b0 + 0.64403710594044452447886365988086b0*%i,
1864   1b-32);
1865 true;
1866
1867 relerror(
1868   gamma_incomplete(0.5,0.25*%i),
1869   1.01109069076165681623650036780470 - 0.64403710594044452447886365988086*%i,
1870   1e-15);
1871 true;
1872
1873 relerror(
1874   gamma_incomplete(0.5,0.25b0*%i),
1875   1.01109069076165681623650036780470b0 - 0.64403710594044452447886365988086b0*%i,
1876   1b-32);
1877 true;
1878
1879 relerror(
1880   gamma_incomplete(0.5,2.5*%i),
1881   -0.59691417904238855062194720247331 - 0.00921495731742953647951029973386*%i,
1882  3.0788e-15);
1883 true;
1884
1885 relerror(
1886   gamma_incomplete(0.5,2.5b0*%i),
1887   -0.59691417904238855062194720247331b0 - 0.00921495731742953647951029973386b0*%i,
1888  1b-32);
1889 true;
1890
1891 relerror(
1892   gamma_incomplete(0.5,5*%i),
1893   0.36441984106355895337750863822965 + 0.24368559063811288395048612967632*%i,
1894   3.0e-15);
1895 true;
1896
1897 relerror(
1898   gamma_incomplete(0.5,5b0*%i),
1899   0.36441984106355895337750863822965b0 + 0.24368559063811288395048612967632b0*%i,
1900   1b-32);
1901 true;
1902
1903 relerror(
1904   gamma_incomplete(0.5,50*%i),
1905   0.123398939447119035626603437814684 - 0.069012601491650179241716316049576*%i,
1906   1e-15);
1907 true;
1908
1909 relerror(
1910   gamma_incomplete(0.5,50b0*%i),
1911   0.123398939447119035626603437814684b0 - 0.069012601491650179241716316049576b0*%i,
1912   1b-32);
1913 true;
1914
1915 relerror(
1916   gamma_incomplete(0.5,100*%i),
1917   0.096899159215326861150517776631041 - 0.024684086404223368298902429118542*%i,
1918   1e-15);
1919 true;
1920
1921 relerror(
1922   gamma_incomplete(0.5,100b0*%i),
1923   0.096899159215326861150517776631041b0 - 0.024684086404223368298902429118542b0*%i,
1924   1b-32);
1925 true;
1926
1927 /* Along the boundary */
1928
1929 relerror(
1930   gamma_incomplete(0.5,-1+1.0*%i),
1931   -0.6460866463446816469727499577965 - 2.2529846180884648591828756025081*%i,
1932   1e-15);
1933 true;
1934
1935 relerror(
1936   gamma_incomplete(0.5,-2+2.0*%i),
1937   -4.6362149776912621823191772006210 - 0.6696404542328766212783095808086*%i,
1938   4.0e-15);
1939 true;
1940
1941 relerror(
1942   gamma_incomplete(0.5,-3+3.0*%i),
1943   -6.2748360165336721010541674791500 + 8.2442936688139837781281990469076*%i,
1944   4.0e-15);
1945 true;
1946
1947 relerror(
1948   gamma_incomplete(0.5,-4+4.0*%i),
1949   9.028176736353114550138088914592 + 22.441213878617781749976458917799*%i,
1950   5.0e-15);
1951 true;
1952
1953 relerror(
1954   gamma_incomplete(0.5,-5+5.0*%i),
1955   57.540711282025537715729492821303 + 9.813229347687737256352386117430*%i,
1956   3.0e-15);
1957 true;
1958
1959 relerror(
1960   gamma_incomplete(0.5,-6+6.0*%i),
1961   95.71761436932787302240057885539 - 107.54220609035943742021479062686*%i,
1962   2.0e-15);
1963 true;
1964
1965 relerror(
1966   gamma_incomplete(0.5,-10+10.0*%i),
1967   921.0687349037275816861912330294 + 5931.0759029700741177372337479958*%i,
1968   1.072e-15);
1969 true;
1970
1971 relerror(
1972   gamma_incomplete(0.5,-15+15.0*%i),
1973   -649132.21386404102825777002406965 + 315145.50427275409146151580354200*%i,
1974   1e-15);
1975 true;
1976
1977 relerror(
1978   gamma_incomplete(0.5,-1+1.0*%i),
1979   -0.6460866463446816469727499577965 - 2.2529846180884648591828756025081*%i,
1980   5.0e-15);
1981 true;
1982
1983 relerror(
1984   gamma_incomplete(0.5,-2+2.0*%i),
1985   -4.6362149776912621823191772006210 - 0.6696404542328766212783095808086*%i,
1986   4.0e-15);
1987 true;
1988
1989 relerror(
1990   gamma_incomplete(0.5,-3+3.0*%i),
1991   -6.2748360165336721010541674791500 + 8.2442936688139837781281990469076*%i,
1992   4.0e-15);
1993 true;
1994
1995 relerror(
1996   gamma_incomplete(0.5,-4+4.0*%i),
1997   9.028176736353114550138088914592 + 22.441213878617781749976458917799*%i,
1998   5.0e-15);
1999 true;
2000
2001 relerror(
2002   gamma_incomplete(0.5,-5+5.0*%i),
2003   57.540711282025537715729492821303 + 9.813229347687737256352386117430*%i,
2004   3.0e-15);
2005 true;
2006
2007 relerror(
2008   gamma_incomplete(0.5,-6+6.0*%i),
2009   95.71761436932787302240057885539 - 107.54220609035943742021479062686*%i,
2010   2.0e-15);
2011 true;
2012
2013 relerror(
2014   gamma_incomplete(0.5,-10+10.0*%i),
2015   921.0687349037275816861912330294 + 5931.0759029700741177372337479958*%i,
2016   1.072e-15);
2017 true;
2018
2019 relerror(
2020   gamma_incomplete(0.5,-15+15.0*%i),
2021   -649132.21386404102825777002406965 + 315145.50427275409146151580354200*%i,
2022   1e-15);
2023 true;
2024
2025 /******************************************************************************
2026   Test gamma_incomplete against expintegral_e
2027 ******************************************************************************/
2028
2029 block([badpoints : [],
2030        ratprint : false,
2031        abserr : 0,
2032        maxerr : -1,
2033        zlimit : 5,
2034        eps    : 1e-12],
2035   for a:1 thru 2 step 0.1 do
2036   (
2037     for z: -zlimit thru zlimit step 1.0 do
2038     (
2039       if is(notequal(z,0.0) and notequal(a,0.0)) then
2040       (
2041         zf : float(z),
2042         af : float(a),
2043         result : gamma_incomplete(af,zf),
2044         answer : rectform(zf^af*expintegral_e(1.0-af,zf)),
2045         abserr : abs(result - answer),
2046         maxerr : max(maxerr, abserr),
2047         if abserr > eps then
2048         (
2049           badpoints : cons([[a, z], result, answer, abserr], badpoints)
2050         )
2051       )
2052     )
2053   ),
2054   /*
2055    * For debugging, if there are any bad points, return the maximum error
2056    * found as the first element.
2057    */
2058   if badpoints # [] then
2059     cons(maxerr, badpoints)
2060   else
2061     badpoints
2062 );
2063 [];
2064
2065
2066 /******************************************************************************
2067   Test Generalized Incomplete Gamma function
2068 ******************************************************************************/
2069
2070 /* Some special values */
2071
2072 (kill(a), assume(a>0), done);
2073 done;
2074
2075 gamma_incomplete_generalized(a,z1,0);
2076 gamma_incomplete(a,z1)-gamma(a);
2077
2078 gamma_incomplete_generalized(a,z1,0.0);
2079 gamma_incomplete(a,z1)-gamma(a);
2080
2081 gamma_incomplete_generalized(a,z1,0.0b0);
2082 gamma_incomplete(a,z1)-gamma(a);
2083
2084 gamma_incomplete_generalized(a,0,z2);
2085 gamma(a)- gamma_incomplete(a,z2);
2086
2087 gamma_incomplete_generalized(a,0.0,z2);
2088 gamma(a)- gamma_incomplete(a,z2);
2089
2090 gamma_incomplete_generalized(a,0.0b0,z2);
2091 gamma(a)- gamma_incomplete(a,z2);
2092
2093 gamma_incomplete_generalized(a,z1,inf);
2094 gamma_incomplete(a,z1);
2095
2096 gamma_incomplete_generalized(a,inf,z2);
2097 -gamma_incomplete(a,z2);
2098
2099 gamma_incomplete_generalized(a,0,inf);
2100 gamma(a);
2101
2102 gamma_incomplete_generalized(a,x,x);
2103 0;
2104
2105 gamma_incomplete_generalized(a,1.0,1.0b0);
2106 0;
2107
2108 /* Mirror symmetry, but not when z1 or z2 on the negative real axis */
2109
2110 declare(a,complex, z1,complex, z2, complex);
2111 done;
2112
2113 conjugate(gamma_incomplete_generalized(a,z1,z2));
2114 conjugate(gamma_incomplete_generalized(a,z1,z2));
2115
2116 conjugate(gamma_incomplete_generalized(x+%i*y,1+%i,1-%i));
2117 gamma_incomplete_generalized(x-%i*y,1-%i,1+%i);
2118
2119 (kill(a,z1,z2),done);
2120 done;
2121
2122 /* Test numerical evaluation for some values */
2123
2124 closeto(
2125   gamma_incomplete_generalized(0.15,0.10,0.90),
2126   1.285210772938496575538196624140369253496313719924712338486508252b0,
2127   1e-14);
2128 true;
2129
2130 fpprec:64;
2131 64;
2132
2133 closeto(
2134   gamma_incomplete_generalized(0.15b0,0.10b0,0.90b0),
2135   1.285210772938496575538196624140369253496313719924712338486508252b0,
2136   1b-62);
2137 true;
2138
2139 closeto(
2140   gamma_incomplete_generalized(0.15+%i,0.10+%i,0.90+%i),
2141   (-0.03956290928621934869542750861441673192206453223955788892863857789b0
2142    -0.13316249485419500645510117515710482169661446536096647384481038655b0*%i),
2143   1e-14);
2144 true;
2145
2146 closeto(
2147   gamma_incomplete_generalized(0.15b0+%i,0.10b0+%i,0.90b0+%i),
2148   (-0.03956290928621934869542750861441673192206453223955788892863857789b0
2149    -0.13316249485419500645510117515710482169661446536096647384481038655b0*%i),
2150   1b-62);
2151 true;
2152
2153 closeto(
2154   gamma_incomplete_generalized(-0.15+%i,0.10+%i,0.90+%i),
2155   (-0.07903699552278027449948116754698066920498863638107044857029927559b0
2156    -0.10615739775378488990365404098165130400907362070260244159331987806b0*%i),
2157   1e-14);
2158 true;
2159
2160 closeto(
2161   gamma_incomplete_generalized(-0.15b0+%i,0.10b0+%i,0.90b0+%i),
2162   (-0.07903699552278027449948116754698066920498863638107044857029927559b0
2163    -0.10615739775378488990365404098165130400907362070260244159331987806b0*%i),
2164   1b-62);
2165 true;
2166
2167 /******************************************************************************
2168
2169   Test Regularized Incomplete Gamma function
2170
2171 ******************************************************************************/
2172
2173 /* Specialized values */
2174
2175 /* Check that gamma_incomplete_regularized (a, 0) just returns a noun
2176    form (the correct one!) */
2177 block ([result_a0: gamma_incomplete_regularized(a,0)],
2178   is (not (atom (result_a0)) and
2179       op (result_a0) = 'gamma_incomplete_regularized and
2180       args (result_a0) = [a, 0]));
2181 true$
2182
2183 (assume(ap>0,am<0),done);
2184 done;
2185
2186 errcatch(gamma_incomplete_regularized(am,0));
2187 [];
2188 gamma_incomplete_regularized(ap,0);
2189 1;
2190 gamma_incomplete_regularized(0,z);
2191 0;
2192 gamma_incomplete_regularized(a,inf);
2193 0;
2194
2195 /* Check that the derivative is correct. */
2196 diff(gamma_incomplete_regularized(a,z),z);
2197 -(z^(a-1)*%e^(-z))/gamma(a);
2198
2199 /* Expand gamma_incomplete_regularized */
2200
2201 gamma_expand:true;
2202 true;
2203
2204 gamma_incomplete_regularized(1,z);
2205 %e^(-z);
2206
2207 gamma_incomplete_regularized(a+1,z);
2208 gamma_incomplete_regularized(a,z)+%e^(-z)*z^a/gamma(a+1);
2209
2210 gamma_incomplete_regularized(a-1,z);
2211 gamma_incomplete_regularized(a,z)-%e^(-z)*z^(a-1)/gamma(a);
2212
2213 gamma_incomplete_regularized(1/2,z);
2214 erfc(sqrt(z));
2215 gamma_incomplete_regularized(-1/2,z);
2216 erfc(sqrt(z))-%e^(-z)/sqrt(%pi)/sqrt(z);
2217 gamma_incomplete_regularized(3/2,z);
2218 erfc(sqrt(z))+2*sqrt(z)*%e^(-z)/sqrt(%pi);
2219 gamma_incomplete_regularized(-3/2,z);
2220 erfc(sqrt(z))-3*(4*z/3-2/3)*%e^(-z)/(4*sqrt(%pi)*z^(3/2));
2221
2222 /* Test expansion for half integral values against expansion of
2223    gamma_incomplete
2224 */
2225
2226 expand(gamma_incomplete_regularized(5/2,z)-gamma_incomplete(5/2,z)/gamma(5/2));
2227 0;
2228 expand(gamma_incomplete_regularized(-5/2,z)-gamma_incomplete(-5/2,z)/gamma(-5/2));
2229 0;
2230 expand(gamma_incomplete_regularized(7/2,z)-gamma_incomplete(7/2,z)/gamma(7/2));
2231 0;
2232 expand(gamma_incomplete_regularized(-7/2,z)-gamma_incomplete(-7/2,z)/gamma(-7/2));
2233 0;
2234 expand(gamma_incomplete_regularized(9/2,z)-gamma_incomplete(9/2,z)/gamma(9/2));
2235 0;
2236 expand(gamma_incomplete_regularized(-9/2,z)-gamma_incomplete(-9/2,z)/gamma(-9/2));
2237 0;
2238
2239 gamma_expand:false;
2240 false;
2241
2242 /* Some numerical tests */
2243
2244 fpprec:64;
2245 64;
2246
2247 closeto(
2248   gamma_incomplete_regularized(0.25,0.15),
2249   0.3331718023153566353128831003164180886983644245472471410932121590b0,
2250   1e-15);
2251 true;
2252
2253 closeto(
2254   gamma_incomplete_regularized(0.25b0,0.15b0),
2255   0.3331718023153566353128831003164180886983644245472471410932121590b0,
2256   1b-62);
2257 true;
2258
2259 closeto(
2260   gamma_incomplete_regularized(-0.25,0.15),
2261    -0.3747953569677745583399657181155178573572870781780605755597341785b0,
2262   1e-15);
2263 true;
2264
2265 closeto(
2266   gamma_incomplete_regularized(-0.25b0,0.15b0),
2267    -0.3747953569677745583399657181155178573572870781780605755597341785b0,
2268   1b-61);
2269 true;
2270
2271 closeto(
2272   gamma_incomplete_regularized(-0.25,-0.15),
2273   (0.1206888313473692669850487605186163406801228067412581029970643626b0
2274   +0.8793111686526307330149512394813836593198771932587418970029356374b0*%i),
2275   1e-15);
2276 true;
2277
2278 closeto(
2279   gamma_incomplete_regularized(-0.25b0,-0.15b0),
2280   (0.1206888313473692669850487605186163406801228067412581029970643626b0
2281   +0.8793111686526307330149512394813836593198771932587418970029356374b0*%i),
2282   1e-15);
2283 true;
2284
2285 closeto(
2286   gamma_incomplete_regularized(0.25,0.15+%i),
2287   (-0.0241833450538703040924257417951024895368614341619005659619193200b0
2288    -0.1759768209797086273285777898669237251900625192446301525551431309b0*%i),
2289   1e-15);
2290 true;
2291
2292 closeto(
2293   gamma_incomplete_regularized(0.25b0,0.15b0+%i),
2294   (-0.0241833450538703040924257417951024895368614341619005659619193200b0
2295    -0.1759768209797086273285777898669237251900625192446301525551431309b0*%i),
2296   1b-62);
2297 true;
2298
2299 closeto(
2300   gamma_incomplete_regularized(0.25,0.15-%i),
2301   (-0.0241833450538703040924257417951024895368614341619005659619193200b0
2302    +0.1759768209797086273285777898669237251900625192446301525551431309b0*%i),
2303   1e-15);
2304 true;
2305
2306 closeto(
2307   gamma_incomplete_regularized(0.25b0,0.15b0-%i),
2308   (-0.0241833450538703040924257417951024895368614341619005659619193200b0
2309    +0.1759768209797086273285777898669237251900625192446301525551431309b0*%i),
2310   1b-62);
2311 true;
2312
2313 /******************************************************************************
2314   Test the Error functions: erf, erfc, erfi and erf_generalized
2315 ******************************************************************************/
2316
2317 /* Specific values */
2318
2319 map(erf, [0,0.0,0.0b0,inf,minf,infinity,und,ind]);
2320 [0,0.0,0.0b0,1,-1,erf(infinity),erf(und),erf(ind)];
2321
2322 limit(erf(z), z, inf);
2323 1;
2324 limit(erf(z), z, minf);
2325 -1;
2326
2327 map(erfc, [0,inf,minf,infinity,und,ind]);
2328 [1,0,2,erfc(infinity),erfc(und),erfc(ind)];
2329
2330 limit(erfc(z), z, inf);
2331 0;
2332 limit(erfc(z), z, minf);
2333 2;
2334
2335 map(erfi, [0,inf,minf,infinity,und,ind]);
2336 [0,inf,minf,erfi(infinity),erfi(und),erfi(ind)];
2337
2338 limit(erfi(z), z, inf);
2339 inf;
2340 limit(erfi(z), z, minf);
2341 minf;
2342
2343 erf_generalized(z1,0);
2344 -erf(z1);
2345
2346 erf_generalized(0,z2);
2347 erf(z2);
2348
2349 erf_generalized(0,0);
2350 0;
2351
2352 erf_generalized(0.0,0.0);
2353 0;
2354
2355 erf_generalized(0.0b0,0.0b0);
2356 0;
2357
2358 erf_generalized(z1,inf);
2359 1-erf(z1);
2360 limit(erf_generalized(z1, z2), z2, inf);
2361 1-erf(z1);
2362
2363 erf_generalized(z1,minf);
2364 -erf(z1)-1;
2365 limit(erf_generalized(z1, z2), z2, minf);
2366 -erf(z1)-1;
2367
2368 erf_generalized(inf,z2);
2369 erf(z2)-1;
2370 limit(erf_generalized(z1, z2), z1, inf);
2371 erf(z2)-1;
2372
2373 erf_generalized(minf,z2);
2374 erf(z2)+1;
2375 limit(erf_generalized(z1, z2), z1, minf);
2376 erf(z2)+1;
2377
2378 /* Parity */
2379
2380 erf(-z);
2381 -erf(z);
2382
2383 erfc(-z);
2384 2-erfc(z);
2385
2386 erfi(-z);
2387 -erfi(z);
2388
2389 erf_generalized(-z1,-z2);
2390 -erf_generalized(z1,z2);
2391
2392 /* Mirror symmetry */
2393
2394 declare(z,complex);
2395 done;
2396
2397 conjugate(erf(z));
2398 erf(conjugate(z));
2399
2400 conjugate(erf(x+%i*y));
2401 erf(x-%i*y);
2402
2403 conjugate(erfc(z));
2404 erfc(conjugate(z));
2405
2406 conjugate(erfc(x+%i*y));
2407 erfc(x-%i*y);
2408
2409 conjugate(erfi(z));
2410 erfi(conjugate(z));
2411
2412 conjugate(erfi(x+%i*y));
2413 erfi(x-%i*y);
2414
2415 declare(z1,complex,z2,complex);
2416 done;
2417
2418 conjugate(erf_generalized(z1,z2));
2419 erf_generalized(conjugate(z1),conjugate(z2));
2420
2421 conjugate(erf_generalized(x1+%i*y1,x2+%i*y2));
2422 erf_generalized(x1-%i*y1,x2-%i*y2);
2423
2424 /* Generalized Erf is antisymmetric */
2425
2426 erf_generalized(x1,x2)+erf_generalized(x2,x1);
2427 0;
2428
2429 /* For a pure real or imaginary argument of the error functions erf and erfi
2430    we get pure real or imaginary result. We test it. */
2431
2432 is(equal(imagpart(erf(1.0)),0));
2433 true;
2434
2435 is(equal(imagpart(erfi(1.0)),0));
2436 true;
2437
2438 is(equal(realpart(erf(1.0*%i)),0));
2439 true;
2440
2441 is(equal(realpart(erfi(1.0*%i)),0));
2442 true;
2443
2444 /* Again for bigfloats */
2445
2446 is(equal(imagpart(erf(1.0b0)),0));
2447 true;
2448
2449 is(equal(imagpart(erfi(1.0b0)),0));
2450 true;
2451
2452 is(equal(realpart(erf(1.0b0*%i)),0));
2453 true;
2454
2455 is(equal(realpart(erfi(1.0b0*%i)),0));
2456 true;
2457
2458 /* Taylor expansion of the Error functions */
2459
2460 taylor(erf(x),x,0,5) - 2/sqrt(%pi)*(x-x^3/3+x^5/10);
2461 0;
2462
2463 (erf(taylor(x,x,0,5)) - 2/sqrt(%pi)*(x-x^3/3+x^5/10));
2464 0;
2465
2466 taylor(erf(x),x,x0,2)
2467   - (erf(x0)+2*%e^(-x0^2)/sqrt(%pi)*(x-x0)-2*x0*%e^(-x0^2)/sqrt(%pi)*(x-x0)^2);
2468 0;
2469
2470 taylor(erfi(x),x,0,5) - 2/sqrt(%pi)*(x+x^3/3+x^5/10);
2471 0;
2472
2473 erfi(taylor(x,x,0,5)) - 2/sqrt(%pi)*(x+x^3/3+x^5/10);
2474 0;
2475
2476 taylor(erfi(x),x,x0,2)
2477  -(erfi(x0)+2*%e^(x0^2)/sqrt(%pi)*(x-x0)+2*x0*%e^(x0^2)/sqrt(%pi)*(x-x0)^2);
2478 0;
2479
2480 /* Numerical test for the Error functions
2481    First check erf in double float precision */
2482
2483 closeto(
2484   erf(-0.50),
2485   -0.5204998778130465376827466538919645287364515757579637000588057256b0,
2486   1e-15);
2487 true;
2488
2489 closeto(
2490   erf(0.50),
2491   0.5204998778130465376827466538919645287364515757579637000588057256b0,
2492   1e-15);
2493 true;
2494
2495 closeto(
2496   erf(0.75),
2497   0.7111556336535151315989378345914107773742059540965372322781333971b0,
2498   5.0e-15);
2499 true;
2500
2501 closeto(
2502   erf(-0.75+%i),
2503   (-1.372897192365736489613456241111589390954675856186764729607156305b0
2504   + 0.539788632227100129936591912063260716699852732091113612337142798b0*%i),
2505   6.0e-15);
2506 true;
2507
2508 closeto(
2509   erf(-0.50+%i),
2510   (-1.204847558314218002702112682097006717296399718277162764595960866b0
2511   + 1.024400881608445881724860454410886676966127193189877583791256132b0*%i),
2512   1e-15);
2513 true;
2514
2515 closeto(
2516   erf(0.50+%i),
2517   ( 1.204847558314218002702112682097006717296399718277162764595960866b0
2518   + 1.024400881608445881724860454410886676966127193189877583791256132b0*%i),
2519   1e-15);
2520 true;
2521
2522 closeto(
2523   erf(0.75+%i),
2524   ( 1.372897192365736489613456241111589390954675856186764729607156305b0
2525   + 0.539788632227100129936591912063260716699852732091113612337142798b0*%i),
2526   6.0e-15);
2527 true;
2528
2529 /* Bug 3587184: erf inaccurate for small float values */
2530 relerror(
2531   erf(1d-10),
2532   float(2d-10/sqrt(%pi)),
2533   1d-15);
2534 true;
2535
2536 /* Erf with bigfloat precision */
2537
2538 fpprec:64;
2539 64;
2540
2541 closeto(
2542   erf(-0.50b0),
2543   -0.5204998778130465376827466538919645287364515757579637000588057256b0,
2544   1b-61);
2545 true;
2546
2547 closeto(
2548   erf(0.50b0),
2549   0.5204998778130465376827466538919645287364515757579637000588057256b0,
2550   1b-61);
2551 true;
2552
2553 closeto(
2554   erf(0.75b0),
2555   0.7111556336535151315989378345914107773742059540965372322781333971b0,
2556   1b-61);
2557 true;
2558
2559 closeto(
2560   erf(-0.75b0+%i),
2561   (-1.372897192365736489613456241111589390954675856186764729607156305b0
2562   + 0.539788632227100129936591912063260716699852732091113612337142798b0*%i),
2563   1b-61);
2564 true;
2565
2566 closeto(
2567   erf(-0.50b0+%i),
2568   (-1.204847558314218002702112682097006717296399718277162764595960866b0
2569   + 1.024400881608445881724860454410886676966127193189877583791256132b0*%i),
2570   1b-61);
2571 true;
2572
2573 closeto(
2574   erf(0.50b0+%i),
2575   ( 1.204847558314218002702112682097006717296399718277162764595960866b0
2576   + 1.024400881608445881724860454410886676966127193189877583791256132b0*%i),
2577   1b-61);
2578 true;
2579
2580 closeto(
2581   erf(0.75b0+%i),
2582   ( 1.372897192365736489613456241111589390954675856186764729607156305b0
2583   + 0.539788632227100129936591912063260716699852732091113612337142798b0*%i),
2584   1b-61);
2585 true;
2586
2587 /* Bug 3587191: erf inaccurate for small bigfloat values */
2588 relerror(
2589   erf(1b-40),
2590   bfloat(2b-40/sqrt(%pi)),
2591   1b-65);
2592 true;
2593
2594 /* Bug 3587304 erfc(x) for x > 6 is wrong */
2595 relerror(
2596   erfc(6.0),
2597   2.15197367124989131659335039918738463047751406e-17,
2598   2.005e-15);
2599 true;
2600
2601 relerror(
2602   erfc(-4.0),
2603   2-erfc(4.0),
2604   1e-15);
2605 true;
2606
2607 relerror(
2608   erfc(6b0),
2609   2.1519736712498913116593350399187384630477514061688542100527892051056337238484927860b-17,
2610   1b-64);
2611 true;
2612
2613 relerror(
2614   erfc(-6b0),
2615   1.9999999999999999784802632875010868834066496008126153695224859383114578994b0,
2616   1b-64);
2617 true;
2618
2619 /* Bug 3587362 inverse_erfc(1d-40) wrong */
2620 relerror(
2621   inverse_erfc(1d-40),
2622   9.448789766720855,
2623   1d-15);
2624 true;
2625
2626 relerror(
2627   inverse_erfc(1b-50),
2628   10.59209016952736518902166392532979911559420645541709912588406440119671044289134079569127583320351428635b0,
2629   1d-64);
2630 true;
2631
2632 /* We have done a check for the numerical algorithm of the Erf function which
2633    calls the Incomplete Gamma function.
2634    We do not do further numerical tests for the other Error functions,
2635    but only check the correct implementation of the numercial routines
2636    against the Erf function
2637 */
2638
2639 /* Check Erfc against Erf for some values */
2640
2641 fpprec:32;
2642 32;
2643
2644 closeto(
2645   erfc(-0.25),
2646   1-erf(-0.25),
2647   1e-15);
2648 true;
2649
2650 closeto(
2651   erfc(0.25),
2652   1-erf(0.25),
2653   1e-15);
2654 true;
2655
2656 closeto(
2657   erfc(-0.25+%i),
2658   1-erf(-0.25+%i),
2659   1e-15);
2660 true;
2661
2662 closeto(
2663   erfc(0.25+%i),
2664   1-erf(0.25+%i),
2665   1e-15);
2666 true;
2667
2668 closeto(
2669   erfc(-0.25b0),
2670   1-erf(-0.25b0),
2671   1.0b-30);
2672 true;
2673
2674 closeto(
2675   erfc(0.25b0),
2676   1-erf(0.25b0),
2677   1b-30);
2678 true;
2679
2680 closeto(
2681   erfc(-0.25b0+%i),
2682   1-erf(-0.25b0+%i),
2683   1b-30);
2684 true;
2685
2686 closeto(
2687   erfc(0.25b0+%i),
2688   1-erf(0.25b0+%i),
2689   1b-30);
2690 true;
2691
2692 /* Check Erfi against Erf for some values */
2693
2694 closeto(
2695   erfi(-0.15),
2696   -%i*erf(-0.15*%i),
2697   1e-15);
2698 true;
2699
2700 closeto(
2701   erfi(0.15),
2702   -%i*erf(0.15*%i),
2703   1e-15);
2704 true;
2705
2706 closeto(
2707   erfi(-0.15+%i),
2708   -%i*erf((-0.15+%i)*%i),
2709   1e-15);
2710 true;
2711
2712 closeto(
2713   erfi(0.15+%i),
2714   -%i*erf((0.15+%i)*%i),
2715   1e-15);
2716 true;
2717
2718 closeto(
2719   erfi(-0.15b0),
2720   -%i*erf(-0.15b0*%i),
2721   5.0b-30);
2722 true;
2723
2724 closeto(
2725   erfi(0.15b0),
2726   -%i*erf(0.15b0*%i),
2727   5.0b-30);
2728 true;
2729
2730 closeto(
2731   erfi(-0.15b0+%i),
2732   expand(-%i*erf((-0.15b0+%i)*%i)),
2733   1b-30);
2734 true;
2735
2736 closeto(
2737   erfi(0.15b0+%i),
2738   expand(-%i*erf((0.15b0+%i)*%i)),
2739   1b-30);
2740 true;
2741
2742 /* Check Generalized Erf against Erf for some values */
2743
2744 closeto(
2745   erf_generalized(0.35,1.25),
2746   erf(1.25)-erf(0.35),
2747   1e-15);
2748 true;
2749
2750 closeto(
2751   erf_generalized(0.35+%i,1.25),
2752   erf(1.25)-erf(0.35+%i),
2753   1e-15);
2754 true;
2755
2756 closeto(
2757   erf_generalized(0.35,1.25+%i),
2758   erf(1.25+%i)-erf(0.35),
2759   1e-15);
2760 true;
2761
2762 closeto(
2763   erf_generalized(0.35+%i,1.25+%i),
2764   erf(1.25+%i)-erf(0.35+%i),
2765   1e-15);
2766 true;
2767
2768 closeto(
2769   erf_generalized(0.35b0,1.25b0),
2770   erf(1.25b0)-erf(0.35b0),
2771   1b-30);
2772 true;
2773
2774 closeto(
2775   erf_generalized(0.35b0+%i,1.25b0),
2776   erf(1.25b0)-erf(0.35b0+%i),
2777   1b-30);
2778 true;
2779
2780 closeto(
2781   erf_generalized(0.35b0,1.25b0+%i),
2782   erf(1.25b0+%i)-erf(0.35b0),
2783   1b-30);
2784 true;
2785
2786 closeto(
2787   erf_generalized(0.35b0+%i,1.25b0+%i),
2788   erf(1.25b0+%i)-erf(0.35b0+%i),
2789   1b-30);
2790 true;
2791
2792 /* Hypergeometric representations of the erf functions */
2793
2794 hypergeometric_representation:true;
2795 true;
2796
2797 erf(z);
2798 2*z*'hypergeometric([1/2],[3/2],-z^2)/sqrt(%pi);
2799
2800 erfc(z);
2801 1-2*z*'hypergeometric([1/2],[3/2],-z^2)/sqrt(%pi);
2802
2803 erfi(z);
2804 2*z*'hypergeometric([1/2],[3/2],z^2)/sqrt(%pi);
2805
2806 erf_generalized(z1,z2);
2807 2*z2*'hypergeometric([1/2],[3/2],-z2^2)/sqrt(%pi) -
2808 2*z1*'hypergeometric([1/2],[3/2],-z1^2)/sqrt(%pi);
2809
2810 hypergeometric_representation:false;
2811 false;
2812
2813 /******************************************************************************
2814    Test the Fresnel Integrals S(z) and C(z)
2815 ******************************************************************************/
2816
2817 /* Specific values for the Fresnel Integrals */
2818
2819 map(fresnel_s,[0,0.0,0.0b0]);
2820 [0,0.0,0.b0];
2821
2822 map(fresnel_c,[0,0.0,0.0b0]);
2823 [0,0.0,0.b0];
2824
2825 limit(fresnel_c(x),x,inf);
2826 1/2;
2827
2828 limit(fresnel_s(x),x,inf);
2829 1/2;
2830
2831 limit(fresnel_c(x),x,minf);
2832 -1/2;
2833
2834 limit(fresnel_s(x),x,minf);
2835 -1/2;
2836
2837 /* Simplification of infinities
2838    The rules for an odd function and the simplifaction for imaginary
2839    arguments are applied too.
2840 */
2841
2842 map(fresnel_s,[inf,-inf,minf,-minf,%i*inf,-%i*inf,%i*minf,-%i*minf]);
2843 [1/2,-1/2,-1/2,1/2,-%i/2,%i/2,%i/2,-%i/2];
2844
2845 map(fresnel_c,[inf,-inf,minf,-minf,%i*inf,-%i*inf,%i*minf,-%i*minf]);
2846 [1/2,-1/2,-1/2,1/2,%i/2,-%i/2,-%i/2,%i/2];
2847
2848 /* No simplification for other infinities and undeterminates */
2849
2850 map(fresnel_s,[infinity,und,ind]);
2851 [fresnel_s(infinity),fresnel_s(und),fresnel_s(ind)];
2852
2853 map(fresnel_c,[infinity,und,ind]);
2854 [fresnel_c(infinity),fresnel_c(und),fresnel_c(ind)];
2855
2856 /* The Fresnel Integrals S(z) and C(z) are odd functions
2857    A reflection rule is given and the rule odd-function-reflect is applied */
2858
2859 map(fresnel_s,[-x, (x-1), (-x+1), (-x-1)]);
2860 [-fresnel_s(x), fresnel_s(x-1), -fresnel_s(x-1),-fresnel_s(1+x)];
2861
2862 map(fresnel_c,[-x, (x-1), (-x+1), (-x-1)]);
2863 [-fresnel_c(x), fresnel_c(x-1), -fresnel_c(x-1),-fresnel_c(1+x)];
2864
2865 /* The Fresnel Integals simplify imaginary arguments */
2866
2867 map(fresnel_s, [%i,%i*x,-%i*x,%i*(x+1)]);
2868 [-%i*fresnel_s(1),-%i*fresnel_s(x),%i*fresnel_s(x),-%i*fresnel_s(x+1)];
2869
2870 map(fresnel_c, [%i,%i*x,-%i*x,%i*(x+1)]);
2871 [%i*fresnel_c(1),%i*fresnel_c(x),-%i*fresnel_c(x),%i*fresnel_c(x+1)];
2872
2873 /* The Fresnel Integrals have Mirror Symmetry */
2874
2875 declare(z,complex);
2876 done;
2877
2878 conjugate(fresnel_s(z));
2879 fresnel_s(conjugate(z));
2880
2881 conjugate(fresnel_s(x+%i*y));
2882 fresnel_s(x-%i*y);
2883
2884 conjugate(fresnel_c(z));
2885 fresnel_c(conjugate(z));
2886
2887 conjugate(fresnel_c(x+%i*y));
2888 fresnel_c(x-%i*y);
2889
2890 /* Taylor expansion of the Fresnel Integrals to order O(z^12) */
2891
2892 /* Expand the function */
2893 taylor(fresnel_s(x),x,0,12);
2894 %pi*x^3/6 - %pi^3*x^7/336 + %pi^5 *x^11/42240;
2895
2896 taylor(fresnel_c(x),x,0,12);
2897 x - %pi^2*x^5/40 + %pi^4 *x^9/3456;
2898
2899 /* Expand the argument and apply the function */
2900 fresnel_s(taylor(x,x,0,12));
2901 %pi*x^3/6 - %pi^3*x^7/336 + %pi^5 *x^11/42240;
2902
2903 fresnel_c(taylor(x,x,0,12));
2904 x - %pi^2*x^5/40 + %pi^4 *x^9/3456;
2905
2906 /* Differentiation of the Fresnel Integrals */
2907
2908 diff(fresnel_s(x),x);
2909 sin(%pi*x^2/2);
2910
2911 diff(fresnel_c(x),x);
2912 cos(%pi*x^2/2);
2913
2914 /* The elementary Integral of the Fresnel Integrals
2915    More complicated integrals can be found in rtest_integrate_special.mac */
2916
2917 integrate(fresnel_s(x),x);
2918 x*fresnel_s(x)+1/%pi*cos(%pi*x^2/2);
2919
2920 integrate(fresnel_c(x),x);
2921 x*fresnel_c(x)-1/%pi*sin(%pi*x^2/2);
2922
2923 /* Representation of the Fresnel Integrals through the Error function Erf */
2924
2925 erf_representation:true;
2926 true;
2927
2928 fresnel_s(x);
2929 (1+%i)/4*(erf((1+%i)/2*sqrt(%pi)*x) - %i* erf((1-%i)/2*sqrt(%pi)*x));
2930
2931 fresnel_c(x);
2932 (1-%i)/4*(erf((1+%i)/2*sqrt(%pi)*x) + %i* erf((1-%i)/2*sqrt(%pi)*x));
2933
2934 erf_representation:false;
2935 false;
2936
2937 /* Representation of the Fresnel Integrals through the
2938    Hypergeometric function */
2939
2940 hypergeometric_representation:true;
2941 true;
2942
2943 fresnel_s(x);
2944 %pi*x^3/6*hypergeometric([3/4],[3/2,7/4],-%pi^2*x^4/16);
2945
2946 fresnel_c(x);
2947 x*hypergeometric([1/4],[1/2,5/4],-%pi^2*x^4/16);
2948
2949 hypergeometric_representation:false;
2950 false;
2951
2952 /* Numerical tests for the Fresnel Integrals */
2953
2954 fpprec:64;
2955 64;
2956
2957 /* Tests for fresnel_s
2958    Real argument 0.25,0.50,1.0,2.0,5.0,10.0 */
2959
2960 relerror(
2961   fresnel_s(0.25),
2962   0.008175600235777755778102308866942774752486734698017086013976457144b0,
2963   1e-15);
2964 true;
2965
2966 relerror(
2967   fresnel_s(0.25b0),
2968   0.008175600235777755778102308866942774752486734698017086013976457144b0,
2969   1b-64);
2970 true;
2971
2972 relerror(
2973   fresnel_s(0.50),
2974   0.06473243285999927761148051223061476765072591849351249278758894565b0,
2975   1e-15);
2976 true;
2977
2978 relerror(
2979   fresnel_s(0.50b0),
2980   0.06473243285999927761148051223061476765072591849351249278758894565b0,
2981   1b-64);
2982 true;
2983
2984 relerror(
2985   fresnel_s(1.0),
2986   0.4382591473903547660767566966251526374937865724524165673344073263b0,
2987   1e-15);
2988 true;
2989
2990 relerror(
2991   fresnel_s(1.0b0),
2992   0.4382591473903547660767566966251526374937865724524165673344073263b0,
2993   2.0b-64);
2994 true;
2995
2996 relerror(
2997   fresnel_s(2.0),
2998   0.3434156783636982421953008159580684568865418122025247675792689204b0,
2999   1e-15);
3000 true;
3001
3002 relerror(
3003   fresnel_s(2.0b0),
3004   0.3434156783636982421953008159580684568865418122025247675792689204b0,
3005   1.5b-63);
3006 true;
3007
3008 relerror(
3009   fresnel_s(5.0),
3010   0.4991913819171168867519283804659916554084319970723881534101411152b0,
3011   1e-15);
3012 true;
3013
3014 relerror(
3015   fresnel_s(5.0b0),
3016   0.4991913819171168867519283804659916554084319970723881534101411152b0,
3017   1b-64);
3018 true;
3019
3020 relerror(
3021   fresnel_s(10.0),
3022   0.4681699785848822404033511108104469460538427245558302799270062272b0,
3023   1e-15);
3024 true;
3025
3026 relerror(
3027   fresnel_s(10.0b0),
3028   0.4681699785848822404033511108104469460538427245558302799270062272b0,
3029   1b-64);
3030 true;
3031
3032 /* Complex argument 0.25+%i,0.50+%i,1.0+%i,2.0+%i,5.0+%i */
3033
3034 relerror(
3035   fresnel_s(0.25+%i),
3036   -0.2762104914409824591766528447060750469693825567583676638192814447b0
3037   - 0.4331061708987372646968782339728450577465559834165609320175098398b0*%i,
3038   1e-15);
3039 true;
3040
3041 relerror(
3042   fresnel_s(0.25b0+%i),
3043   -0.2762104914409824591766528447060750469693825567583676638192814447b0
3044   - 0.4331061708987372646968782339728450577465559834165609320175098398b0*%i,
3045   1b-64);
3046 true;
3047
3048 relerror(
3049   fresnel_s(0.50+%i),
3050    -0.7169788451833594258616872412575572495423663980107873580716959051b0
3051   - 0.3393082523853171783825858066689424560943675693298848574923503884b0*%i,
3052   1e-15);
3053 true;
3054
3055 relerror(
3056   fresnel_s(0.50b0+%i),
3057    -0.7169788451833594258616872412575572495423663980107873580716959051b0
3058   - 0.3393082523853171783825858066689424560943675693298848574923503884b0*%i,
3059   1b-64);
3060 true;
3061
3062 relerror(
3063   fresnel_s(1.0+%i),
3064   -2.061888219194840468080716536685708600815908323737868052048638806b0
3065   + 2.061888219194840468080716536685708600815908323737868052048638806b0*%i,
3066   1e-15);
3067 true;
3068
3069 relerror(
3070   fresnel_s(1.0b0+%i),
3071   -2.061888219194840468080716536685708600815908323737868052048638806b0
3072   + 2.061888219194840468080716536685708600815908323737868052048638806b0*%i,
3073   2.0b-64);
3074 true;
3075
3076 relerror(
3077   fresnel_s(2.0+%i),
3078   -15.58775110440458732748278797797881643198730378904101846898562610b0
3079   - 36.72546488399143842838787627677917885752587065976755449373500438b0*%i,
3080   4.0e-13);
3081 true;
3082
3083 relerror(
3084   fresnel_s(2.0b0+%i),
3085   -15.58775110440458732748278797797881643198730378904101846898562610b0
3086   - 36.72546488399143842838787627677917885752587065976755449373500438b0*%i,
3087   4.0b-62);
3088 true;
3089
3090 relerror(
3091   fresnel_s(5.0+%i),
3092   -204452.5505063210590493330126425293361797143144299005524393297869b0
3093   + 38438.9093777759513198736581956757227475265101347572827364567521b0*%i,
3094   3.0e-14);
3095 true;
3096
3097 relerror(
3098   fresnel_s(5.0b0+%i),
3099   -204452.5505063210590493330126425293361797143144299005524393297869b0
3100   + 38438.9093777759513198736581956757227475265101347572827364567521b0*%i,
3101   4.0b-63);
3102 true;
3103
3104 /* Tests for fresnel_c
3105    Real argument 0.25,0.50,1.0,2.0,5.0,10.0 */
3106
3107 relerror(
3108   fresnel_c(0.25),
3109   0.2497591503565431834592215178008857243781399770549380697377810451b0,
3110   1e-15);
3111 true;
3112
3113 relerror(
3114   fresnel_c(0.25b0),
3115   0.2497591503565431834592215178008857243781399770549380697377810451b0,
3116   2.0b-64);
3117 true;
3118
3119 relerror(
3120   fresnel_c(0.50),
3121   0.4923442258714463928788436651566816377660951457715012532946526193b0,
3122   1e-15);
3123 true;
3124
3125 relerror(
3126   fresnel_c(0.50b0),
3127   0.4923442258714463928788436651566816377660951457715012532946526193b0,
3128   1b-64);
3129 true;
3130
3131 relerror(
3132   fresnel_c(1.0),
3133   0.7798934003768228294742064136526901366306257081363209601031335832b0,
3134   1e-15);
3135 true;
3136
3137 relerror(
3138   fresnel_c(1.0b0),
3139   0.7798934003768228294742064136526901366306257081363209601031335832b0,
3140   2.0b-64);
3141 true;
3142
3143 relerror(
3144   fresnel_c(2.0),
3145   0.4882534060753407545002235033572610376883671545092153829475964427b0,
3146   1e-15);
3147 true;
3148
3149 relerror(
3150   fresnel_c(2.0b0),
3151   0.4882534060753407545002235033572610376883671545092153829475964427b0,
3152   5.1b-64);
3153 true;
3154
3155 relerror(
3156   fresnel_c(5.0),
3157   0.5636311887040122311021074044130139641207537623099921078616593412b0,
3158   1e-15);
3159 true;
3160
3161 relerror(
3162   fresnel_c(5.0b0),
3163   0.5636311887040122311021074044130139641207537623099921078616593412b0,
3164   2.0b-64);
3165 true;
3166
3167 relerror(
3168   fresnel_c(10.0),
3169   0.4998986942055157236141518477356211143923468402262626572074674093b0,
3170   1.73e-15);
3171 true;
3172
3173 relerror(
3174   fresnel_c(10.0b0),
3175   0.4998986942055157236141518477356211143923468402262626572074674093b0,
3176   2.47b-64);
3177 true;
3178
3179 /* Complex argument 0.25+%i,0.50+%i,1.0+%i,2.0+%i,5.0+%i */
3180
3181 relerror(
3182   fresnel_c(0.25+%i),
3183   0.0097446563393801078320153856258158947458946246448139394089371651b0
3184   + 0.8830495953515316267486626148763161568578011755347662219418704558b0*%i,
3185   1e-15);
3186 true;
3187
3188 relerror(
3189   fresnel_c(0.25b0+%i),
3190   0.0097446563393801078320153856258158947458946246448139394089371651b0
3191   + 0.8830495953515316267486626148763161568578011755347662219418704558b0*%i,
3192   1b-64);
3193 true;
3194
3195 relerror(
3196   fresnel_c(0.50+%i),
3197   0.1199549363708813724035204126626172258713764410185526201803481040b0
3198   + 1.2468579809337107891237368150539206578767217615436096318514136378b0*%i,
3199   1e-15);
3200 true;
3201
3202 relerror(
3203   fresnel_c(0.50b0+%i),
3204   0.1199549363708813724035204126626172258713764410185526201803481040b0
3205   + 1.2468579809337107891237368150539206578767217615436096318514136378b0*%i,
3206   1b-64);
3207 true;
3208
3209 relerror(
3210   fresnel_c(1.0+%i),
3211   2.555793778102439024634522388352195842156623604203584296357752992b0
3212   + 2.555793778102439024634522388352195842156623604203584296357752992b0*%i,
3213   1e-15);
3214 true;
3215
3216 relerror(
3217   fresnel_c(1.0b0+%i),
3218   2.555793778102439024634522388352195842156623604203584296357752992b0
3219   + 2.555793778102439024634522388352195842156623604203584296357752992b0*%i,
3220   2.0b-64);
3221 true;
3222
3223 relerror(
3224   fresnel_c(2.0+%i),
3225   -36.22568799288165021578757830205090012554103292231420092205271528b0
3226   + 16.08787137412548041729488986757594802326683407694887086741793640b0*%i,
3227   4.0e-13);
3228 true;
3229
3230 relerror(
3231   fresnel_c(2.0b0+%i),
3232   -36.22568799288165021578757830205090012554103292231420092205271528b0
3233   + 16.08787137412548041729488986757594802326683407694887086741793640b0*%i,
3234   4.0b-62);
3235 true;
3236
3237 relerror(
3238   fresnel_c(5.0+%i),
3239   38439.4093777740143202918713550472184235160647415045418329908291b0
3240   + 204453.0505063119004499308897858846121088557663560705550579996093b0*%i,
3241   3.0e-14);
3242 true;
3243
3244 relerror(
3245   fresnel_c(5.0b0+%i),
3246   38439.4093777740143202918713550472184235160647415045418329908291b0
3247   + 204453.0505063119004499308897858846121088557663560705550579996093b0*%i,
3248   4.0b-63);
3249 true;
3250
3251 /* Bug #2509 fresnel_s incorrect for small values */
3252 relerror(
3253   fresnel_c(1d-20),
3254   1d-20,
3255   1d-15);
3256 true;
3257
3258 relerror(
3259   fresnel_c(1b-40),
3260   1b-40,
3261   5b-63);
3262 true;
3263
3264 relerror(
3265   fresnel_s(1d-20),
3266   float(%pi/6*(1d-20)^3),
3267   1d-15);
3268 true;
3269
3270 relerror(
3271   fresnel_s(1b-40),
3272   bfloat(%pi/6*(1b-40)^3),
3273   5b-63);
3274 true;
3275
3276 /******************************************************************************
3277    Test the Beta incomplete function
3278 ******************************************************************************/
3279
3280 /* Specialized values */
3281
3282 (assume(am<0,ap>0,b>0),done);
3283 done;
3284
3285 beta_incomplete(a,b,1);
3286 beta(a,b);
3287
3288 beta_incomplete(ap,b,0);
3289 0;
3290
3291 errcatch(beta_incomplete(am,b,0));
3292 [];
3293
3294 (forget(am<0, ap>0,b>0),done);
3295 done;
3296
3297 /* b is a positive integer */
3298
3299 beta_incomplete(a,1,z);
3300 z^a/a;
3301 beta_incomplete(a,2,z);
3302 (a*(1-z)+1)*z^a/(a*(a+1));
3303 beta_incomplete(1,1,z);
3304 z;
3305
3306 /* b is a positive integer.
3307    Check the handling of float and bigfloat numbers */
3308 fpprec:16;
3309 16;
3310
3311 beta_incomplete(1.0,1,z);
3312 1.0*z;
3313 /* Unfortunate, but we don't get exactly 1.0b0 */
3314 (closeto(beta_incomplete(1.0b0,1,z)/z,1b0,1b-15));
3315 true;
3316 beta_incomplete(1.0,1,1/2);
3317 0.5;
3318 beta_incomplete(1.0b0,1,1/2);
3319 0.5b0;
3320 beta_incomplete(1,1,1/2+%i);
3321 %i+1/2;
3322 beta_incomplete(1.0,1,1/2+%i);
3323 1.0*%i+0.5;
3324 beta_incomplete(1.0b0,1,1/2+%i);
3325 1.0b0*%i+0.5b0;
3326 beta_incomplete(1,2,1/2);
3327 3/8;
3328 beta_incomplete(1.0,2,1/2);
3329 0.375;
3330 beta_incomplete(1.0b0,2,1/2);
3331 0.375b0;
3332 beta_incomplete(1,2,1/2+%i);
3333 (3/2-%i)*(1/2+%i)/2;
3334 beta_incomplete(1.0,2,1/2+%i);
3335 0.875+0.5*%i;
3336 beta_incomplete(1.0b0,2,1/2+%i);
3337 0.875b0+0.5b0*%i;
3338
3339 /* We check the expansion for b positive against the numerical evaluation */
3340 closeto(
3341   float(beta_incomplete(1,1,1/2)) - beta_incomplete(1.0,1.0,0.5),
3342   0.0,
3343   1e-15);
3344 true;
3345 closeto(
3346   float(beta_incomplete(3/2,2,1/2)) - beta_incomplete(1.5,2.0,0.5),
3347   0.0,
3348   1e-15);
3349 true;
3350 closeto(
3351   float(beta_incomplete(2,3,1/2)) - beta_incomplete(2.0,3.0,0.5),
3352   0.0,
3353   1e-15);
3354 true;
3355 closeto(
3356   float(beta_incomplete(5/2,4,1/2)) - beta_incomplete(2.5,4.0,0.5),
3357   0.0,
3358   1e-15);
3359 true;
3360 closeto(
3361   float(beta_incomplete(3,5,1/2)) - beta_incomplete(3.0,5.0,0.5),
3362   0.0,
3363   1e-15);
3364 true;
3365 closeto(
3366   float(beta_incomplete(2,6,1/2+%i)) - beta_incomplete(2.0,6.0,0.5+%i),
3367   0.0,
3368   1e-15);
3369 true;
3370 closeto(
3371   /* We do this extra rectform, because of a bug in abs for expressions
3372      with a complex exponent */
3373   rectform(float(beta_incomplete(2+%i,7,1/2+%i))
3374                - beta_incomplete(2.0+%i,7.0,0.5+%i)),
3375   0.0,
3376   1e-15);
3377 true;
3378
3379 /* a is a positive integer we expand */
3380 beta_incomplete(1,b,z);
3381 (1-(1-z)^b)/b;
3382 beta_incomplete(2,b,z);
3383  (1-(1-z)^b*(b*z+1))/(b*(b+1));
3384 beta_incomplete(3,b,z);
3385  2*(1-(1-z)^b*(b*(b+1)*z^2/2+b*z+1))/(b*(b+1)*(b+2));
3386
3387 /* a is a positive integer.
3388    Check the expansion against numerically evaluation.
3389    First we test for b a positive value. */
3390
3391 closeto(
3392   float(beta_incomplete(1,2,3/2))-beta_incomplete(1.0,2.0,1.5),
3393   0.0,
3394   1e-15);
3395 true;
3396 closeto(
3397   float(beta_incomplete(2,5/2,3/2))-beta_incomplete(2.0,2.5,1.5),
3398   0.0,
3399   1e-15);
3400 true;
3401 closeto(
3402   float(beta_incomplete(3,3,3/2))-beta_incomplete(3.0,3.0,1.5),
3403   0.0,
3404   1e-15);
3405 true;
3406 closeto(
3407   float(beta_incomplete(2,7/2,3/2))-beta_incomplete(2.0,3.5,1.5),
3408   0.0,
3409   1e-15);
3410 true;
3411 closeto(
3412   float(beta_incomplete(2,5/2,3/2+%i))-beta_incomplete(2.0,2.5,1.5+%i),
3413   0.0,
3414   1e-15);
3415 true;
3416 /* I (rtoy) think the limiting accuracy of the following test is the conversion of
3417  * the exact answer to floating point.  If the exact answer is converted to a
3418  * bfloat 32 digits, the difference is less than 1.6e-15 or so.
3419 */
3420 closeto(
3421   float(rectform(float(beta_incomplete(2,5/2,-3/2+%i))))
3422                 -beta_incomplete(2.0,2.5,-1.5+%i),
3423   0.0,
3424   4.35e-15); /* we have lost accuracy. More tests necessary? */
3425 true;
3426
3427 /* Incomplete Beta is defined for negative integers a and b >= (-a)
3428    Add this point we have a problem. functions.wolfram.com gives different
3429    numerical results for this cases. When b not equal -a Maxima and
3430    functions.wolfram.com differ by a factor 2. For other valid integer b
3431    functions.wolfram.com gives ComplexInfinity and not an expected result.
3432    What is wrong?
3433 */
3434
3435 beta_incomplete(-1,1,z);
3436 -1/z;
3437 beta_incomplete(-2,1,z);
3438 -1/(2*z^2);
3439 beta_incomplete(-2,2,z);
3440 (z-1/2)/z^2;
3441 beta_incomplete(-3,1,z);
3442 -1/(3*z^3);
3443 beta_incomplete(-3,2,z);
3444 (z/2-1/3)/z^3;
3445 beta_incomplete(-3,3,z);
3446 (-z^2+z-1/3)/z^3;
3447
3448 /* Some numerical tests in double float precision */
3449
3450 closeto(
3451   beta_incomplete(0.5,1.0,0.10),
3452   0.6324555320336758663997787088865437067439110278650433653715009706b0,
3453   1e-15);
3454 true;
3455
3456 closeto(
3457   beta_incomplete(0.5,1.0,0.15),
3458   0.7745966692414833770358530799564799221665843410583181653175147532b0,
3459   1e-15);
3460 true;
3461
3462 closeto(
3463   beta_incomplete(0.5,1.0,0.20),
3464   0.8944271909999158785636694674925104941762473438446102897083588982b0,
3465   1e-15);
3466 true;
3467
3468 closeto(
3469   beta_incomplete(0.5,1.0,0.25),
3470   1.000000000000000000000000000000000000000000000000000000000000000b0,
3471   1e-15);
3472 true;
3473
3474 closeto(
3475   beta_incomplete(0.5,1.0,0.50),
3476   1.414213562373095048801688724209698078569671875376948073176679738b0,
3477   1.416e-15);
3478 true;
3479
3480 closeto(
3481   beta_incomplete(0.5,1.0,0.75),
3482   1.732050807568877293527446341505872366942805253810380628055806979b0,
3483   2.11e-15);
3484 true;
3485
3486 closeto(
3487   beta_incomplete(0.5,1.0,1.00),
3488   2.000000000000000000000000000000000000000000000000000000000000000b0,
3489   2.0e-15);
3490 true;
3491
3492 closeto(
3493   beta_incomplete(0.5,1.0,1.25),
3494   2.236067977499789696409173668731276235440618359611525724270897245b0,
3495   2.0e-15);
3496 true;
3497
3498 closeto(
3499   beta_incomplete(0.5,1.0,1.50),
3500   2.449489742783178098197284074705891391965947480656670128432692567b0,
3501   2.0e-15);
3502 true;
3503
3504 closeto(
3505   beta_incomplete(0.5,1.0,1.75),
3506   2.645751311064590590501615753639260425710259183082450180368334459b0,
3507   1.446e-15);
3508 true;
3509
3510 closeto(
3511   beta_incomplete(0.5,1.0,2.00),
3512   2.828427124746190097603377448419396157139343750753896146353359476b0,
3513   1.503e-15);
3514 true;
3515
3516 /* Some numerical tests in bigfloat precision */
3517
3518 fpprec:64;
3519 64;
3520
3521 closeto(
3522   beta_incomplete(0.5,1.0,0.10b0),
3523   0.6324555320336758663997787088865437067439110278650433653715009706b0,
3524   1b-60);
3525 true;
3526
3527 closeto(
3528   beta_incomplete(0.5,1.0,0.15b0),
3529   0.7745966692414833770358530799564799221665843410583181653175147532b0,
3530   1b-60);
3531 true;
3532
3533 closeto(
3534   beta_incomplete(0.5,1.0,0.20b0),
3535   0.8944271909999158785636694674925104941762473438446102897083588982b0,
3536   1b-60);
3537 true;
3538
3539 closeto(
3540   beta_incomplete(0.5,1.0,0.25b0),
3541   1.000000000000000000000000000000000000000000000000000000000000000b0,
3542   1b-60);
3543 true;
3544
3545 closeto(
3546   beta_incomplete(0.5,1.0,0.50b0),
3547   1.414213562373095048801688724209698078569671875376948073176679738b0,
3548   1b-60);
3549 true;
3550
3551 closeto(
3552   beta_incomplete(0.5,1.0,0.75b0),
3553   1.732050807568877293527446341505872366942805253810380628055806979b0,
3554   1b-16);
3555 true;
3556
3557 closeto(
3558   beta_incomplete(0.5,1.0,1.00b0),
3559   2.000000000000000000000000000000000000000000000000000000000000000b0,
3560   1b-60);
3561 true;
3562
3563 closeto(
3564   beta_incomplete(0.5,1.0,1.25b0),
3565   2.236067977499789696409173668731276235440618359611525724270897245b0,
3566   1b-60);
3567 true;
3568
3569 closeto(
3570   beta_incomplete(0.5,1.0,1.50b0),
3571   2.449489742783178098197284074705891391965947480656670128432692567b0,
3572   1b-60);
3573 true;
3574
3575 closeto(
3576   beta_incomplete(0.5,1.0,1.75b0),
3577   2.645751311064590590501615753639260425710259183082450180368334459b0,
3578   1b-60);
3579 true;
3580
3581 closeto(
3582   beta_incomplete(0.5,1.0,2.00b0),
3583   2.828427124746190097603377448419396157139343750753896146353359476b0,
3584   1b-60);
3585 true;
3586
3587 /* See Bug 3220128, but this isn't really that bug */
3588 closeto(
3589   gamma_incomplete(0, 200b0*%i),
3590   0.00437844609302782567916569771749325524128345091344187598851110680706344144459295b0
3591   - %i*.00241398745542678587253611621620491057595401709907514761094360488114169654741b0,
3592   3.4b-67);
3593 true;
3594
3595 /* Make sure we don't overflow unnecessarily in gamma_incomplete_regularized */
3596 closeto(
3597   gamma_incomplete_regularized(45001d0, 45000d0),
3598   0.5012537510700691773177801688515861486945632498553288,
3599   6.4e-12);
3600 true;
3601
3602 /* Check accuracy on the negative real axis.
3603  * See Bug  3526359 - gamma_incomplete(1/5,-32.0) not accurate
3604  */
3605 relerror(
3606   gamma_incomplete(1/5,-32.0),
3607   -4.0986398326716649399d12 - %i*2.9778361454160762231d12,
3608   2.545e-15);
3609 true;
3610
3611 relerror(
3612   gamma_incomplete(10,-100d0),
3613   subst(x=-100d0, block([gamma_expand : true], gamma_incomplete(10,x))),
3614   8.78e-15);
3615 true;
3616
3617 relerror(
3618   gamma_incomplete(10,-100b0),
3619   subst(x=-100b0, block([gamma_expand : true], gamma_incomplete(10,x))),
3620   3.57b-64);
3621 true;
3622
3623 relerror(
3624   gamma_incomplete(10,-100d0+%i),
3625   subst(x=-100d0+%i, block([gamma_expand : true], gamma_incomplete(10,x))),
3626   1.046e-14);
3627 true;
3628
3629 relerror(
3630   gamma_incomplete(10,-100b0+%i),
3631   subst(x=-100b0+%i, block([gamma_expand : true], gamma_incomplete(10,x))),
3632   2.04b-64);
3633 true;
3634
3635 relerror(
3636  erf(inverse_erf(.5)),
3637  0.5,
3638  1.4d-14);
3639 true;
3640
3641 relerror(
3642  erf(inverse_erf(-.5)),
3643  -0.5,
3644  1.4e-14);
3645 true;
3646
3647 relerror(
3648  erf(inverse_erf(.5b0)),
3649  0.5b0,
3650  2b-64);
3651 true;
3652
3653 relerror(
3654  erf(inverse_erf(-.5b0)),
3655  -0.5b0,
3656  2b-64);
3657 true;
3658
3659 relerror(
3660  erf(inverse_erf(2.0)),
3661  2.0,
3662  7.9d-15);
3663 true;
3664
3665 relerror(
3666  erf(inverse_erf(-2.0)),
3667  -2.0,
3668  7.9d-15);
3669 true;
3670
3671 /* These are a bit slow if fpprec is large.  Make fpprec a bit smaller */
3672 fpprec:24;
3673 24;
3674
3675 relerror(
3676  erf(inverse_erf(2b0)),
3677  2b0,
3678  5.9b-24);
3679 true;
3680
3681 relerror(
3682  erf(inverse_erf(-2b0)),
3683  -2b0,
3684  5.9b-24);
3685 true;
3686
3687 relerror(
3688  erf(inverse_erf(2b0+2b0*%i)),
3689  2b0+2b0*%i,
3690  6.6b-25);
3691 true;
3692
3693 relerror(
3694  erf(inverse_erf(-2b0-2b0*%i)),
3695  -2b0-2b0*%i,
3696  6.6b-25);
3697 true;
3698
3699 relerror(
3700  erf(inverse_erf(10b0+1000b0*%i)),
3701  10b0+1000b0*%i,
3702  2b-23);
3703 true;
3704
3705
3706 /* From the mailing list.  This should not signal an error and should give a simplifed answer */
3707 relerror(
3708  expand(bfloat(erf((sqrt(2)*%i-sqrt(2))/2))),
3709  4.74147636640994245161681b-1 * %i - 9.69264211944215930381491b-1,
3710  1b-23);
3711 true;
3712
3713 /* Bug #2619
3714  * Function inverse_erf - error in numerical evaluation
3715  */
3716 relerror(inverse_erf(0.7715),
3717  .8515204685911937,
3718  1d-15);
3719 true;
3720
3721 relerror(inverse_erf(- .9763545580611441 ),
3722  -1.600070795398459,
3723  1d-15);
3724 true;
3725
3726 /* Bug #2668 Bigfloat Gamma inaccurae for small inputs
3727  *
3728  * For these small values, gamma(z) = gamma(z+1)/z = 1/z
3729  * since 1+z = 1 and gamma(1) = 1.
3730  */
3731 relerror(gamma(2.0^-256), 2.0^256, 1d-14);
3732 true;
3733
3734 relerror(gamma(2b0^-256), 2b0^256, 10^(-fpprec+1));
3735 true;
3736
3737 relerror(
3738   gamma(2.0^-256 + 1d-200*%i),
3739   rectform(1/(2.0^-256 + 1d-200*%i)),
3740   2d-15);
3741 true;
3742
3743 closeto(
3744   gamma(2b0^-256 + 1b-500*%i),
3745   rectform(1/(2b0^-256 + 1b-500*%i)),
3746   2d-15);
3747 true;
3748
3749 (fpprec:oldfpprec,done);
3750 done;
3751
3752 /* SF bug #3090: "erfi switches sign at approximately -0.476 and 0.476" */
3753
3754 block
3755   (kill (x, L1, L2),
3756    x : makelist (i/10.0, i, -20, 20),
3757    L1 : makelist (erfi (x1), x1, x),
3758    L2 : makelist (-%i*erf (%i*x1), x1, x),
3759    every (lambda ([x1, x2], x1 = x2), L1, L2));
3760 true;
3761
3762 every (lambda ([x1, e1], sign(x1) = sign(e1)), x, L1);
3763 true;
3764
3765 every (lambda ([x1, e1], sign(x1) = sign(e1)), x, L2);
3766 true;
3767
3768 /* SF bug #3277: "no numerical evaluation for gamma_greek" */
3769
3770 block([gamma_expand: true], gamma_incomplete_lower(10, 1));
3771 362880-986410*%e^-1;
3772
3773 /* numerical value not very accurate; perhaps it is being computed naively via preceding formula? */
3774 closeto (ev (gamma_incomplete_lower(10, 1), numer), 4.043407757955496e-2, 1e-10);
3775 true;
3776
3777 /* this one is way off the mark ... sheesh. */
3778 closeto (gamma_incomplete_lower(10, 1.0), 4.043407757955496e-2, 1e-8);
3779 true;
3780
3781 (reset (fpprec),
3782  closeto (gamma_incomplete_lower(10, 1b0), 4.043407757955496b-2, 1b-10));
3783 true;
3784
3785 block([gamma_expand: true], gamma_incomplete_lower(1, 10));
3786 1 - %e^-10;
3787
3788 closeto (ev (gamma_incomplete_lower(1, 10), numer), 0.9999546000702375e0, 1e-15);
3789 true;
3790
3791 closeto (gamma_incomplete_lower(1.0, 10), 0.9999546000702375e0, 1e-15);
3792 true;
3793
3794 closeto (gamma_incomplete_lower(1b0, 10), 0.9999546000702375b0, 1b-15);
3795 true;
3796
3797 /* complex float and complex bfloat */
3798
3799 /* yikes */
3800 closeto (ev (gamma_incomplete_lower(10, 1 + %i), numer), 7.992250288513467e-1*%i + 1.011121671734864e0, 1e-7);
3801 true;
3802
3803 /* ditto */
3804 closeto (gamma_incomplete_lower(10, 1.0 + %i), 7.992250288513467e-1*%i + 1.011121671734864e0, 1e-7);
3805 true;
3806
3807 closeto (gamma_incomplete_lower(10, 1b0 + %i), 7.992250288513467b-1*%i + 1.011121671734864b0, 1b-10);
3808 true;
3809
3810 /* mailing list 2018-02-04: "Bug in beta_incomplete_regularized." */
3811
3812 (kill(f, g, x, y),
3813  f(x) := beta_incomplete_regularized(2,2,x/(1+x)),
3814  g(x) := beta_incomplete(2,2,x),
3815  should_be_equal(a, b) := if equal(a, b) then true else notequal(a, b),
3816  0);
3817 0;
3818
3819 f(1/y);
3820 (2*(1-1/((1/y+1)*y))+1)/((1/y+1)^2*y^2);
3821
3822 f(1/x);
3823 (2*(1-1/((1/x+1)*x))+1)/((1/x+1)^2*x^2);
3824
3825 should_be_equal (f(1/x), subst (y = 1/x, f(y)));
3826 true;
3827
3828 g(-y);
3829 (y^2*(2*(y+1)+1))/6;
3830
3831 g(-x);
3832 (x^2*(2*(x+1)+1))/6;
3833
3834 should_be_equal (g(-x), subst (y = -x, g(y)));
3835 true;
3836
3837 (kill(a, h, i, j, k, u, v, w),
3838  h(u) := gamma_incomplete(3, u),
3839  i(v, w) := gamma_incomplete_generalized(3 + a, v, w),
3840  j(u) := gamma_incomplete_regularized(4, u),
3841  k(v, w) := beta_incomplete_regularized(2, v, w),
3842  0);
3843 0;
3844
3845 gamma_expand : true;
3846 true;
3847
3848 should_be_equal (h(u), integrate (t^2*exp(-t),t,u,inf));
3849 true;
3850
3851 [h(u), h(-u)];
3852 [2*(u^2/2+u+1)*%e^-u,2*(u^2/2-u+1)*%e^u];
3853
3854 should_be_equal (h(-u), subst (v = -u, h(v)));
3855 true;
3856
3857 should_be_equal (i(2*v, 4*w), a*(a+1)*(a+2)
3858                                *((-((4^(a+2)*w^(a+2))/(a*(a+1)*(a+2))+(4^(a+1)*w^(a+1))/(a*(a+1))
3859                                                                      +(4^a*w^a)/a)
3860                                 *%e^-(4*w))
3861                                 +((2^(a+2)*v^(a+2))/(a*(a+1)*(a+2))+(2^(a+1)*v^(a+1))/(a*(a+1))+(2^a*v^a)/a)
3862                                  *%e^-(2*v)+gamma_incomplete_generalized(a,2*v,4*w))); /* FOR SURE ?? */
3863 true;
3864
3865 should_be_equal (i(2*v, 4*w), subst ([x = 2*v, y = 4*w], i(x, y)));
3866 true;
3867
3868 should_be_equal (j(u), integrate (t^3*exp(-t),t,u,inf)/gamma(4));
3869 true;
3870
3871 should_be_equal (j(1 + u), ((u+1)^3/6+(u+1)^2/2+u+2)*%e^((-u)-1));
3872 true;
3873
3874 should_be_equal (j(1 + u), subst (v = 1 + u, j(v)));
3875 true;
3876
3877 block ([foo, bar, l:[0.1, 0.25, 0.5, 0.8, 0.9, 0.95]],
3878        foo : makelist (k(4, x), x, l),
3879        bar : ev (makelist (first (quad_qags (t^(a-1)*(1-t)^(b-1), t, 0, x))/beta(a, b), x, l), a=2, b=4),
3880        if lmax (abs (foo - bar)) < 1e-8 then true else notequal (foo, bar));
3881 true;
3882
3883 k(1 - w, 2*v);
3884 1-(1-2*v)^(1-w)*(2*v*(1-w)+1); /* assuming k(v, w) = 1-(1-w)^v*(v*w+1) */
3885
3886 should_be_equal (k(1 - v, 2*w), subst ([x = 1 - v, y = 2*w], k(x, y)));
3887 true;
3888
3889 /*
3890  * From the mailing list, 2020/01/20, Richard Fateman.  The integrand here is
3891  * the derivative of exp(-x^3)/x^4.  Without gamma_expand, we get an expression
3892  * involving gamma_incomplete of orders -1/3 and -4/3.  By having
3893  * gamma_expand:true, these are expressed in terms of order 2/3.  These nicely
3894  * cancel to produce the original expression.
3895  */
3896 expand(integrate((-(3*%e^-x^3)/x^2)-(4*%e^-x^3)/x^5, x));
3897 exp(-x^3)/x^4;
3898
3899 /*
3900  * Test expansion of gamma_incomplete_lower
3901  */
3902 gamma_incomplete_lower(a+2,z);
3903 a*(a+1)*gamma_incomplete_lower(a,z)-z^a*(z+a+1)*%e^-z;
3904
3905 gamma_incomplete_lower(1/3+2,z);
3906 (4*gamma_incomplete_lower(1/3,z))/9-z^(1/3)*(z+4/3)*%e^-z;
3907
3908 gamma_incomplete_lower(a-2,z);
3909 (1/((a-2)*z)+1/((a-2)*(a-1)))*z^(a-1)*%e^-z +gamma_incomplete_lower(a,z)/((a-2)*(a-1));
3910
3911 gamma_incomplete_lower(1/3-2,z);
3912 ((9/10-3/(5*z))*%e^-z)/z^(2/3)+(9*''gamma_incomplete_lower(1/3,z))/10;
3913
3914 /*
3915  * Test expansion of gamma_incomplete_regularized.
3916  */
3917 gamma_incomplete_regularized(3,z) - gamma_incomplete(3,z)/gamma(3);
3918 0;
3919
3920 gamma_incomplete_regularized(3/2,z);
3921 (2*sqrt(z)*%e^-z)/sqrt(%pi)+erfc(sqrt(z));
3922
3923 gamma_incomplete_regularized(-3/2,z);
3924 erfc(sqrt(z))-(3*((4*z)/3-2/3)*%e^-z)/(4*sqrt(%pi)*z^(3/2));
3925
3926 gamma_incomplete_regularized(2+a,z);
3927 (z^a*(z+a+1)*%e^-z)/(a*(a+1)*gamma(a))+gamma_incomplete_regularized(a,z);
3928
3929 gamma_incomplete_regularized(-2+a,z);
3930 gamma_incomplete_regularized(a,z)-(z^(a-2)*(z+a-1)*%e^-z)/gamma(a);
3931
3932 gamma_incomplete_regularized(2+1/3,z);
3933 (9*z^(1/3)*(z+4/3)*%e^-z)/(4*gamma(1/3))+gamma_incomplete_regularized(1/3,z);
3934
3935 gamma_incomplete_regularized(-2+1/3,z);
3936 gamma_incomplete_regularized(1/3,z)-((z-2/3)*%e^-z)/(gamma(1/3)*z^(5/3));
3937
3938 integrate(log_gamma(z),z);
3939 psi[-2](z)$
3940
3941 integrate(log_gamma(1-s),s);
3942 -psi[-2](1-s)$
3943
3944 integrate(2*x*log_gamma(x^2),x);
3945 psi[-2](x^2)$
3946
3947 reset (gamma_expand);
3948 [gamma_expand];
3949
3950 /* tests for derivative of binomial */
3951 ratsimp(makegamma(diff(binomial(a,b) - makegamma(binomial(a,b)),a)));
3952 0$
3953
3954 ratsimp(makegamma(diff(binomial(a,b) - makegamma(binomial(a,b)),b)));
3955 0$
3956
3957 ratsimp(makegamma(diff(binomial(a^2,a) - makegamma(binomial(a^2,a)),a,2)));
3958 0$
3959
3960 /* tests for sign of gamma expressions */
3961 sign(gamma(exp(x)));
3962 pos$
3963
3964 /* since x might be zero, making gamma(x) undefined, return pnz */
3965 sign(gamma(x));
3966 pnz$
3967
3968 /* since sign(1/x) = pn, we have sign(gamma(1/x)) = pn */
3969 sign(gamma(1/x));
3970 pn$
3971
3972 sign(gamma(x^2 +1));
3973 pos$
3974
3975 sign(gamma(-(1 + x^2)));
3976 pn$
3977
3978 /* since sign(sqrt(x)) = pz, we have sign(gamma(sqrt(x)) = pnz */
3979 sign(gamma(sqrt(x)));
3980 pnz$
3981
3982 /* tests for errors noted while reviewing
3983  * SF bug #2999: "bug in cdf_binomial caused by beta_incomplete_regularized"
3984  * problem was that beta_incomplete_regularized and beta_incomplete_generalized
3985  * are no longer implemented as Lisp functions (i.e., not FBOUNDP).
3986  * I guess they're implemented solely via simplification now.
3987  *
3988  * NOTE: expected results look plausible to me,
3989  * but I did not check them carefully.
3990  */
3991
3992 /* examples from #2999 */
3993
3994 beta_incomplete_regularized(4830,171,0.8333333333333334);
3995 1.134595788994921E-194;
3996
3997 beta_incomplete_regularized (1000.0, 1000.0, 0.5);
3998 0.5;
3999
4000 /* other cases discovered while fixing preceding error */
4001
4002 beta_incomplete_regularized (m + 2, n, 1/2), beta_expand = true;
4003 beta_incomplete_regularized(m,n,1/2)
4004  -((n+m)*((((-m)-1)*2^((-n)-m))/((-n)-m)+2^((-n)-m-1)))/(m*(m+1)*beta(m,n));
4005
4006 beta_incomplete_regularized (m - 2, n, 1/2), beta_expand = true;
4007 beta_incomplete_regularized(m,n,1/2)
4008  -(((1-m)*2^((-n)-m+2))/((-n)-m+2)+2^((-n)-m+1))/(beta(m,n)*(n+m-1));
4009
4010 beta_incomplete_generalized (m + 2, n, x, y), beta_expand = true;
4011 ((-(1-y)^n*y^(m+1))+(((-m)-1)*((1-x)^n*x^m-(1-y)^n*y^m))/((-n)-m)
4012                          +(1-x)^n*x^(m+1))
4013  /(n+m+1)
4014  +(m*(m+1)*beta_incomplete_generalized(m,n,x,y))/((n+m)*(n+m+1));
4015
4016 beta_incomplete_generalized (m - 2, n, x, y), beta_expand = true;
4017 (beta_incomplete_generalized(m,n,x,y)*((-n)-m+1)*((-n)-m+2))/((1-m)
4018                                                                    *(2-m))
4019  -(((-n)-m+2)*((1-y)^n*y^(m-1)+((1-m)*((1-y)^n*y^(m-2)-(1-x)^n*x^(m-2)))
4020                                /((-n)-m+2)-(1-x)^n*x^(m-1)))
4021   /((1-m)*(2-m));
4022
4023 /* verify derivatives of beta_incomplete and friends (unreported bug) */
4024
4025 (kill (a, b, x, y, e),
4026 finite_difference_epsilon: 1e-8,
4027
4028 finite_difference_jacobian (ee, vars, vars0) :=
4029     genmatrix (lambda ([i, j], finite_difference_partial (ee[i], vars, vars0, j)), length(ee), length(vars)),
4030
4031 finite_difference_partial (e, vars, vars0, j) :=
4032     block ([e0, e1, vars1],
4033            subst (vars0, e),
4034            e0: ev (%%, nouns, numer),
4035            vars1: copy (vars0),
4036            vars1[j]: (lhs(vars0[j]) = rhs(vars0[j]) + finite_difference_epsilon),
4037            subst (vars1, e),
4038            e1: ev (%%, nouns, numer),
4039            (e1 - e0)/finite_difference_epsilon),
4040
4041 finite_difference_gradient (e, vars, vars0) :=
4042     makelist (finite_difference_partial (e, vars, vars0, j), j, 1, length(vars)),
4043
4044 /* definition of hypergeometric_regularized can go away when it is built in ... */
4045 hypergeometric_regularized (p, q, x) :=
4046     hypergeometric (p, q, x)/product (gamma(q[i]), i, 1, length(q)),
4047
4048 test_derivative (expr, grad_expr, vars, vars0) :=
4049     block ([grad_via_diff, grad_via_finite_difference, norm_error],
4050            subst (vars0, grad_expr),
4051            grad_via_diff: ev (%%, numer),
4052            grad_via_finite_difference: finite_difference_gradient (expr, vars, vars0),
4053            norm_error: sqrt (lsum (e^2, e, grad_via_diff - grad_via_finite_difference)),
4054            if norm_error < sqrt (finite_difference_epsilon)
4055                then true
4056                else FAILED ('expr = expr, 'grad_via_diff = grad_via_diff, 'grad_via_finite_difference = grad_via_finite_difference, 'norm_error = norm_error)),
4057
4058 ab0_list:
4059     cartesian_product_list
4060         ([a = 2, a = 3, a = 4],
4061          [b = 4, b = 5, b = 6]),
4062
4063 abx0_list:
4064     cartesian_product_list
4065         ([a = 2, a = 3, a = 4],
4066          [b = 4, b = 5, b = 6],
4067          [x = 0.25, x = 0.5, x = 0.75]),
4068
4069 abxy0_list:
4070     cartesian_product_list
4071         ([a = 2, a = 3, a = 4],
4072          [b = 4, b = 5, b = 6],
4073          [x = 0.1, x = 0.25, x = 0.4],
4074          [y = 0.6, y = 0.75, y = 0.9]),
4075 0);
4076 0;
4077
4078 grad_beta: makelist (diff (beta (a, b), v), v, [a, b]);
4079 [-beta(a,b)*(psi[0](b+a)-psi[0](a)),
4080  -beta(a,b)*(psi[0](b+a)-psi[0](b))];
4081
4082 grad_beta_incomplete: makelist (diff (beta_incomplete (a, b, x), v), v, [a, b, x]);
4083 [beta_incomplete(a,b,x)*log(x)
4084   -gamma(a)^2*hypergeometric_regularized([a,a,1-b],[a+1,a+1],x)*x^a,
4085
4086  gamma(b)^2*hypergeometric_regularized([b,b,1-a],[b+1,b+1],1-x)*(1-x)^b
4087   -beta_incomplete(b,a,1-x)*log(1-x)+beta(a,b)*(psi[0](b)-psi[0](b+a)),
4088
4089  (1-x)^(b-1)*x^(a-1)];
4090
4091 grad_beta_incomplete_regularized: makelist (diff (beta_incomplete_regularized (a, b, x), v), v, [a, b, x]);
4092 [beta_incomplete_regularized(a,b,x)*(log(x)+psi[0](b+a)-psi[0](a))
4093   -(gamma(a)*hypergeometric_regularized([a,a,1-b],[a+1,a+1],x)*gamma(b+a)*x^a)/gamma(b),
4094
4095  (gamma(b)*gamma(b+a)*hypergeometric_regularized([b,b,1-a],[b+1,b+1],1-x)*(1-x)^b)/gamma(a)
4096    +beta_incomplete_regularized(b,a,1-x)*((-log(1-x))-psi[0](b+a)+psi[0](b)),
4097
4098  ((1-x)^(b-1)*x^(a-1))/beta(a,b)];
4099
4100 grad_beta_incomplete_generalized: makelist (diff (beta_incomplete_generalized (a, b, x, y), v), v, [a, b, x, y]);
4101 [beta_incomplete(a,b,y)*log(y)
4102   +gamma(a)^2*(hypergeometric_regularized([a,a,1-b],[a+1,a+1],x)*x^a
4103                 -hypergeometric_regularized([a,a,1-b],[a+1,a+1],y)*y^a)
4104   -beta_incomplete(a,b,x)*log(x),
4105
4106  (-gamma(b)^2*(hypergeometric_regularized([b,b,1-a],[b+1,b+1],1-x)*(1-x)^b
4107                -hypergeometric_regularized([b,b,1-a],[b+1,b+1],1-y)*(1-y)^b))
4108   -beta_incomplete(b,a,1-y)*log(1-y)+beta_incomplete(b,a,1-x)*log(1-x),
4109
4110  -(1-x)^(b-1)*x^(a-1),
4111
4112  (1-y)^(b-1)*y^(a-1)];
4113
4114 map (lambda ([ab0], test_derivative (beta (a, b), grad_beta, [a, b], ab0)), ab0_list);
4115 ''(makelist (true, length (ab0_list)));
4116
4117 map (lambda ([abx0], test_derivative (beta_incomplete (a, b, x), grad_beta_incomplete, [a, b, x], abx0)), abx0_list);
4118 ''(makelist (true, length (abx0_list)));
4119
4120 map (lambda ([abx0], test_derivative (beta_incomplete_regularized (a, b, x), grad_beta_incomplete_regularized, [a, b, x], abx0)), abx0_list);
4121 ''(makelist (true, length (abx0_list)));
4122
4123 map (lambda ([abxy0], test_derivative (beta_incomplete_generalized (a, b, x, y), grad_beta_incomplete_generalized, [a, b, x, y], abxy0)), abxy0_list);
4124 ''(makelist (true, length (abxy0_list)));
4125
4126 /******************************************************************************
4127
4128   Test the Beta function
4129
4130   Numerical values are taken from functions/wolfram.com.
4131 ******************************************************************************/
4132
4133 /* Derivatives of beta function */
4134 diff(beta(a,b),a);
4135 -(beta(a,b)*(psi[0](b+a)-psi[0](a)));
4136
4137 diff(beta(a,b),b);
4138 -(beta(a,b)*(psi[0](b+a)-psi[0](b)));
4139
4140 diff(beta(a,b),a,2);
4141 beta(a,b)*(psi[0](b+a)-psi[0](a))^2-beta(a,b)*(psi[1](b+a)-psi[1](a));
4142
4143 diff(beta(a,b),b,2);
4144 beta(a,b)*(psi[0](b+a)-psi[0](b))^2-beta(a,b)*(psi[1](b+a)-psi[1](b));
4145
4146 diff(beta(a,b),a,1,b,1);
4147 beta(a,b)*(psi[0](b+a)-psi[0](a))*(psi[0](b+a)-psi[0](b))
4148  -beta(a,b)*psi[1](b+a);