reduce-historical/r36/XMPL/SPECFN.TST

%
%  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;