% Floating point tests. A C Norman and Codemist Ltd. 1994
%
% The code here is for use by people who have read and
% enjoyed "Software manual for the elementary functions" by
% Cody and Waite (Prentice Hall, 1980). It attempts to
% measure the accuracy of some of the elementary functions
% in the implementation of REDUCE that it is run on. Good
% implementations should show maximum errors of only a few
% units in the last place, and errors reasonably evenly
% distributed between +ve and -ve. Various special cases are
% tested too. The layout for presenting results from this code
% is somewhat crude at present.
%
% call fptest() to perform complete set of tests;
symbolic;
on comp;
fluid '(ibeta it irnd ngrd machep negep iexp minexp
maxexp eps epsneg xmin xmax iy one zero x);
symbolic procedure fptest();
begin
scalar ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,
maxexp,eps,epsneg,xmin,xmax;
iy:=100001; % seed for random number gereration;
machar();
% only expt is left after this lot. hurrah;
test!-arcsin();
test!-atan();
test!-tan();
test!-sin();
test!-sqrt();
test!-log();
test!-exp();
princ "fp tests over";
terpri()
end;
symbolic procedure machar();
begin
scalar one,zero,a,b,beta,betam1,betain,i,k,x,y,z,mx;
one:=float 1;
zero:=float 0;
a:=one;
repeat a:=a+a until not (((a+one)-a)-one = zero);
b:=one;
repeat b:=b+b until not ((a+b)-a = zero);
ibeta:=fix((a+b)-a);
princ list("ibeta = ",ibeta); terpri();
beta:=float ibeta;
it:=0;
b:=one;
repeat << it:=it+1; b:=b*beta >> until not (((b+one)-b)-one = zero);
irnd:=0;
betam1:=beta-one;
if ((a+betam1)-a) neq zero then irnd:=1;
princ list("it=",it,"irnd=",irnd); terpri();
negep:=it+3;
betain:=one/beta;
a:=one;
for i:=1:negep do a:=a*betain;
b:=a;
while ((one-a)-one) = zero do << a:=a*beta; negep:=negep-1 >>;
negep:=-negep;
epsneg:=a;
if not ((ibeta=2) or (irnd=0)) then <<
a := (a*(one+a))/(one+one);
if ((one-a)-one) neq zero then epsneg:=a >>;
princ list("negep=",negep,"epsneg=",epsneg); terpri();
machep:=-it-3;
a:=b;
while ((one+a)-one) = zero do <<
a:=a*beta;
machep:=machep+1 >>;
eps:=a;
if not (ibeta=2 or irnd=0) then <<
a:=(a*(one+a))/(one+one);
if (one+a)-one neq zero then eps:=a >>;
princ list("machep=",machep,"eps=",eps); terpri();
ngrd:=0;
if irnd=0 and ((one+eps)*one-one) neq zero then ngrd:=1;
princ list("ngrd=",ngrd); terpri();
i:=0;
k:=1;
z:=betain;
l400: y:=z;
z:=y*y;
a:=z*one;
if (a+a=zero) or (abs z >= y) then go to l410;
i:=i+1;
k:=k+k;
goto l400;
l410: % assume it is not a decimal machine! ;
iexp:=i+1;
mx:=k+k;
l450: xmin:=y;
y:=y*betain;
a:=y*one;
if (a+a=zero) or (abs y >= xmin) or
((a * (one + eps)) = a) then goto l460;
k:=k+1;
goto l450;
l460: minexp:=-k;
princ list("iexp=",iexp,"minexp=",minexp,"xmin=",xmin); terpri();
if not (mx>k+k-3 or ibeta=10) then <<
mx:=mx+mx;
iexp:=iexp+1 >>;
maxexp:=mx+minexp;
maxexp:=maxexp-2; % allow for 2 reserved exponents in ieee format;
i:=maxexp+minexp;
if (ibeta=2) and (i=0) then maxexp:=maxexp-1;
if i>20 then maxexp:=maxexp-1;
if a neq y then maxexp:=maxexp-2;
princ list("maxexp=",maxexp); terpri();
xmax:=one-epsneg;
if xmax*one neq xmax then xmax:=one-beta*epsneg;
xmax:=xmax/(beta*beta*beta*xmin);
i:=maxexp+minexp+3;
if i>0 then <<
for j:=1:i do
if ibeta=2 then xmax:=xmax+xmax else xmax:=xmax*beta >>;
princ list("xmax=",xmax); terpri();
princ "end of machar"; terpri()
end;
symbolic procedure rnd();
% random floating point number in range 0 to 1;
begin
iy:=iy*125;
iy:=remainder(iy,2796203);
return float(iy)/2796203.0
end;
symbolic procedure randl(x);
exp(x*rnd());
symbolic procedure tab_to_column n;
while posn() < n do princ " ";
symbolic procedure heading(name,trials,from,to);
<< terpri();
princ name;
tab_to_column 20; princ trials; princ " trials";
tab_to_column 40; princ "from "; princ from;
princ " to "; princ to; terpri()
>>;
symbolic procedure errsigns(k1,k2,k3);
<< princ "+ve: "; princ k1;
tab_to_column 20; princ "zero: "; princ k2;
tab_to_column 40; princ "-ve: "; princ k3; terpri()
>>;
symbolic procedure maxrel(w,x1,d);
<< princ "mre was "; princ w;
princ " at "; princ x1;
princ " lost "; princ d; princ " digits"; terpri();
>>;
symbolic procedure rmserr(w,d);
<< princ "rms was "; princ w;
princ " lost "; princ d; princ " digits"; terpri()
>>;
symbolic procedure specials(a,b,c);
<< princ a;
princ b;
princ " ";
princ c; terpri()
>>;
symbolic procedure errortest u;
begin
scalar w;
w := errorset(u, t, t);
if atom w then <<
princ "Evaluation of "; prin1 u; princ " failed"; terpri() >>
else <<
prin1 u; princ " = "; print car w >>;
end;
symbolic procedure test!-arcsin();
begin
scalar a,ait,albeta,b,beta,betap,c1,c2,del,half,pi,r6,r7,one,xsq,two,l,
k,em,ob32,w,x,xl,xn,x1,y,z,zero,zz,i,i1,j,k1,k2,k3,n,expon,zsum,
ysq,xm,s;
beta:=float ibeta;
albeta:=log beta;
ait:=float it;
zero:=float 0;
one:=float 1;
half:=0.5;
two:=float 2;
k:=fix(log10(beta**it))+1;
c1:=201.0/128.0;
c2:=4.8382679489661923132e-4;
a:=-0.125;
b:=-a;
n:=2000;
xn:=float n;
i1:=0;
l:=-1;
for j:=1:5 do <<
princ list("start of test ",j," of arcsin/arccos"); terpri();
k1:=0;
k3:=0;
l:=-l;
x1:=zero;
r6:=zero;
r7:=zero;
del:=(b-a)/xn;
xl:=a;
for i:=1:n do <<
x:=del*rnd()+xl;
if j>2 then <<
ysq := half - half*abs x;
x:=(half-(ysq+ysq)) + half;
if j=5 then x:=-x;
y:=sqrt(ysq);
y:=y+y >>
else << y:=x; ysq:=y*y >>;
zsum:=zero;
xm:=float(k+k+1);
if l>0 then z:=asin x else z:=acos x;
for m:=1:k do <<
zsum:=ysq*(zsum+1.0/xm);
xm:=xm-2.0;
zsum:=zsum*(xm/(xm+1.0)) >>;
zsum:=zsum*y;
if j=1 or j=4 then <<
zz:=y+zsum;
zsum:=(y-zz)+zsum;
if irnd neq 1 then zz:=zz + (zsum+zsum) >>
else <<
s:=c1+c2;
zsum:=((c1-s) + c2) - zsum;
zz:=s+zsum;
zsum:=((s-zz)+zsum)-y;
s:=zz;
zz:=s+zsum;
zsum:=(s-zz)+zsum;
if irnd neq 1 then zz:=zz+(zsum+zsum)>>;
w:=1.0;
if z neq zero then w:=(z-zz)/z;
if w>zero then k1:=k1+1;
if w<zero then k3:=k3+1;
w:=abs w;
if w>r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w;
xl:=xl+del >>;
k2:=n-k1-k3;
r7:=sqrt(r7/xn);
heading("arcsin/arccos",n,a,b);
errsigns(k1,k2,k3);
w:=-999.0;
if r6 neq zero then w:=log abs r6/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log abs r7/albeta;
rmserr(w,max(ait+w,zero));
if j=2 then << a:=0.75; b:=1.0 >>
else if j=4 then << b:=-a; a:=-1.0; c1:=c1+c1; c2:=c2+c2; l:=-l >> >>;
for i:=1:5 do <<
x:=rnd()*a;
z:=asin x + asin(-x);
specials("asin(x)+asin(-x) ",x,z) >>;
betap:=beta**it;
x:=rnd()/betap;
for i:=1:5 do <<
z:=x-asin x;
specials("small args to asin ",x,z);
x := x/beta >>;
expon:=float minexp * 0.75;
x:=beta**expon;
y:=asin x;
specials("tiny arg ",x,y);
x:=1.2;
errortest list('asin, x);
princ "end of test on arcsin/arccos"; terpri()
end;
symbolic procedure test!-atan();
begin
scalar a,ait,albeta,b,beta,betap,c1,c2,del,half,pi,r6,r7,one,xsq,two,
em,ob32,w,x,xl,xn,x1,y,z,zero,zz,i,i1,j,k1,k2,k3,n,expon,zsum;
beta:=float ibeta;
albeta:=log beta;
ait:=float it;
zero:=float 0;
one:=float 1;
half:=0.5;
two:=float 2;
a:=-0.0625;
b:=-a;
ob32:=b*half;
n:=2000;
xn:=float n;
i1:=0;
for j:=1:4 do <<
princ list("start of test ",j," of atan"); terpri();
k1:=0;
k3:=0;
x1:=zero;
r6:=zero;
r7:=zero;
del:=(b-a)/xn;
xl:=a;
for i:=1:n do <<
x:=del*rnd()+xl;
if j=2 then x:=((one+x*a)-one)*16.0;
z:=atan x;
if j=1 then <<
xsq:=x*x;
em:=17.0;
zsum:=xsq/em;
for i1:=1:7 do <<
em := em-two;
zsum:=(one/em - zsum)*xsq >>;
zsum:=-x*zsum;
zz:=x+zsum;
zsum:=(x-zz)+zsum;
if irnd=0 then zz:=zz+(zsum+zsum) >>
else if j=2 then <<
y:=x-0.0625;
y:=y/(one+x*a);
zz:=(atan y - 8.1190004042651526021/100000.0)+ob32;
zz:=zz+ob32 >>
else <<
z:=z+z;
y:=x/((half+x*half)*((half-x)+half));
zz := atan y >>;
w:=1.0;
if z neq zero then w:=(z-zz)/z;
if w>zero then k1:=k1+1;
if w<zero then k3:=k3+1;
w:=abs w;
if w>r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w;
xl:=xl+del >>;
k2:=n-k1-k3;
r7:=sqrt(r7/xn);
heading("atan",n,a,b);
errsigns(k1,k2,k3);
w:=-999.0;
if r6 neq zero then w:=log abs r6/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log abs r7/albeta;
rmserr(w,max(ait+w,zero));
a:=b;
if j=1 then b:=two-sqrt(3.0);
if j=2 then b:=sqrt(two)-one;
if j=3 then b:=one >>;
a:=5.0;
for i:=1:5 do <<
x:=rnd()*a;
z:=atan x + atan(-x);
specials("atan(x)+atan(-x) ",x,z) >>;
betap:=beta**it;
x:=rnd()/betap;
for i:=1:5 do <<
z:=x-atan x;
specials("small args to atan ",x,z);
x := x/beta >>;
a:=-two;
b:=4.0;
for i:=1:5 do <<
x:=rnd()*b+a;
y:=rnd();
w:=-y;
z:=atan(x/y) - atan2(x,y);
zz:=atan(x/w) - atan2(x,w);
princ list("atan vs atan2 ",x,y,w," ",z,zz); terpri() >>;
expon:=float minexp * 0.75;
x:=beta**expon;
y:=atan x;
specials("tiny arg ",x,y);
specials("xmax ",xmax,atan xmax);
specials("atan(xmax,xmin) ",xmax,atan2(xmax,xmin));
princ "end of test on atan"; terpri();
end;
symbolic procedure test!-sqrt();
begin
scalar beta,sqbeta,albeta,ait,one,zero,n,xn,c,k1,k2,k3,
x1,r6,r7,a,b,x,y,z,w,w1;
princ "start of tests on sqrt"; terpri();
beta:=float(ibeta);
sqbeta:=sqrt(beta);
albeta:=log(beta);
ait:=float(it);
one:=float(1);
zero:=float(0);
a:=one/sqbeta;
b:=one;
n:=2000;
xn:=float(n);
for j:=1:2 do <<
c:=log(b/a);
k1:=0;
k3:=0;
x1:=zero;
r6:=zero;
r7:=zero;
for i:=1:n do <<
x:=a*randl(c);
y:=x*x;
z:=sqrt(y);
w:=(z-x)/x;
if w>zero then k1:=k1+1;
if w<zero then k3:=k3+1;
w:=abs w;
if w>=r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w >>;
k2:=n-k1-k3;
heading("sqrt",n,a,b);
errsigns(k1,k2,k3);
r7:=sqrt(r7/xn);
w:=-999.0;
if r6 neq zero then w:=log(abs r6)/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log(r7)/albeta;
rmserr(w,max(ait+w,zero));
a:=one;
b:=sqbeta >>;
specials("sqrt(xmin)",xmin,sqrt(xmin));
specials("1-epsneg ",one-epsneg,sqrt(one-epsneg));
specials("1 ",one,sqrt(one));
specials("1+eps ",one+eps,sqrt(one+eps));
specials("xmax ",xmax,sqrt(xmax));
specials("zero ",zero,sqrt(zero));
errortest list('sqrt, - one);
end;
symbolic procedure sign(a,b);
if minusp b then if minusp a then a else -a
else if minusp a then -a else a;
symbolic procedure test!-log();
begin
scalar a,ait,albeta,b,beta,c,del,eight,half,one,
r6,r7,tenth,w,x,xl,xn,x1,y,z,zero,zz,i,i1,j,k1,k2,k3,n;
princ "start of tests on log"; terpri();
beta:=float ibeta;
albeta:=log beta;
ait:=float it;
j:=it/3;
zero:=float 0;
one:=float 1;
eight:=float 8;
half:=one/float 2;
tenth:=one/float 10;
c:=one;
for i:=1:j do c:=c/beta;
b:=one+c;
a:=one-c;
n:=2000;
xn:=float n;
i1:=0;
for j:=1:4 do <<
princ list("start of test ",j); terpri();
k1:=0;
k3:=0;
x1:=zero;
r6:=zero;
r7:=zero;
del:=(b-a)/xn;
xl:=a;
for i:=1:n do <<
x:=del*rnd()+xl;
if j=1 then <<
y:=(x-half)-half;
zz:=log x;
z:=one/3.0;
z:=y*(z-y/4.0);
z:=(z-half)*y*y+y >>
else if j=2 then <<
x:=(x+eight)-eight;
y:=x+x/16.0;
z:=log x;
zz:=log y - 7.7746816434842581/100000.0;
zz:=zz-31.0/512.0 >>
else if j=3 then <<
x:=(x+eight)-eight;
y:=x+x*tenth;
z:=log10 x;
zz:=log10 y-3.7706015822504075/10000.0;
zz:=zz-21.0/512.0 >>
else <<
z:=log(x*x);
zz:=log x;
zz:=zz+zz >>;
w:=one;
if z neq zero then w:=(z-zz)/z;
z:=sign(w,z);
if z>zero then k1:=k1+1;
if z<zero then k3:=k3+1;
w:=abs w;
if w>=r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w;
xl:=xl+del >>;
k2:=n-k1-k3;
r7:=sqrt(r7/xn);
if j=1 then princ "log vs taylor expansion"
else if j=2 then princ "log vs 17/16 on"
else if j=3 then princ "log10 vs 11/10 on"
else princ "log vs squared arg";
terpri();
if not (j=1) then heading("log",n,a,b)
else << princ list("c=",c); terpri() >>;
errsigns(k1,k2,k3);
w:=-999.0;
if r6 neq zero then w:=log(abs r6)/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log(abs r7)/albeta;
rmserr(w,max(ait+w,zero));
if j=1 then << a:=sqrt(half); b:=15.0/16.0 >>
else if j=2 then << a:=sqrt(tenth); b:=0.9 >>
else << a:=16.0; b:=240.0 >> >>;
for i:=1:5 do <<
x:=rnd();
x:=x+x+15.0;
y:=one/x;
z:=log x+log y;
specials("log(1/x)+log(x) ",x,z) >>;
specials( "log(1) ",one,log(one));
specials( "xmin ",xmin,log(xmin));
specials( "xmax ",xmax,log(xmax));
errortest list('log, -2.0);
errortest list('log, zero);
princ "end of tests on log"; terpri()
end;
symbolic procedure test!-exp();
begin
scalar a,ait,albeta,b,beta,d,del,one,r6,r7,two,ten,v,w,
x,xl,xn,x1,y,z,zero,zz,i1,j,k1,k2,k3,n;
beta:=float ibeta;
albeta:=log beta;
ait:=float it;
one:=float 1;
two:=float 2;
ten:=float 10;
zero:=float 0;
v:=0.0625;
a:=two;
b:=log(a)*0.5;
a:=-b+v;
d:=log(0.9*xmax);
n:=2000;
xn:=float n;
i1:=0;
for j:=1:3 do <<
k1:=0;
k3:=0;
x1:=zero;
r6:=zero;
r7:=zero;
del:=(b-a)/xn;
xl:=a;
for i:=1:n do <<
x:=del*rnd()+xl;
y:=x-v;
if y<zero then x:=y+v;
z:=exp x;
zz:=exp y;
if j neq 1 then
z:=z*0.0625 - z*2.4453321046920570389/1000.0
else z:=z-z*6.058693718652421388/100.0;
w:=one;
if zz neq zero then w:=(z-zz)/zz;
if w < zero then k1:=k1+1;
if w>zero then k3:=k3+1;
w:=abs w;
if w>r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w;
xl:=xl+del >>;
k2:=n-k1-k3;
r7:=sqrt(r7/xn);
princ list("exp with v=",v); terpri();
heading("exp",n,a,b);
errsigns(k1,k2,k3);
w:=-999.0;
if r6 neq zero then w:=log abs r6/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log abs r7/albeta;
rmserr(w,max(ait+w,zero));
if j neq 2 then <<
v:=45.0 / 16.0;
a:=-10*b;
b:=4.0 * xmin * (beta**it);
b:=log b >>
else << a:=-two * a;
b:=ten * a;
if b<d then b:=d >> >>;
for i:=1:5 do <<
x:=rnd()*beta;
y:=-x;
z:=exp x * exp y - one;
specials("exp(x)*exp(-x)-1 ",x,z) >>;
specials("exp(0)-1 ",zero,exp(zero)-one);
x:=float fix log xmin;
specials("xmin ",x,exp(x));
x:=float fix log xmax;
specials("xmax ",x,exp x);
x:=x / two;
v:=x/two;
y:=exp x;
z:=exp v;
z:=z*z;
princ list(x,y,v,z); terpri();
x:=-one/sqrt(xmin);
errortest list('exp, x);
x:=-x;
errortest list('exp, x);
princ "end of test on exp"; terpri()
end;
symbolic procedure test!-sin();
begin
scalar a,ait,albeta,b,beta,betap,c,del,one,r6,r7,three,w,
x,xl,xn,x1,y,z,zero,zz,expon,i1,j,k1,k2,k3,n;
beta:=float ibeta;
albeta:=log beta;
ait:=float it;
one:=float 1;
three:=float 3;
zero:=float 0;
a:=zero;
b:=1.570796327;
c:=b;
n:=2000;
xn:=float n;
i1:=0;
for j:=1:3 do <<
princ list("start of test ",j," of sin/cos"); terpri();
k1:=0;
k3:=0;
x1:=zero;
r6:=zero;
r7:=zero;
del:=(b-a)/xn;
xl:=a;
for i:=1:n do <<
x:=del*rnd()+xl;
y:=x/three;
y:=(x+y)-x;
x:=three*y;
if j neq 3 then <<
z:=sin x;
zz:=sin y;
w:=one;
if z neq zero then w:=(z-zz*(three-4.0*zz*zz))/z >>
else <<
z:=cos x;
zz:=cos y;
w:=one;
if z neq zero then w:=(z+zz*(three-4.0*zz*zz))/z >>;
if w < zero then k1:=k1+1;
if w>zero then k3:=k3+1;
w:=abs w;
if w>r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w;
xl:=xl+del >>;
k2:=n-k1-k3;
r7:=sqrt(r7/xn);
heading("sin/cos",n,a,b);
errsigns(k1,k2,k3);
w:=-999.0;
if r6 neq zero then w:=log abs r6/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log abs r7/albeta;
rmserr(w,max(ait+w,zero));
a:=18.84955592;
if j=2 then a:=b+c;
b:=a+c >>;
c:=one/beta**(it/2);
z:=(sin(a+c)-sin(a-c))/(c+c);
princ list("this should be about 1.0 ",z); terpri();
for i:=1:5 do <<
x:=rnd()*a;
z:=sin x + sin(-x);
specials("sin(x)+sin(-x) ",x,z) >>;
betap:=beta**it;
x:=rnd()/betap;
for i:=1:5 do <<
z:=x-sin x;
specials("small args to sin ",x,z);
x := x/beta >>;
for i:=1:5 do <<
x:=rnd()*a;
z:=cos(x)-cos(-x);
specials("cos(x)-cos(-x) ",x,z) >>;
expon:=float minexp * 0.75;
x:=beta**expon;
y:=sin x;
specials("tiny arg ",x,y);
%-- z:=sqrt betap;
%-- x:=z*(one-epsneg);
%-- y:=sin x;
%-- specials("?? ",x,y);
%-- y:=sin z;
%-- specials("??+ ",z,y);
%-- x:=z*(one+eps);
%-- y:=sin x;
%-- specials("??++",x,y);
x:=betap;
errortest list('sin, x);
princ "end of test on sin/cos"; terpri()
end;
symbolic procedure test!-tan();
begin
scalar a,ait,albeta,b,beta,betap,c1,c2,del,half,pi,r6,r7,
w,x,xl,xn,x1,y,z,zero,zz,i,i1,j,k1,k2,k3,n,expon;
beta:=float ibeta;
albeta:=log beta;
ait:=float it;
zero:=float 0;
half:=0.5;
pi:=3.1415926;
a:=zero;
b:=pi*0.25;
n:=2000;
xn:=float n;
i1:=0;
for j:=1:4 do <<
princ list("start of test ",j," of tan/cot"); terpri();
k1:=0;
k3:=0;
x1:=zero;
r6:=zero;
r7:=zero;
del:=(b-a)/xn;
xl:=a;
for i:=1:n do <<
x:=del*rnd()+xl;
y:=x*half;
x:=y+y;
if not (j=4) then <<
z:=tan x;
zz:=tan y;
w:=1.0;
if z neq zero then <<
w:=((half-zz)+half)*((half+zz)+half);
w:=(z-(zz+zz)/w)/z >> >>
else <<
z:=cot x;
zz:=cot y;
w:=1.0;
if z neq zero then <<
w:=((half-zz)+half)*((half+zz)+half);
w:=(z+(w/(zz+zz)))/z >> >>;
if w>zero then k1:=k1+1;
if w<zero then k3:=k3+1;
w:=abs w;
if w>r6 then << r6:=w; x1:=x >>;
r7:=r7+w*w;
xl:=xl+del >>;
k2:=n-k1-k3;
r7:=sqrt(r7/xn);
heading("tan/cot",n,a,b);
errsigns(k1,k2,k3);
w:=-999.0;
if r6 neq zero then w:=log abs r6/albeta;
maxrel(w,x1,max(ait+w,zero));
w:=-999.0;
if r7 neq zero then w:=log abs r7/albeta;
rmserr(w,max(ait+w,zero));
if j=1 then <<
a:=pi*0.875;
b:=pi*1.125 >>
else <<
a:=6.0*pi;
b:=a+pi*0.25 >> >>;
for i:=1:5 do <<
x:=rnd()*a;
z:=tan x + tan(-x);
specials("tan(x)+tan(-x) ",x,z) >>;
betap:=beta**it;
x:=rnd()/betap;
for i:=1:5 do <<
z:=x-tan x;
specials("small args to tan ",x,z);
x := x/beta >>;
expon:=float minexp * 0.75;
x:=beta**expon;
y:=tan x;
specials("tiny arg ",x,y);
c1:=-225.0;
c2:=-0.950846454195142026;
x:=11.0;
y:=tan x;
w:=((c1-y)+c2)/(c1+c2);
z:=log abs w/albeta;
princ list("the relative error in tan 11 is ",z,max(ait+z,zero));
terpri();
x:=beta**(it/2);
errortest list('tan, x);
x:=betap;
errortest list('tan, x);
princ "end of test on tan/cot"; terpri()
end;
on echo, time;
out "fptest.log";
fptest();
out t;
end;