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 (1-(-a-1)/z)*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   2.81e-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  2.81e-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 /******************************************************************************
2793    Test the Fresnel Integrals S(z) and C(z)
2794 ******************************************************************************/
2795
2796 /* Specific values for the Fresnel Integrals */
2797
2798 map(fresnel_s,[0,0.0,0.0b0]);
2799 [0,0.0,0.b0];
2800
2801 map(fresnel_c,[0,0.0,0.0b0]);
2802 [0,0.0,0.b0];
2803
2804 limit(fresnel_c(x),x,inf);
2805 1/2;
2806
2807 limit(fresnel_s(x),x,inf);
2808 1/2;
2809
2810 limit(fresnel_c(x),x,minf);
2811 -1/2;
2812
2813 limit(fresnel_s(x),x,minf);
2814 -1/2;
2815
2816 /* Simplification of infinities
2817    The rules for an odd function and the simplifaction for imaginary
2818    arguments are applied too.
2819 */
2820
2821 map(fresnel_s,[inf,-inf,minf,-minf,%i*inf,-%i*inf,%i*minf,-%i*minf]);
2822 [1/2,-1/2,-1/2,1/2,-%i/2,%i/2,%i/2,-%i/2];
2823
2824 map(fresnel_c,[inf,-inf,minf,-minf,%i*inf,-%i*inf,%i*minf,-%i*minf]);
2825 [1/2,-1/2,-1/2,1/2,%i/2,-%i/2,-%i/2,%i/2];
2826
2827 /* No simplification for other infinities and undeterminates */
2828
2829 map(fresnel_s,[infinity,und,ind]);
2830 [fresnel_s(infinity),fresnel_s(und),fresnel_s(ind)];
2831
2832 map(fresnel_c,[infinity,und,ind]);
2833 [fresnel_c(infinity),fresnel_c(und),fresnel_c(ind)];
2834
2835 /* The Fresnel Integrals S(z) and C(z) are odd functions
2836    A reflection rule is given and the rule odd-function-reflect is applied */
2837
2838 map(fresnel_s,[-x, (x-1), (-x+1), (-x-1)]);
2839 [-fresnel_s(x), fresnel_s(x-1), -fresnel_s(x-1),-fresnel_s(1+x)];
2840
2841 map(fresnel_c,[-x, (x-1), (-x+1), (-x-1)]);
2842 [-fresnel_c(x), fresnel_c(x-1), -fresnel_c(x-1),-fresnel_c(1+x)];
2843
2844 /* The Fresnel Integals simplify imaginary arguments */
2845
2846 map(fresnel_s, [%i,%i*x,-%i*x,%i*(x+1)]);
2847 [-%i*fresnel_s(1),-%i*fresnel_s(x),%i*fresnel_s(x),-%i*fresnel_s(x+1)];
2848
2849 map(fresnel_c, [%i,%i*x,-%i*x,%i*(x+1)]);
2850 [%i*fresnel_c(1),%i*fresnel_c(x),-%i*fresnel_c(x),%i*fresnel_c(x+1)];
2851
2852 /* The Fresnel Integrals have Mirror Symmetry */
2853
2854 declare(z,complex);
2855 done;
2856
2857 conjugate(fresnel_s(z));
2858 fresnel_s(conjugate(z));
2859
2860 conjugate(fresnel_s(x+%i*y));
2861 fresnel_s(x-%i*y);
2862
2863 conjugate(fresnel_c(z));
2864 fresnel_c(conjugate(z));
2865
2866 conjugate(fresnel_c(x+%i*y));
2867 fresnel_c(x-%i*y);
2868
2869 /* Taylor expansion of the Fresnel Integrals to order O(z^12) */
2870
2871 /* Expand the function */
2872 taylor(fresnel_s(x),x,0,12);
2873 %pi*x^3/6 - %pi^3*x^7/336 + %pi^5 *x^11/42240;
2874
2875 taylor(fresnel_c(x),x,0,12);
2876 x - %pi^2*x^5/40 + %pi^4 *x^9/3456;
2877
2878 /* Expand the argument and apply the function */
2879 fresnel_s(taylor(x,x,0,12));
2880 %pi*x^3/6 - %pi^3*x^7/336 + %pi^5 *x^11/42240;
2881
2882 fresnel_c(taylor(x,x,0,12));
2883 x - %pi^2*x^5/40 + %pi^4 *x^9/3456;
2884
2885 /* Differentiation of the Fresnel Integrals */
2886
2887 diff(fresnel_s(x),x);
2888 sin(%pi*x^2/2);
2889
2890 diff(fresnel_c(x),x);
2891 cos(%pi*x^2/2);
2892
2893 /* The elementary Integral of the Fresnel Integrals
2894    More complicated integrals can be found in rtest_integrate_special.mac */
2895
2896 integrate(fresnel_s(x),x);
2897 x*fresnel_s(x)+1/%pi*cos(%pi*x^2/2);
2898
2899 integrate(fresnel_c(x),x);
2900 x*fresnel_c(x)-1/%pi*sin(%pi*x^2/2);
2901
2902 /* Representation of the Fresnel Integrals through the Error function Erf */
2903
2904 erf_representation:true;
2905 true;
2906
2907 fresnel_s(x);
2908 (1+%i)/4*(erf((1+%i)/2*sqrt(%pi)*x) - %i* erf((1-%i)/2*sqrt(%pi)*x));
2909
2910 fresnel_c(x);
2911 (1-%i)/4*(erf((1+%i)/2*sqrt(%pi)*x) + %i* erf((1-%i)/2*sqrt(%pi)*x));
2912
2913 erf_representation:false;
2914 false;
2915
2916 /* Representation of the Fresnel Integrals through the
2917    Hypergeometric function */
2918
2919 hypergeometric_representation:true;
2920 true;
2921
2922 fresnel_s(x);
2923 %pi*x^3/6*hypergeometric([3/4],[3/2,7/4],-%pi^2*x^4/16);
2924
2925 fresnel_c(x);
2926 x*hypergeometric([1/4],[1/2,5/4],-%pi^2*x^4/16);
2927
2928 hypergeometric_representation:false;
2929 false;
2930
2931 /* Numerical tests for the Fresnel Integrals */
2932
2933 fpprec:64;
2934 64;
2935
2936 /* Tests for fresnel_s
2937    Real argument 0.25,0.50,1.0,2.0,5.0,10.0 */
2938
2939 relerror(
2940   fresnel_s(0.25),
2941   0.008175600235777755778102308866942774752486734698017086013976457144b0,
2942   1e-15);
2943 true;
2944
2945 relerror(
2946   fresnel_s(0.25b0),
2947   0.008175600235777755778102308866942774752486734698017086013976457144b0,
2948   1b-64);
2949 true;
2950
2951 relerror(
2952   fresnel_s(0.50),
2953   0.06473243285999927761148051223061476765072591849351249278758894565b0,
2954   1e-15);
2955 true;
2956
2957 relerror(
2958   fresnel_s(0.50b0),
2959   0.06473243285999927761148051223061476765072591849351249278758894565b0,
2960   1b-64);
2961 true;
2962
2963 relerror(
2964   fresnel_s(1.0),
2965   0.4382591473903547660767566966251526374937865724524165673344073263b0,
2966   1e-15);
2967 true;
2968
2969 relerror(
2970   fresnel_s(1.0b0),
2971   0.4382591473903547660767566966251526374937865724524165673344073263b0,
2972   2.0b-64);
2973 true;
2974
2975 relerror(
2976   fresnel_s(2.0),
2977   0.3434156783636982421953008159580684568865418122025247675792689204b0,
2978   1e-15);
2979 true;
2980
2981 relerror(
2982   fresnel_s(2.0b0),
2983   0.3434156783636982421953008159580684568865418122025247675792689204b0,
2984   1.5b-63);
2985 true;
2986
2987 relerror(
2988   fresnel_s(5.0),
2989   0.4991913819171168867519283804659916554084319970723881534101411152b0,
2990   1e-15);
2991 true;
2992
2993 relerror(
2994   fresnel_s(5.0b0),
2995   0.4991913819171168867519283804659916554084319970723881534101411152b0,
2996   1b-64);
2997 true;
2998
2999 relerror(
3000   fresnel_s(10.0),
3001   0.4681699785848822404033511108104469460538427245558302799270062272b0,
3002   1e-15);
3003 true;
3004
3005 relerror(
3006   fresnel_s(10.0b0),
3007   0.4681699785848822404033511108104469460538427245558302799270062272b0,
3008   1b-64);
3009 true;
3010
3011 /* Complex argument 0.25+%i,0.50+%i,1.0+%i,2.0+%i,5.0+%i */
3012
3013 relerror(
3014   fresnel_s(0.25+%i),
3015   -0.2762104914409824591766528447060750469693825567583676638192814447b0
3016   - 0.4331061708987372646968782339728450577465559834165609320175098398b0*%i,
3017   1e-15);
3018 true;
3019
3020 relerror(
3021   fresnel_s(0.25b0+%i),
3022   -0.2762104914409824591766528447060750469693825567583676638192814447b0
3023   - 0.4331061708987372646968782339728450577465559834165609320175098398b0*%i,
3024   1b-64);
3025 true;
3026
3027 relerror(
3028   fresnel_s(0.50+%i),
3029    -0.7169788451833594258616872412575572495423663980107873580716959051b0
3030   - 0.3393082523853171783825858066689424560943675693298848574923503884b0*%i,
3031   1e-15);
3032 true;
3033
3034 relerror(
3035   fresnel_s(0.50b0+%i),
3036    -0.7169788451833594258616872412575572495423663980107873580716959051b0
3037   - 0.3393082523853171783825858066689424560943675693298848574923503884b0*%i,
3038   1b-64);
3039 true;
3040
3041 relerror(
3042   fresnel_s(1.0+%i),
3043   -2.061888219194840468080716536685708600815908323737868052048638806b0
3044   + 2.061888219194840468080716536685708600815908323737868052048638806b0*%i,
3045   1e-15);
3046 true;
3047
3048 relerror(
3049   fresnel_s(1.0b0+%i),
3050   -2.061888219194840468080716536685708600815908323737868052048638806b0
3051   + 2.061888219194840468080716536685708600815908323737868052048638806b0*%i,
3052   2.0b-64);
3053 true;
3054
3055 relerror(
3056   fresnel_s(2.0+%i),
3057   -15.58775110440458732748278797797881643198730378904101846898562610b0
3058   - 36.72546488399143842838787627677917885752587065976755449373500438b0*%i,
3059   4.0e-13);
3060 true;
3061
3062 relerror(
3063   fresnel_s(2.0b0+%i),
3064   -15.58775110440458732748278797797881643198730378904101846898562610b0
3065   - 36.72546488399143842838787627677917885752587065976755449373500438b0*%i,
3066   4.0b-62);
3067 true;
3068
3069 relerror(
3070   fresnel_s(5.0+%i),
3071   -204452.5505063210590493330126425293361797143144299005524393297869b0
3072   + 38438.9093777759513198736581956757227475265101347572827364567521b0*%i,
3073   3.0e-14);
3074 true;
3075
3076 relerror(
3077   fresnel_s(5.0b0+%i),
3078   -204452.5505063210590493330126425293361797143144299005524393297869b0
3079   + 38438.9093777759513198736581956757227475265101347572827364567521b0*%i,
3080   4.0b-63);
3081 true;
3082
3083 /* Tests for fresnel_c
3084    Real argument 0.25,0.50,1.0,2.0,5.0,10.0 */
3085
3086 relerror(
3087   fresnel_c(0.25),
3088   0.2497591503565431834592215178008857243781399770549380697377810451b0,
3089   1e-15);
3090 true;
3091
3092 relerror(
3093   fresnel_c(0.25b0),
3094   0.2497591503565431834592215178008857243781399770549380697377810451b0,
3095   2.0b-64);
3096 true;
3097
3098 relerror(
3099   fresnel_c(0.50),
3100   0.4923442258714463928788436651566816377660951457715012532946526193b0,
3101   1e-15);
3102 true;
3103
3104 relerror(
3105   fresnel_c(0.50b0),
3106   0.4923442258714463928788436651566816377660951457715012532946526193b0,
3107   1b-64);
3108 true;
3109
3110 relerror(
3111   fresnel_c(1.0),
3112   0.7798934003768228294742064136526901366306257081363209601031335832b0,
3113   1e-15);
3114 true;
3115
3116 relerror(
3117   fresnel_c(1.0b0),
3118   0.7798934003768228294742064136526901366306257081363209601031335832b0,
3119   2.0b-64);
3120 true;
3121
3122 relerror(
3123   fresnel_c(2.0),
3124   0.4882534060753407545002235033572610376883671545092153829475964427b0,
3125   1e-15);
3126 true;
3127
3128 relerror(
3129   fresnel_c(2.0b0),
3130   0.4882534060753407545002235033572610376883671545092153829475964427b0,
3131   5.1b-64);
3132 true;
3133
3134 relerror(
3135   fresnel_c(5.0),
3136   0.5636311887040122311021074044130139641207537623099921078616593412b0,
3137   1e-15);
3138 true;
3139
3140 relerror(
3141   fresnel_c(5.0b0),
3142   0.5636311887040122311021074044130139641207537623099921078616593412b0,
3143   2.0b-64);
3144 true;
3145
3146 relerror(
3147   fresnel_c(10.0),
3148   0.4998986942055157236141518477356211143923468402262626572074674093b0,
3149   1.73e-15);
3150 true;
3151
3152 relerror(
3153   fresnel_c(10.0b0),
3154   0.4998986942055157236141518477356211143923468402262626572074674093b0,
3155   2.47b-64);
3156 true;
3157
3158 /* Complex argument 0.25+%i,0.50+%i,1.0+%i,2.0+%i,5.0+%i */
3159
3160 relerror(
3161   fresnel_c(0.25+%i),
3162   0.0097446563393801078320153856258158947458946246448139394089371651b0
3163   + 0.8830495953515316267486626148763161568578011755347662219418704558b0*%i,
3164   1e-15);
3165 true;
3166
3167 relerror(
3168   fresnel_c(0.25b0+%i),
3169   0.0097446563393801078320153856258158947458946246448139394089371651b0
3170   + 0.8830495953515316267486626148763161568578011755347662219418704558b0*%i,
3171   1b-64);
3172 true;
3173
3174 relerror(
3175   fresnel_c(0.50+%i),
3176   0.1199549363708813724035204126626172258713764410185526201803481040b0
3177   + 1.2468579809337107891237368150539206578767217615436096318514136378b0*%i,
3178   1e-15);
3179 true;
3180
3181 relerror(
3182   fresnel_c(0.50b0+%i),
3183   0.1199549363708813724035204126626172258713764410185526201803481040b0
3184   + 1.2468579809337107891237368150539206578767217615436096318514136378b0*%i,
3185   1b-64);
3186 true;
3187
3188 relerror(
3189   fresnel_c(1.0+%i),
3190   2.555793778102439024634522388352195842156623604203584296357752992b0
3191   + 2.555793778102439024634522388352195842156623604203584296357752992b0*%i,
3192   1e-15);
3193 true;
3194
3195 relerror(
3196   fresnel_c(1.0b0+%i),
3197   2.555793778102439024634522388352195842156623604203584296357752992b0
3198   + 2.555793778102439024634522388352195842156623604203584296357752992b0*%i,
3199   2.0b-64);
3200 true;
3201
3202 relerror(
3203   fresnel_c(2.0+%i),
3204   -36.22568799288165021578757830205090012554103292231420092205271528b0
3205   + 16.08787137412548041729488986757594802326683407694887086741793640b0*%i,
3206   4.0e-13);
3207 true;
3208
3209 relerror(
3210   fresnel_c(2.0b0+%i),
3211   -36.22568799288165021578757830205090012554103292231420092205271528b0
3212   + 16.08787137412548041729488986757594802326683407694887086741793640b0*%i,
3213   4.0b-62);
3214 true;
3215
3216 relerror(
3217   fresnel_c(5.0+%i),
3218   38439.4093777740143202918713550472184235160647415045418329908291b0
3219   + 204453.0505063119004499308897858846121088557663560705550579996093b0*%i,
3220   3.0e-14);
3221 true;
3222
3223 relerror(
3224   fresnel_c(5.0b0+%i),
3225   38439.4093777740143202918713550472184235160647415045418329908291b0
3226   + 204453.0505063119004499308897858846121088557663560705550579996093b0*%i,
3227   4.0b-63);
3228 true;
3229
3230 /* Bug #2509 fresnel_s incorrect for small values */
3231 relerror(
3232   fresnel_c(1d-20),
3233   1d-20,
3234   1d-15);
3235 true;
3236
3237 relerror(
3238   fresnel_c(1b-40),
3239   1b-40,
3240   5b-63);
3241 true;
3242
3243 relerror(
3244   fresnel_s(1d-20),
3245   float(%pi/6*(1d-20)^3),
3246   1d-15);
3247 true;
3248
3249 relerror(
3250   fresnel_s(1b-40),
3251   bfloat(%pi/6*(1b-40)^3),
3252   5b-63);
3253 true;
3254
3255 /******************************************************************************
3256    Test the Beta incomplete function
3257 ******************************************************************************/
3258
3259 /* Specialized values */
3260
3261 (assume(am<0,ap>0,b>0),done);
3262 done;
3263
3264 beta_incomplete(a,b,1);
3265 beta(a,b);
3266
3267 beta_incomplete(ap,b,0);
3268 0;
3269
3270 errcatch(beta_incomplete(am,b,0));
3271 [];
3272
3273 (forget(am<0, ap>0,b>0),done);
3274 done;
3275
3276 /* b is a positive integer */
3277
3278 beta_incomplete(a,1,z);
3279 z^a/a;
3280 beta_incomplete(a,2,z);
3281 (a*(1-z)+1)*z^a/(a*(a+1));
3282 beta_incomplete(1,1,z);
3283 z;
3284
3285 /* b is a positive integer.
3286    Check the handling of float and bigfloat numbers */
3287 fpprec:16;
3288 16;
3289
3290 beta_incomplete(1.0,1,z);
3291 1.0*z;
3292 /* Unfortunate, but we don't get exactly 1.0b0 */
3293 (closeto(beta_incomplete(1.0b0,1,z)/z,1b0,1b-15));
3294 true;
3295 beta_incomplete(1.0,1,1/2);
3296 0.5;
3297 beta_incomplete(1.0b0,1,1/2);
3298 0.5b0;
3299 beta_incomplete(1,1,1/2+%i);
3300 %i+1/2;
3301 beta_incomplete(1.0,1,1/2+%i);
3302 1.0*%i+0.5;
3303 beta_incomplete(1.0b0,1,1/2+%i);
3304 1.0b0*%i+0.5b0;
3305 beta_incomplete(1,2,1/2);
3306 3/8;
3307 beta_incomplete(1.0,2,1/2);
3308 0.375;
3309 beta_incomplete(1.0b0,2,1/2);
3310 0.375b0;
3311 beta_incomplete(1,2,1/2+%i);
3312 (3/2-%i)*(1/2+%i)/2;
3313 beta_incomplete(1.0,2,1/2+%i);
3314 0.875+0.5*%i;
3315 beta_incomplete(1.0b0,2,1/2+%i);
3316 0.875b0+0.5b0*%i;
3317
3318 /* We check the expansion for b positive against the numerical evaluation */
3319 closeto(
3320   float(beta_incomplete(1,1,1/2)) - beta_incomplete(1.0,1.0,0.5),
3321   0.0,
3322   1e-15);
3323 true;
3324 closeto(
3325   float(beta_incomplete(3/2,2,1/2)) - beta_incomplete(1.5,2.0,0.5),
3326   0.0,
3327   1e-15);
3328 true;
3329 closeto(
3330   float(beta_incomplete(2,3,1/2)) - beta_incomplete(2.0,3.0,0.5),
3331   0.0,
3332   1e-15);
3333 true;
3334 closeto(
3335   float(beta_incomplete(5/2,4,1/2)) - beta_incomplete(2.5,4.0,0.5),
3336   0.0,
3337   1e-15);
3338 true;
3339 closeto(
3340   float(beta_incomplete(3,5,1/2)) - beta_incomplete(3.0,5.0,0.5),
3341   0.0,
3342   1e-15);
3343 true;
3344 closeto(
3345   float(beta_incomplete(2,6,1/2+%i)) - beta_incomplete(2.0,6.0,0.5+%i),
3346   0.0,
3347   1e-15);
3348 true;
3349 closeto(
3350   /* We do this extra rectform, because of a bug in abs for expressions
3351      with a complex exponent */
3352   rectform(float(beta_incomplete(2+%i,7,1/2+%i))
3353                - beta_incomplete(2.0+%i,7.0,0.5+%i)),
3354   0.0,
3355   1e-15);
3356 true;
3357
3358 /* a is a positive integer we expand */
3359 beta_incomplete(1,b,z);
3360 (1-(1-z)^b)/b;
3361 beta_incomplete(2,b,z);
3362  (1-(1-z)^b*(b*z+1))/(b*(b+1));
3363 beta_incomplete(3,b,z);
3364  2*(1-(1-z)^b*(b*(b+1)*z^2/2+b*z+1))/(b*(b+1)*(b+2));
3365
3366 /* a is a positive integer.
3367    Check the expansion against numerically evaluation.
3368    First we test for b a positive value. */
3369
3370 closeto(
3371   float(beta_incomplete(1,2,3/2))-beta_incomplete(1.0,2.0,1.5),
3372   0.0,
3373   1e-15);
3374 true;
3375 closeto(
3376   float(beta_incomplete(2,5/2,3/2))-beta_incomplete(2.0,2.5,1.5),
3377   0.0,
3378   1e-15);
3379 true;
3380 closeto(
3381   float(beta_incomplete(3,3,3/2))-beta_incomplete(3.0,3.0,1.5),
3382   0.0,
3383   1e-15);
3384 true;
3385 closeto(
3386   float(beta_incomplete(2,7/2,3/2))-beta_incomplete(2.0,3.5,1.5),
3387   0.0,
3388   1e-15);
3389 true;
3390 closeto(
3391   float(beta_incomplete(2,5/2,3/2+%i))-beta_incomplete(2.0,2.5,1.5+%i),
3392   0.0,
3393   1e-15);
3394 true;
3395 /* I (rtoy) think the limiting accuracy of the following test is the conversion of
3396  * the exact answer to floating point.  If the exact answer is converted to a
3397  * bfloat 32 digits, the difference is less than 1.6e-15 or so.
3398 */
3399 closeto(
3400   float(rectform(float(beta_incomplete(2,5/2,-3/2+%i))))
3401                 -beta_incomplete(2.0,2.5,-1.5+%i),
3402   0.0,
3403   4.35e-15); /* we have lost accuracy. More tests necessary? */
3404 true;
3405
3406 /* Incomplete Beta is definied for negative integers a and b >= (-a)
3407    Add this point we have a problem. functions.wolfram.com gives different
3408    numerical results for this cases. When b not equal -a Maxima and
3409    functions.wolfram.com differ by a factor 2. For other valid integer b
3410    functions.wolfram.com gives ComplexInfinity and not an expected result.
3411    What is wrong?
3412 */
3413
3414 beta_incomplete(-1,1,z);
3415 -1/z;
3416 beta_incomplete(-2,1,z);
3417 -1/(2*z^2);
3418 beta_incomplete(-2,2,z);
3419 (z-1/2)/z^2;
3420 beta_incomplete(-3,1,z);
3421 -1/(3*z^3);
3422 beta_incomplete(-3,2,z);
3423 (z/2-1/3)/z^3;
3424 beta_incomplete(-3,3,z);
3425 (-z^2+z-1/3)/z^3;
3426
3427 /* Some numerical tests in double float precision */
3428
3429 closeto(
3430   beta_incomplete(0.5,1.0,0.10),
3431   0.6324555320336758663997787088865437067439110278650433653715009706b0,
3432   1e-15);
3433 true;
3434
3435 closeto(
3436   beta_incomplete(0.5,1.0,0.15),
3437   0.7745966692414833770358530799564799221665843410583181653175147532b0,
3438   1e-15);
3439 true;
3440
3441 closeto(
3442   beta_incomplete(0.5,1.0,0.20),
3443   0.8944271909999158785636694674925104941762473438446102897083588982b0,
3444   1e-15);
3445 true;
3446
3447 closeto(
3448   beta_incomplete(0.5,1.0,0.25),
3449   1.000000000000000000000000000000000000000000000000000000000000000b0,
3450   1e-15);
3451 true;
3452
3453 closeto(
3454   beta_incomplete(0.5,1.0,0.50),
3455   1.414213562373095048801688724209698078569671875376948073176679738b0,
3456   1.416e-15);
3457 true;
3458
3459 closeto(
3460   beta_incomplete(0.5,1.0,0.75),
3461   1.732050807568877293527446341505872366942805253810380628055806979b0,
3462   2.11e-15);
3463 true;
3464
3465 closeto(
3466   beta_incomplete(0.5,1.0,1.00),
3467   2.000000000000000000000000000000000000000000000000000000000000000b0,
3468   2.0e-15);
3469 true;
3470
3471 closeto(
3472   beta_incomplete(0.5,1.0,1.25),
3473   2.236067977499789696409173668731276235440618359611525724270897245b0,
3474   2.0e-15);
3475 true;
3476
3477 closeto(
3478   beta_incomplete(0.5,1.0,1.50),
3479   2.449489742783178098197284074705891391965947480656670128432692567b0,
3480   2.0e-15);
3481 true;
3482
3483 closeto(
3484   beta_incomplete(0.5,1.0,1.75),
3485   2.645751311064590590501615753639260425710259183082450180368334459b0,
3486   1.446e-15);
3487 true;
3488
3489 closeto(
3490   beta_incomplete(0.5,1.0,2.00),
3491   2.828427124746190097603377448419396157139343750753896146353359476b0,
3492   1.503e-15);
3493 true;
3494
3495 /* Some numerical tests in bigfloat precision */
3496
3497 fpprec:64;
3498 64;
3499
3500 closeto(
3501   beta_incomplete(0.5,1.0,0.10b0),
3502   0.6324555320336758663997787088865437067439110278650433653715009706b0,
3503   1b-60);
3504 true;
3505
3506 closeto(
3507   beta_incomplete(0.5,1.0,0.15b0),
3508   0.7745966692414833770358530799564799221665843410583181653175147532b0,
3509   1b-60);
3510 true;
3511
3512 closeto(
3513   beta_incomplete(0.5,1.0,0.20b0),
3514   0.8944271909999158785636694674925104941762473438446102897083588982b0,
3515   1b-60);
3516 true;
3517
3518 closeto(
3519   beta_incomplete(0.5,1.0,0.25b0),
3520   1.000000000000000000000000000000000000000000000000000000000000000b0,
3521   1b-60);
3522 true;
3523
3524 closeto(
3525   beta_incomplete(0.5,1.0,0.50b0),
3526   1.414213562373095048801688724209698078569671875376948073176679738b0,
3527   1b-60);
3528 true;
3529
3530 closeto(
3531   beta_incomplete(0.5,1.0,0.75b0),
3532   1.732050807568877293527446341505872366942805253810380628055806979b0,
3533   1b-16);
3534 true;
3535
3536 closeto(
3537   beta_incomplete(0.5,1.0,1.00b0),
3538   2.000000000000000000000000000000000000000000000000000000000000000b0,
3539   1b-60);
3540 true;
3541
3542 closeto(
3543   beta_incomplete(0.5,1.0,1.25b0),
3544   2.236067977499789696409173668731276235440618359611525724270897245b0,
3545   1b-60);
3546 true;
3547
3548 closeto(
3549   beta_incomplete(0.5,1.0,1.50b0),
3550   2.449489742783178098197284074705891391965947480656670128432692567b0,
3551   1b-60);
3552 true;
3553
3554 closeto(
3555   beta_incomplete(0.5,1.0,1.75b0),
3556   2.645751311064590590501615753639260425710259183082450180368334459b0,
3557   1b-60);
3558 true;
3559
3560 closeto(
3561   beta_incomplete(0.5,1.0,2.00b0),
3562   2.828427124746190097603377448419396157139343750753896146353359476b0,
3563   1b-60);
3564 true;
3565
3566 /* See Bug 3220128, but this isn't really that bug */
3567 closeto(
3568   gamma_incomplete(0, 200b0*%i),
3569   0.00437844609302782567916569771749325524128345091344187598851110680706344144459295b0
3570   - %i*.00241398745542678587253611621620491057595401709907514761094360488114169654741b0,
3571   3.4b-67);
3572 true;
3573
3574 /* Make sure we don't overflow unnecessarily in gamma_incomplete_regularized */
3575 closeto(
3576   gamma_incomplete_regularized(45001d0, 45000d0),
3577   0.5012537510700691773177801688515861486945632498553288,
3578   6.4e-12);
3579 true;
3580
3581 /* Check accuracy on the negative real axis.
3582  * See Bug  3526359 - gamma_incomplete(1/5,-32.0) not accurate
3583  */
3584 relerror(
3585   gamma_incomplete(1/5,-32.0),
3586   -4.0986398326716649399d12 - %i*2.9778361454160762231d12,
3587   2.545e-15);
3588 true;
3589
3590 relerror(
3591   gamma_incomplete(10,-100d0),
3592   subst(x=-100d0, block([gamma_expand : true], gamma_incomplete(10,x))),
3593   8.78e-15);
3594 true;
3595
3596 relerror(
3597   gamma_incomplete(10,-100b0),
3598   subst(x=-100b0, block([gamma_expand : true], gamma_incomplete(10,x))),
3599   3.57b-64);
3600 true;
3601
3602 relerror(
3603   gamma_incomplete(10,-100d0+%i),
3604   subst(x=-100d0+%i, block([gamma_expand : true], gamma_incomplete(10,x))),
3605   1.046e-14);
3606 true;
3607
3608 relerror(
3609   gamma_incomplete(10,-100b0+%i),
3610   subst(x=-100b0+%i, block([gamma_expand : true], gamma_incomplete(10,x))),
3611   2.04b-64);
3612 true;
3613
3614 relerror(
3615  erf(inverse_erf(.5)),
3616  0.5,
3617  1.4d-14);
3618 true;
3619
3620 relerror(
3621  erf(inverse_erf(-.5)),
3622  -0.5,
3623  1.4e-14);
3624 true;
3625
3626 relerror(
3627  erf(inverse_erf(.5b0)),
3628  0.5b0,
3629  2b-64);
3630 true;
3631
3632 relerror(
3633  erf(inverse_erf(-.5b0)),
3634  -0.5b0,
3635  2b-64);
3636 true;
3637
3638 relerror(
3639  erf(inverse_erf(2.0)),
3640  2.0,
3641  7.9d-15);
3642 true;
3643
3644 relerror(
3645  erf(inverse_erf(-2.0)),
3646  -2.0,
3647  7.9d-15);
3648 true;
3649
3650 /* These are a bit slow if fpprec is large.  Make fpprec a bit smaller */
3651 fpprec:24;
3652 24;
3653
3654 relerror(
3655  erf(inverse_erf(2b0)),
3656  2b0,
3657  5.9b-24);
3658 true;
3659
3660 relerror(
3661  erf(inverse_erf(-2b0)),
3662  -2b0,
3663  5.9b-24);
3664 true;
3665
3666 relerror(
3667  erf(inverse_erf(2b0+2b0*%i)),
3668  2b0+2b0*%i,
3669  6.6b-25);
3670 true;
3671
3672 relerror(
3673  erf(inverse_erf(-2b0-2b0*%i)),
3674  -2b0-2b0*%i,
3675  6.6b-25);
3676 true;
3677
3678 relerror(
3679  erf(inverse_erf(10b0+1000b0*%i)),
3680  10b0+1000b0*%i,
3681  2b-23);
3682 true;
3683
3684
3685 /* From the mailing list.  This should not signal an error and should give a simplifed answer */
3686 relerror(
3687  expand(bfloat(erf((sqrt(2)*%i-sqrt(2))/2))),
3688  4.74147636640994245161681b-1 * %i - 9.69264211944215930381491b-1,
3689  1b-23);
3690 true;
3691
3692 /* Bug #2619
3693  * Function inverse_erf - error in numerical evaluation
3694  */
3695 relerror(inverse_erf(0.7715),
3696  .8515204685911937,
3697  1d-15);
3698 true;
3699
3700 relerror(inverse_erf(- .9763545580611441 ),
3701  -1.600070795398459,
3702  1d-15);
3703 true;
3704
3705 /* Bug #2668 Bigfloat Gamma inaccurae for small inputs
3706  *
3707  * For these small values, gamma(z) = gamma(z+1)/z = 1/z
3708  * since 1+z = 1 and gamma(1) = 1.
3709  */
3710 relerror(gamma(2.0^-256), 2.0^256, 1d-14);
3711 true;
3712
3713 relerror(gamma(2b0^-256), 2b0^256, 10^(-fpprec+1));
3714 true;
3715
3716 relerror(
3717   gamma(2.0^-256 + 1d-200*%i),
3718   rectform(1/(2.0^-256 + 1d-200*%i)),
3719   2d-15);
3720 true;
3721
3722 closeto(
3723   gamma(2b0^-256 + 1b-500*%i),
3724   rectform(1/(2b0^-256 + 1b-500*%i)),
3725   2d-15);
3726 true;
3727
3728 (fpprec:oldfpprec,done);
3729 done;
3730
3731 /* SF bug #3090: "erfi switches sign at approximately -0.476 and 0.476" */
3732
3733 block
3734   (kill (x, L1, L2),
3735    x : makelist (i/10.0, i, -20, 20),
3736    L1 : makelist (erfi (x1), x1, x),
3737    L2 : makelist (-%i*erf (%i*x1), x1, x),
3738    every (lambda ([x1, x2], x1 = x2), L1, L2));
3739 true;
3740
3741 every (lambda ([x1, e1], sign(x1) = sign(e1)), x, L1);
3742 true;
3743
3744 every (lambda ([x1, e1], sign(x1) = sign(e1)), x, L2);
3745 true;
3746
3747 /* SF bug #3277: "no numerical evaluation for gamma_greek" */
3748
3749 expand (gamma_incomplete_lower(10, 1));
3750 362880-986410*%e^-1;
3751
3752 /* numerical value not very accurate; perhaps it is being computed naively via preceding formula? */
3753 closeto (ev (gamma_incomplete_lower(10, 1), numer), 4.043407757955496e-2, 1e-10);
3754 true;
3755
3756 /* this one is way off the mark ... sheesh. */
3757 closeto (gamma_incomplete_lower(10, 1.0), 4.043407757955496e-2, 1e-8);
3758 true;
3759
3760 (reset (fpprec),
3761  closeto (gamma_incomplete_lower(10, 1b0), 4.043407757955496b-2, 1b-10));
3762 true;
3763
3764 gamma_incomplete_lower(1, 10);
3765 1 - %e^-10;
3766
3767 closeto (ev (gamma_incomplete_lower(1, 10), numer), 0.9999546000702375e0, 1e-15);
3768 true;
3769
3770 closeto (gamma_incomplete_lower(1.0, 10), 0.9999546000702375e0, 1e-15);
3771 true;
3772
3773 closeto (gamma_incomplete_lower(1b0, 10), 0.9999546000702375b0, 1b-15);
3774 true;
3775
3776 /* complex float and complex bfloat */
3777
3778 /* yikes */
3779 closeto (ev (gamma_incomplete_lower(10, 1 + %i), numer), 7.992250288513467e-1*%i + 1.011121671734864e0, 1e-7);
3780 true;
3781
3782 /* ditto */
3783 closeto (gamma_incomplete_lower(10, 1.0 + %i), 7.992250288513467e-1*%i + 1.011121671734864e0, 1e-7);
3784 true;
3785
3786 closeto (gamma_incomplete_lower(10, 1b0 + %i), 7.992250288513467b-1*%i + 1.011121671734864b0, 1b-10);
3787 true;
3788
3789 /* mailing list 2018-02-04: "Bug in beta_incomplete_regularized." */
3790
3791 (kill(f, g, x, y),
3792  f(x) := beta_incomplete_regularized(2,2,x/(1+x)),
3793  g(x) := beta_incomplete(2,2,x),
3794  should_be_equal(a, b) := if equal(a, b) then true else notequal(a, b),
3795  0);
3796 0;
3797
3798 f(1/y);
3799 (2*(1-1/((1/y+1)*y))+1)/((1/y+1)^2*y^2);
3800
3801 f(1/x);
3802 (2*(1-1/((1/x+1)*x))+1)/((1/x+1)^2*x^2);
3803
3804 should_be_equal (f(1/x), subst (y = 1/x, f(y)));
3805 true;
3806
3807 g(-y);
3808 (y^2*(2*(y+1)+1))/6;
3809
3810 g(-x);
3811 (x^2*(2*(x+1)+1))/6;
3812
3813 should_be_equal (g(-x), subst (y = -x, g(y)));
3814 true;
3815
3816 (kill(a, h, i, j, k, u, v, w),
3817  h(u) := gamma_incomplete(3, u),
3818  i(v, w) := gamma_incomplete_generalized(3 + a, v, w),
3819  j(u) := gamma_incomplete_regularized(4, u),
3820  k(v, w) := beta_incomplete_regularized(2, v, w),
3821  0);
3822 0;
3823
3824 gamma_expand : true;
3825 true;
3826
3827 should_be_equal (h(u), integrate (t^2*exp(-t),t,u,inf));
3828 true;
3829
3830 [h(u), h(-u)];
3831 [2*(u^2/2+u+1)*%e^-u,2*(u^2/2-u+1)*%e^u];
3832
3833 should_be_equal (h(-u), subst (v = -u, h(v)));
3834 true;
3835
3836 should_be_equal (i(2*v, 4*w), a*(a+1)*(a+2)
3837                                *((-((4^(a+2)*w^(a+2))/(a*(a+1)*(a+2))+(4^(a+1)*w^(a+1))/(a*(a+1))
3838                                                                      +(4^a*w^a)/a)
3839                                 *%e^-(4*w))
3840                                 +((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)
3841                                  *%e^-(2*v)+gamma_incomplete_generalized(a,2*v,4*w))); /* FOR SURE ?? */
3842 true;
3843
3844 should_be_equal (i(2*v, 4*w), subst ([x = 2*v, y = 4*w], i(x, y)));
3845 true;
3846
3847 should_be_equal (j(u), integrate (t^3*exp(-t),t,u,inf)/gamma(4));
3848 true;
3849
3850 should_be_equal (j(1 + u), ((u+1)^3/6+(u+1)^2/2+u+2)*%e^((-u)-1));
3851 true;
3852
3853 should_be_equal (j(1 + u), subst (v = 1 + u, j(v)));
3854 true;
3855
3856 block ([foo, bar, l:[0.1, 0.25, 0.5, 0.8, 0.9, 0.95]],
3857        foo : makelist (k(4, x), x, l),
3858        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),
3859        if lmax (abs (foo - bar)) < 1e-8 then true else notequal (foo, bar));
3860 true;
3861
3862 k(1 - w, 2*v);
3863 1-(1-2*v)^(1-w)*(2*v*(1-w)+1); /* assuming k(v, w) = 1-(1-w)^v*(v*w+1) */
3864
3865 should_be_equal (k(1 - v, 2*w), subst ([x = 1 - v, y = 2*w], k(x, y)));
3866 true;
3867
3868 reset (gamma_expand);
3869 [gamma_expand];