1% 2% Testing file for REDUCE Special Functions Package 3% 4% Chris Cannam, ZIB Berlin 5% October 1992 -> Feb 1993 6% (only some of the time, of course) 7% 8% Corrections and comments to neun@sc.zib-berlin.de 9% 10 11 12on savesfs; % just in case it's off for some reason 13 14off bfspace; % to provide more similarity between runs 15 % with old & new bigfloats 16 17let {sinh (~x) => (exp(x) - exp (-x))/2 }; 18 % this will improve some results 19 20 21% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 22 23% 1. Bernoulli numbers 24 25% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 26 27 28off rounded; 29 30procedure do!*one!*bern(x); 31 write ("Bernoulli ", x, " is ", bernoulli x); 32 33do!*one!*bern(1); 34do!*one!*bern(2); 35do!*one!*bern(3); 36do!*one!*bern(13); 37do!*one!*bern(14); 38do!*one!*bern(300); 39do!*one!*bern(-2); 40do!*one!*bern(0); 41 42for n := 2 step 2 until 100 do do!*one!*bern n; 43 44on rounded; 45precision 100; 46 47do!*one!*bern(1); 48do!*one!*bern(2); 49do!*one!*bern(3); 50do!*one!*bern(13); 51do!*one!*bern(14); 52do!*one!*bern(300); 53do!*one!*bern(-2); 54do!*one!*bern(0); 55do!*one!*bern(38); 56do!*one!*bern(400); 57 58 59% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 60 61% 2. Gamma function 62 63% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 64 65 66on rounded; 67off complex; 68precision 40; 69 70algebraic procedure wg(x); 71 write ("gamma (", x, ") ==> ", gamma x); 72 73algebraic procedure wp(x); 74 write ("-- precision ", x, ", from ", precision(x)); 75 76wg (1/2); 77wg (3/2); 78 79write ("sqrt(pi)/2 ==> ", sqrt(pi)/2); 80 81wp(10); 82 83for x := 0 step 5 until 100 do 84 << wg (1 + x/1000); 85 wg (-1 - x/13); 86 wp (8+floor(x/4)) >>; 87 88wg(1/1000000003); 89 90off rounded; 91 92gamma(17/2); 93gamma(-17/2); 94gamma(4); 95gamma(0); 96gamma(-4); 97gamma(-17/3); 98 99p := gamma(x**2) * gamma(x-y**gamma(y)) - (1/(gamma(4*(x-y)))); 100y := 1/4; 101p; 102x := 3; 103p; 104y := -3/8; 105p; 106 107on rounded, complex; 108precision 50; 109 110p; 111 112off rounded, complex; 113clear y; 114 115p; 116 117 118 119% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 120 121% 3. Beta function. Not very interesting 122 123% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 124 125 126algebraic procedure do!*one!*beta(x,y); 127 write ("Beta (", x, ",", y, ") = ", beta(x,y)); 128 129do!*one!*beta(0,1); 130do!*one!*beta(2,-3); 131do!*one!*beta(3,2); 132do!*one!*beta(a+b,(c+d)**(b-a)); 133do!*one!*beta(-3,4); 134do!*one!*beta(-3,2); 135do!*one!*beta(-3,-7.5); 136do!*one!*beta((pi * 10), exp(5)); 137 138on rounded; 139precision 30; 140 141do!*one!*beta(0,1); 142do!*one!*beta(2,-3); 143do!*one!*beta(3,2); 144do!*one!*beta(a+b,(c+d)**(b-a)); 145do!*one!*beta(-3,4); 146do!*one!*beta(-3,2); 147do!*one!*beta(-3,-7.5); 148do!*one!*beta((pi * 10), exp(5)); 149 150 151 152% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 153 154% 4. Pochhammer notation 155 156% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 157 158 159off rounded; 160 161pochhammer(4,5); 162pochhammer(-4,5); 163pochhammer(4,-5); 164pochhammer(-4,-5); 165pochhammer(17/2,12); 166pochhammer(-17/2,12); 167pochhammer(1/3,14)*pochhammer(2/3,15); 168 169q := pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)* 170 pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)* 171 pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11); 172 173on complex; 174 175pochhammer(a+b*i,c)*pochhammer(a-b*i,c); 176 177a := 2; 178b := 3; 179c := 5; 180 181pochhammer(a+b*i,c)*pochhammer(a-b*i,c); 182 183off complex; 184on rounded; 185 186pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)* 187 pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)* 188 pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11); 189 190q; 191 192pochhammer(pi,floor (pi**8)); 193pochhammer(-pi,floor (pi**7)); 194pochhammer(1.5,floor (pi**8)); 195 196 197 198% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 199 200% 5. Digamma function 201 202% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 203 204 205procedure do!*one!*psi(x); 206 << precision (precision(0) + 4)$ 207 write("Psi of ", x, " is ", psi(x) ) >> ; 208 209clear x, y; 210 211z := x * ((x+y)**2 + (x**y)); 212 213off rounded; 214 215do!*one!*psi(3); 216do!*one!*psi(pi); 217do!*one!*psi(1.005); 218do!*one!*psi(1.995); 219do!*one!*psi(74); 220do!*one!*psi(-1/2); 221do!*one!*psi(-3); 222do!*one!*psi(z); 223 224on rounded; 225precision 100; 226 227do!*one!*psi(3); 228do!*one!*psi(pi); 229do!*one!*psi(1.005); 230do!*one!*psi(1.995); 231do!*one!*psi(74); 232do!*one!*psi(-1/2); 233do!*one!*psi(-3); 234do!*one!*psi(z); 235 236precision 15; 237 238x := 8/3; 239y := 7/1000; 240do!*one!*psi(z); 241 242off rounded; 243 244clear x, y; 245 246df(psi(z), x); 247df(df(psi(z), y),x); 248int(psi(z), z); 249 250on rounded; 251 252for k := 1 step 0.1 until 2 do do!*one!*psi(k); 253 254off rounded; 255 256% PSI_SIMP.TST F.J.Wright, 2 July 1993 257 258on evallhseqp; 259factor psi; on rat, intstr, div; % for neater output 260% Do not try using "off mcd"! 261 262psi(x+m) - psi(x+m-1) = 1/(x+m-1); 263psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x); 264psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1); 265 266psi(x + 1) = psi(x) + 1/x; 267psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2); 268psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2); 269psi((x + 3a)/a); psi(x/y + 3); 270 271off rat, intstr, div; on rational; 272 273psi(x+m) - psi(x+m-1) = 1/(x+m-1); 274psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x); 275psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1); 276 277psi(x + 1) = psi(x) + 1/x; 278psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2); 279psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2); 280psi((x + 3a)/a); psi(x/y + 3); 281 282off rational; 283 284% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 285 286% 6. Polygamma functions 287 288% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 289 290 291procedure do!*one!*pg(n,x); 292 write ("Polygamma (", n, ") of ", x, " is ", polygamma(n,x)); 293 294off rounded; 295 296do!*one!*pg(1,1/2); 297do!*one!*pg(1,1); 298do!*one!*pg(1,3/2); 299do!*one!*pg(1,1.005); 300do!*one!*pg(1,1.995); 301do!*one!*pg(1,1e-10); 302do!*one!*pg(2,1.45); 303do!*one!*pg(3,1.99); 304do!*one!*pg(4,-8.2); 305do!*one!*pg(5,0); 306do!*one!*pg(6,-5); 307do!*one!*pg(7,200); 308 309on rounded; 310precision 100; 311 312do!*one!*pg(1,1/2); 313do!*one!*pg(1,1); 314do!*one!*pg(1,3/2); 315do!*one!*pg(1,1.005); 316do!*one!*pg(1,1.995); 317do!*one!*pg(1,1e-10); 318do!*one!*pg(2,1.45); 319do!*one!*pg(3,1.99); 320do!*one!*pg(4,-8.2); 321do!*one!*pg(5,0); 322do!*one!*pg(6,-5); 323do!*one!*pg(7,200); 324 325off rounded; 326clear x; 327 328 329% Polygamma differentiation has already 330% been tried a bit in the psi section 331 332df(int(int(int(polygamma(3,x),x),x),x),x); 333 334clear w, y, z; 335 336 337 338% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 339 340% 7. Zeta function 341 342% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 343 344 345procedure do!*one!*zeta(n); 346 write ("Zeta of ", n, " is ", zeta n); 347 348off rounded; 349clear x, y, z; 350 351z := x * ((x+y)**5 + (x**y)); 352 353do!*one!*zeta(0); 354 355for k := 4 step 2 until 35 do 356 do!*one!*zeta(k); 357 358for k := 1 step 2 until 31 do 359 do!*one!*zeta(-k); 360 361do!*one!*zeta(-17/3); 362do!*one!*zeta(190); 363do!*one!*zeta(300); 364do!*one!*zeta(0); 365do!*one!*zeta(-44); 366 367on rounded; 368clear x, y; 369 370for k := 3 step 3 until 36 do << 371 precision (31+k*3); 372 do!*one!*zeta(k) >>; 373 374precision 20; 375 376do!*one!*zeta(-17/3); 377do!*one!*zeta(z); 378 379y := 3; 380x := pi; 381 382do!*one!*zeta(z); 383do!*one!*zeta(190); 384do!*one!*zeta(300); 385do!*one!*zeta(0); 386do!*one!*zeta(-44); 387 388off rounded; 389 390 391 392% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 393 394% 8. Kummer functions 395 396% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 397 398 399off rounded; 400 401t!*kummer!*a := { 1, 2.4, -1397/10 }$ 402t!*kummer!*b := { 0, 1, pi, -pi, 26 }$ 403 404for each a in t!*kummer!*a do 405 for each b in t!*kummer!*a do 406 for each z in t!*kummer!*a do 407 << write "KummerM(", a, ",", b, ",", z, ") = ", 408 kummerm(a,b,z); 409 write "KummerU(", a, ",", b, ",", z, ") = ", 410 kummeru(a,b,z) >>; 411 412on rounded; 413precision 30; 414t!*k!*c := 7; 415 416% To test each and every possible combination of 417% three arguments from t!*kummer!*b would take too 418% long, but we want the possibility of trying most 419% special cases. Compromise: test every seventh 420% possibility. 421 422for each a in t!*kummer!*b do 423 for each b in t!*kummer!*b do 424 for each z in t!*kummer!*b do 425 << if t!*k!*c = 7 426 then << write "KummerM(", a, ",", b, ",", z, ") = ", 427 kummerm(a,b,z); 428 write "KummerU(", a, ",", b, ",", z, ") = ", 429 kummeru(a,b,z); 430 t!*k!*c := 0 >>; 431 t!*k!*c := t!*k!*c + 1 >>; 432 433off rounded; 434clear x, y, z, t!*k!*c; 435 436df(df(kummerM(a,b,z),z),z); 437df(kummerU(a,b,z),z); 438 439z := ((x^2 + y)^5) + (x^(x+y)); 440 441df(df(kummerM(a,b,z),y),x); 442df(kummerU(a,b,z),x); 443 444 445 446% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 447 448% 9. Bessel functions 449 450% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 451 452 453% Lengthy test of the Bessel functions. This isn't even 454% remotely exhaustive of the special cases -- though a 455% real person with lots of time could no doubt come up 456% with a better lot of tests than this automated rubbish. 457% Again, compromise by only actually doing one in five or 458% nine. If you want a really thorough test, you can 459% easily change this to get it; but it'll take hours to 460% run. 461 462clear p, q; 463 464hankel1(p,q); 465r := df(ws,q); 466 467on complex; 468 469r; 470p:=3/8; 471r; 472q := pi; 473r; 474 475on rounded; 476 477r; 478 479off complex, rounded; 480 481df(df(besselj(pp,qq)+rr * hankel1(pp*2,qq) * 482 bessely(pp-qq,qq),qq),qq); 483 484% Possible values for real args 485t!*bes!*vr := { 1, pi, -pi, 26 }$ 486 487% Possible values for real and imaginary parts of complex args 488t!*bes!*vc := { 0, 3, -41/2 }$ 489 490array s!*bes(4)$ 491 492s!*bes(1) := "BesselJ"$ 493s!*bes(2) := "BesselY"$ 494s!*bes(3) := "BesselI"$ 495s!*bes(4) := "BesselK"$ 496 497pre := 16; 498precision pre; 499preord := 10**pre; 500t!*b!*c := 3; 501 502algebraic procedure do!*one!*bessel(s,n,z); 503 (if s = 1 then besselj(n,z) 504 else if s = 2 then bessely(n,z) 505 else if s = 3 then besseli(n,z) 506 else besselk(n,z)); 507 508algebraic procedure pr!*bessel(s,n,z,k); 509 << if t!*b!*c = k 510 then 511 << on rounded; 512 bes1 := do!*one!*bessel(s,n,z); 513 precision(pre+5); 514 bes2 := do!*one!*bessel(s,n,z); 515 if bes1 neq 0 516 then disc := floor abs(100*(bes2-bes1)*preord/bes1) 517 else disc := 0; 518 precision pre; 519 write s!*bes(s), "(", n, ",", z, ") = ", bes1; 520 if not numberp disc then 521 << precom := !*complex; 522 on complex; 523 disc := disc; 524 if not precom then off complex >>; 525 if disc neq 0 526 then write " (discrepancy ", disc, "% of one s.f.)"; 527 if numberp disc and disc > 200 then 528 << write "***** WARNING Significant Inaccuracy."; 529 write " Lower precision result:"; 530 write " ", bes1; 531 write " Higher precision result:"; 532 precision(pre+5); write " ", bes2; precision pre >>; 533 off rounded; 534 t!*b!*c := 0 >>; 535 t!*b!*c := t!*b!*c + 1 >>; 536 537% About to begin Bessel test. We have a list of possible 538% values, and we test every Bessel, with every value on the 539% list as both order and argument. Every Bessel is computed 540% twice, to different precisions (differing by 3), and any 541% discrepancy is reported. The value reported is the diff- 542% erence between the two computed values, expressed as a 543% percentage of the unit of the least significant displayed 544% digit. A discrepancy between 100% and 200% means that the 545% last digit of the displayed value was found to differ at 546% higher precision; values greater than 200% are cause for 547% concern. An ideal discrepancy would be between about 1% 548% and 20%; if the value is found to be zero, no discrepancy 549% is reported. 550 551off msg; 552 553for s := 1:4 do 554 << write(" ... Testing ", s!*bes(s), " for real domains ... "); 555 for each n in t!*bes!*vr do 556 for each z in t!*bes!*vr do 557 pr!*bessel(s, n, z, 5) >>; 558 559on complex; 560 561for s := 1:3 do 562 << write (" ... Testing ", s!*bes(s), " for complex domains ... "); 563 for each nr in t!*bes!*vc do 564 for each ni in t!*bes!*vc do 565 for each zr in t!*bes!*vc do 566 for each zi in t!*bes!*vc do 567 pr!*bessel(s, nr+ni*i, zr+zi*i, 9) >>; 568 569off complex; 570 571on msg; 572 573write (" ..."); 574write ("Bessel test complete."); 575 576 577% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 578 579% 10. Incomplete Gamma and Beta functions (regularized) 580 581% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 582 583 584igamma(3,0); 585 586igamma(1,1); 587 588igamma(1,4); 589 590on rounded; 591 592igamma(1,1); 593 594igamma(1,4); 595 596igamma(2,4); 597 598igamma(0.5,4); 599 600beta(1,1,1); 601 602ibeta(1,2,1); 603 604ibeta(1,4,1); 605 606ibeta(2,4,1); 607 608ibeta(2,4,0.5); 609 610ibeta(0.12,0.43,0.9); 611 612precision 50; 613 614ibeta(0.12,0.43,0.9); 615 616precision 20; 617 618ibeta(0.12,0.43,0.9); 619 620on complex; 621 622ibeta(1+i,1,1.5*i); 623 624off rounded,complex; 625 626ibeta(3,2,x); 627 628 629% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 630 631% 11. Dilogarithm, polylogartihm and Lerch_phi 632 633% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 634 635 636polylog(n,0); 637 638polylog(2,1); 639 640polylog(3,1); 641 642polylog(2,i); 643 644df(polylog(a,x),x); 645 646polylog(1,x); 647 648precision reset; 649on rounded; 650 651polylog(2,1/3); 652 653off rounded; 654 655 656lerch_phi(3,4,1); 657 658lerch_phi(4,0,3); 659 660lerch_phi(x,a,1); 661 662lerch_phi(1,x,1); 663 664df(lerch_phi(x,3,4),x); 665 666 667% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 668 669% 12. Constants 670 671% =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 672 673on rounded; 674off complex; 675 676precision 50; 677 678euler_gamma; 679 680Khinchin; 681 682golden_ratio; 683 684Catalan; 685 686on complex; 687 688euler_gamma; 689 690Khinchin; 691 692golden_ratio; 693 694Catalan; 695 696 697end; 698 699 700