Artifact d6da60a18c73e4e005d4e192a77f52a734826e023d8c4af6110ea0bb292a86cf:
- Executable file
r37/packages/specfn/specfn.tst
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 12495) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/specfn.tst
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 12495) [annotate] [blame] [check-ins using]
% % Testing file for REDUCE Special Functions Package % % Chris Cannam, ZIB Berlin % October 1992 -> Feb 1993 % (only some of the time, of course) % % Corrections and comments to neun@sc.zib-berlin.de % on savesfs; % just in case it's off for some reason off bfspace; % to provide more similarity between runs % with old & new bigfloats let {sinh (~x) => (exp(x) - exp (-x))/2 }; % this will improve some results % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 1. Bernoulli numbers % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= off rounded; procedure do!*one!*bern(x); write ("Bernoulli ", x, " is ", bernoulli x); do!*one!*bern(1); do!*one!*bern(2); do!*one!*bern(3); do!*one!*bern(13); do!*one!*bern(14); do!*one!*bern(300); do!*one!*bern(-2); do!*one!*bern(0); for n := 2 step 2 until 100 do do!*one!*bern n; on rounded; precision 100; do!*one!*bern(1); do!*one!*bern(2); do!*one!*bern(3); do!*one!*bern(13); do!*one!*bern(14); do!*one!*bern(300); do!*one!*bern(-2); do!*one!*bern(0); do!*one!*bern(38); do!*one!*bern(400); % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 2. Gamma function % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= on rounded; off complex; precision 40; algebraic procedure wg(x); write ("gamma (", x, ") ==> ", gamma x); algebraic procedure wp(x); write ("-- precision ", x, ", from ", precision(x)); wg (1/2); wg (3/2); write ("sqrt(pi)/2 ==> ", sqrt(pi)/2); wp(10); for x := 0 step 5 until 100 do << wg (1 + x/1000); wg (-1 - x/13); wp (8+floor(x/4)) >>; wg(1/1000000003); off rounded; gamma(17/2); gamma(-17/2); gamma(4); gamma(0); gamma(-4); gamma(-17/3); p := gamma(x**2) * gamma(x-y**gamma(y)) - (1/(gamma(4*(x-y)))); y := 1/4; p; x := 3; p; y := -3/8; p; on rounded, complex; precision 50; p; off rounded, complex; clear y; p; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 3. Beta function. Not very interesting % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= algebraic procedure do!*one!*beta(x,y); write ("Beta (", x, ",", y, ") = ", beta(x,y)); do!*one!*beta(0,1); do!*one!*beta(2,-3); do!*one!*beta(3,2); do!*one!*beta(a+b,(c+d)**(b-a)); do!*one!*beta(-3,4); do!*one!*beta(-3,2); do!*one!*beta(-3,-7.5); do!*one!*beta((pi * 10), exp(5)); on rounded; precision 30; do!*one!*beta(0,1); do!*one!*beta(2,-3); do!*one!*beta(3,2); do!*one!*beta(a+b,(c+d)**(b-a)); do!*one!*beta(-3,4); do!*one!*beta(-3,2); do!*one!*beta(-3,-7.5); do!*one!*beta((pi * 10), exp(5)); % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 4. Pochhammer notation % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= off rounded; pochhammer(4,5); pochhammer(-4,5); pochhammer(4,-5); pochhammer(-4,-5); pochhammer(17/2,12); pochhammer(-17/2,12); pochhammer(1/3,14)*pochhammer(2/3,15); q := pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)* pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)* pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11); on complex; pochhammer(a+b*i,c)*pochhammer(a-b*i,c); a := 2; b := 3; c := 5; pochhammer(a+b*i,c)*pochhammer(a-b*i,c); off complex; on rounded; pochhammer(1/5,11)*pochhammer(2/5,11)*pochhammer(3/5,11)* pochhammer(1-1/5,11)*pochhammer(1,11)*pochhammer(6/5,11)* pochhammer(70/50,11)*pochhammer(8/5,11)*pochhammer(9/5,11); q; pochhammer(pi,floor (pi**8)); pochhammer(-pi,floor (pi**7)); pochhammer(1.5,floor (pi**8)); % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 5. Digamma function % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= procedure do!*one!*psi(x); << precision (precision(0) + 4)$ write("Psi of ", x, " is ", psi(x) ) >> ; clear x, y; z := x * ((x+y)**2 + (x**y)); off rounded; do!*one!*psi(3); do!*one!*psi(pi); do!*one!*psi(1.005); do!*one!*psi(1.995); do!*one!*psi(74); do!*one!*psi(-1/2); do!*one!*psi(-3); do!*one!*psi(z); on rounded; precision 100; do!*one!*psi(3); do!*one!*psi(pi); do!*one!*psi(1.005); do!*one!*psi(1.995); do!*one!*psi(74); do!*one!*psi(-1/2); do!*one!*psi(-3); do!*one!*psi(z); precision 15; x := 8/3; y := 7/1000; do!*one!*psi(z); off rounded; clear x, y; df(psi(z), x); df(df(psi(z), y),x); int(psi(z), z); on rounded; for k := 1 step 0.1 until 2 do do!*one!*psi(k); off rounded; % PSI_SIMP.TST F.J.Wright, 2 July 1993 on evallhseqp; factor psi; on rat, intstr, div; % for neater output % Do not try using "off mcd"! psi(x+m) - psi(x+m-1) = 1/(x+m-1); psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x); psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1); psi(x + 1) = psi(x) + 1/x; psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2); psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2); psi((x + 3a)/a); psi(x/y + 3); off rat, intstr, div; on rational; psi(x+m) - psi(x+m-1) = 1/(x+m-1); psi(x+2) - psi(x+1) + 2*psi(x) = 1/(x+1) + 2*psi(x); psi(x+2) + 3*psi(x) = 4*psi(x) + 1/x + 1/(x+1); psi(x + 1) = psi(x) + 1/x; psi(x + 3/2) = psi(x + 1/2) + 1/(x + 1/2); psi(x - 1/2) = psi(x + 1/2) - 1/(x - 1/2); psi((x + 3a)/a); psi(x/y + 3); off rational; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 6. Polygamma functions % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= procedure do!*one!*pg(n,x); write ("Polygamma (", n, ") of ", x, " is ", polygamma(n,x)); off rounded; do!*one!*pg(1,1/2); do!*one!*pg(1,1); do!*one!*pg(1,3/2); do!*one!*pg(1,1.005); do!*one!*pg(1,1.995); do!*one!*pg(1,1e-10); do!*one!*pg(2,1.45); do!*one!*pg(3,1.99); do!*one!*pg(4,-8.2); do!*one!*pg(5,0); do!*one!*pg(6,-5); do!*one!*pg(7,200); on rounded; precision 100; do!*one!*pg(1,1/2); do!*one!*pg(1,1); do!*one!*pg(1,3/2); do!*one!*pg(1,1.005); do!*one!*pg(1,1.995); do!*one!*pg(1,1e-10); do!*one!*pg(2,1.45); do!*one!*pg(3,1.99); do!*one!*pg(4,-8.2); do!*one!*pg(5,0); do!*one!*pg(6,-5); do!*one!*pg(7,200); off rounded; clear x; % Polygamma differentiation has already % been tried a bit in the psi section df(int(int(int(polygamma(3,x),x),x),x),x); clear w, y, z; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 7. Zeta function % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= procedure do!*one!*zeta(n); write ("Zeta of ", n, " is ", zeta n); off rounded; clear x, y, z; z := x * ((x+y)**5 + (x**y)); do!*one!*zeta(0); for k := 4 step 2 until 35 do do!*one!*zeta(k); do!*one!*zeta(-17/3); do!*one!*zeta(190); do!*one!*zeta(300); do!*one!*zeta(0); do!*one!*zeta(-44); on rounded; clear x, y; for k := 3 step 3 until 36 do << precision (31+k*3); do!*one!*zeta(k) >>; precision 20; do!*one!*zeta(-17/3); do!*one!*zeta(z); y := 3; x := pi; do!*one!*zeta(z); do!*one!*zeta(190); do!*one!*zeta(300); do!*one!*zeta(0); do!*one!*zeta(-44); off rounded; % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 8. Kummer functions % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= off rounded; t!*kummer!*a := { 1, 2.4, -1397/10 }$ t!*kummer!*b := { 0, 1, pi, -pi, 26 }$ for each a in t!*kummer!*a do for each b in t!*kummer!*a do for each z in t!*kummer!*a do << write "KummerM(", a, ",", b, ",", z, ") = ", kummerm(a,b,z); write "KummerU(", a, ",", b, ",", z, ") = ", kummeru(a,b,z) >>; on rounded; precision 30; t!*k!*c := 7; % To test each and every possible combination of % three arguments from t!*kummer!*b would take too % long, but we want the possibility of trying most % special cases. Compromise: test every seventh % possibility. for each a in t!*kummer!*b do for each b in t!*kummer!*b do for each z in t!*kummer!*b do << if t!*k!*c = 7 then << write "KummerM(", a, ",", b, ",", z, ") = ", kummerm(a,b,z); write "KummerU(", a, ",", b, ",", z, ") = ", kummeru(a,b,z); t!*k!*c := 0 >>; t!*k!*c := t!*k!*c + 1 >>; off rounded; clear x, y, z, t!*k!*c; df(df(kummerM(a,b,z),z),z); df(kummerU(a,b,z),z); z := ((x^2 + y)^5) + (x^(x+y)); df(df(kummerM(a,b,z),y),x); df(kummerU(a,b,z),x); % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % 9. Bessel functions % =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= % Lengthy test of the Bessel functions. This isn't even % remotely exhaustive of the special cases -- though a % real person with lots of time could no doubt come up % with a better lot of tests than this automated rubbish. % Again, compromise by only actually doing one in five or % nine. If you want a really thorough test, you can % easily change this to get it; but it'll take hours to % run. clear p, q; hankel1(p,q); r := df(ws,q); on complex; r; p:=3/8; r; q := pi; r; on rounded; r; off complex, rounded; df(df(besselj(pp,qq)+rr * hankel1(pp*2,qq) * bessely(pp-qq,qq),qq),qq); % Possible values for real args t!*bes!*vr := { 1, pi, -pi, 26 }$ % Possible values for real and imaginary parts of complex args t!*bes!*vc := { 0, 3, -41/2 }$ array s!*bes(4)$ s!*bes(1) := "BesselJ"$ s!*bes(2) := "BesselY"$ s!*bes(3) := "BesselI"$ s!*bes(4) := "BesselK"$ pre := 16; precision pre; preord := 10**pre; t!*b!*c := 3; algebraic procedure do!*one!*bessel(s,n,z); (if s = 1 then besselj(n,z) else if s = 2 then bessely(n,z) else if s = 3 then besseli(n,z) else besselk(n,z)); algebraic procedure pr!*bessel(s,n,z,k); << if t!*b!*c = k then << on rounded; bes1 := do!*one!*bessel(s,n,z); precision(pre+5); bes2 := do!*one!*bessel(s,n,z); if bes1 neq 0 then disc := floor abs(100*(bes2-bes1)*preord/bes1) else disc := 0; precision pre; write s!*bes(s), "(", n, ",", z, ") = ", bes1; if not numberp disc then << precom := !*complex; on complex; disc := disc; if not precom then off complex >>; if disc neq 0 then write " (discrepancy ", disc, "% of one s.f.)"; if numberp disc and disc > 200 then << write "***** WARNING Significant Inaccuracy."; write " Lower precision result:"; write " ", bes1; write " Higher precision result:"; precision(pre+5); write " ", bes2; precision pre >>; off rounded; t!*b!*c := 0 >>; t!*b!*c := t!*b!*c + 1 >>; % About to begin Bessel test. We have a list of possible % values, and we test every Bessel, with every value on the % list as both order and argument. Every Bessel is computed % twice, to different precisions (differing by 3), and any % discrepancy is reported. The value reported is the diff- % erence between the two computed values, expressed as a % percentage of the unit of the least significant displayed % digit. A discrepancy between 100% and 200% means that the % last digit of the displayed value was found to differ at % higher precision; values greater than 200% are cause for % concern. An ideal discrepancy would be between about 1% % and 20%; if the value is found to be zero, no discrepancy % is reported. off msg; for s := 1:4 do << write(" ... Testing ", s!*bes(s), " for real domains ... "); for each n in t!*bes!*vr do for each z in t!*bes!*vr do pr!*bessel(s, n, z, 5) >>; on complex; for s := 1:3 do << write (" ... Testing ", s!*bes(s), " for complex domains ... "); for each nr in t!*bes!*vc do for each ni in t!*bes!*vc do for each zr in t!*bes!*vc do for each zi in t!*bes!*vc do pr!*bessel(s, nr+ni*i, zr+zi*i, 9) >>; off complex; on msg; write (" ..."); write ("Bessel test complete."); end;