module arith; % Header module for real arith package.
create!-package('(arith smlbflot bfauxil paraset math rounded comprd
rdelem crelem bfelem), nil);
endmodule;
module smlbflot; % Basic support for bigfloat arithmetic.
% Authors: S.L. Kameny and T. Sasaki.
% Mod 1.92, 3 July 1990.
fluid '(!*bfspace !*nat !:prec!:);
global '(!!nfpd);
switch bfspace; symbolic(!*bfspace := t);
symbolic smacro procedure mt!: u;
% This function selects the mantissa of U, a bigfloat representation
% of a number.
cadr u;
symbolic smacro procedure ep!: u;
% This function selects the exponent of U, a bigfloat representation
% of a number.
cddr u;
symbolic procedure bflerrmsg u;
% Revised error message for BFLOAT module, using error, not rederr.
error(0,list("Invalid argument to",u));
symbolic procedure i2bf!: u; '!:bf!: . u . 0;
symbolic procedure !*rn2bf u;
apply2(get('!:bf!:,'quotient),apply1(x,cadr u),apply1(x,cddr u))
where x=get('!:bf!:,'i2d);
symbolic procedure !*ft2bf u; conv!:a2bf cdr u;
symbolic procedure bfzerop!: u;
% This is possibly too restricted a definition.
mt!: u = 0;
symbolic smacro procedure preci!: nmbr;
% This function counts the precision of a number "n". NMBR is a
% bigfloat representation of "n".
length explode abs mt!: nmbr;
symbolic procedure bfprin!: u;
if preci!: u>!!nfpd then bfprin0 u
else (bfprin0 u where !*bfspace=nil);
symbolic procedure bfprin0 nmbr;
% Prints a big-float in a variety of formats. Still needs work
% for Fortran output.
begin integer j,k; scalar u,v;
if bfzerop!:(nmbr) then return prin2!* '!0!.0;
u := explode abs(j := mt!: nmbr);
k := ep!: nmbr;
if k>=0 then if k>5 then go to etype
else <<v := list('!.,'!0);
while (k := k-1)>=0 do v := '!0 . v;
u := nconc(u,v)>>
else if (k := order!:(nmbr)+1)>0
then <<v := u;
while (k := k-1)>0 do v := cdr v;
rplacd(v,'!. . cdr v)>>
else if k<-10 then go to etype
else <<while (k := k+1)<=0 do u := '!0 . u;
u := '!0 . '!. . u>>;
bfprin1(u,j);
return nmbr;
etype:
if null( cdr(u)) then rplacd(u , list('!0));
u:= car u . '!. . cdr u;
j := bfprin1(u,j);
if not(!*nat and !*bfspace) then prin2!* "E" else
if j=0 then <<prin2!*("E " ); j:=2>> else
if j=1 then <<prin2!*(" E " ); j:=4>> else
if j=2 then <<prin2!*(" E "); j:=0>> else
if j=3 then <<prin2!*(" E " ); j:=0>> else
if j=4 then <<prin2!*(" E "); j:=2>>;
u:=explode( k:=order!:(nmbr));
if k>=0 then u:=cons('!+,u);
while u do <<prin2!*( car(u)); u:=cdr(u); j:=j+1;
if j=5 and !*nat and !*bfspace then <<prin2!*(" "); j:=0>>>>;
return nmbr
end;
symbolic procedure bfprin1(u,j);
begin
if j<0 then u := '!- . u;
j := 0;
for each char in u do <<prin2!* char; j := j+1;
if j=5 then <<if !*nat and !*bfspace then prin2!* '! ;
j := 0>>>>;
return j
end;
symbolic smacro procedure make!:bf(mt,ep); '!:bf!: . (mt . ep);
symbolic procedure bfp!:(x);
% This function returns true if X is a bigfloat representation,
% else NIL. X is any Lisp entity.
eqcar(x,'!:bf!:);
symbolic procedure order!: nmbr;
% This function counts the order of a number "n". NMBR is a bigfloat
% representation of "n".
% **** ORDER(n)=k if 10**k <= ABS(n) < 10**(k+1)
% **** when n is not 0, and ORDER(0)=0.
if mt!: nmbr = 0 then 0 else preci!: nmbr + ep!: nmbr - 1;
symbolic procedure conv!:a2bf(n);
if bfp!: n then n
else if fixp n then make!:bf(n, 0)
else if floatp n then read!:num n
else if stringp n then read!:num n
else if ratnump!: n
then conv!:r2bf(n, if !:prec!: then !:prec!: else 50)
else if not atom n and idp car n and get(car n,'dname)
then apply1(get(car n,'!:bf!:),n)
else bflerrmsg 'conv!:a2bf;
symbolic procedure decprec!:(nmbr, k);
make!:bf( mt!: nmbr / 10**k, ep!: nmbr + k);
symbolic procedure incprec!:(nmbr, k);
make!:bf( mt!: nmbr * 10**k, ep!: nmbr - k);
symbolic procedure conv!:mt(nmbr, k);
% This function converts a number "n" to an equivalent number of
% precision K by rounding "n" or adding "0"s to "n". NMBR is a
% bigfloat representation of "n". K is a positive integer.
if bfp!: nmbr and fixp k and k > 0
then if (k := preci!: nmbr - k) = 0 then nmbr
else if k < 0 then incprec!:(nmbr, -k)
else round!:last(decprec!:(nmbr, k - 1))
else bflerrmsg 'conv!:mt;
symbolic procedure round!:mt(nmbr, k);
% This function rounds a number "n" at the (K+1)th place and returns
% an equivalent number of precision K if the precision of "n" is
% greater than K, else it returns the given number unchanged. NMBR
% is a bigfloat representation of "n". K is a positive integer.
if bfp!: nmbr and fixp k and k > 0 then
if (k := preci!: nmbr - k - 1) < 0 then nmbr
else if k = 0 then round!:last nmbr
else round!:last decprec!:(nmbr, k)
else bflerrmsg 'round!:mt;
symbolic procedure round!:last nmbr;
% This function rounds a number "n" at its last place. NMBR is a
% bigfloat representation of "n".
begin scalar n;
n := divide(abs mt!: nmbr, 10);
if cdr n < 5 then n := car n else n := car n + 1;
if mt!: nmbr < 0 then n := -n;
return make!:bf(n, ep!: nmbr + 1)
end;
symbolic procedure read!:num(n);
% This function reads a number or a number-like entity N
% and constructs a bigfloat representation of it.
% N is an integer, a floating-point number, or a string
% representing a number.
% **** If the system does not accept or may incorrectly
% **** accept the floating-point numbers, you can
% **** input them as strings such as "1.234E-56",
% **** "-78.90 D+12" , "+3456 B -78", or "901/234".
% **** A rational number in a string form is converted
% **** to a bigfloat of precision !:PREC!: if
% **** !:PREC!: is not NIL, else the precision of
% **** the result is set 50.
% **** Some systems set the maximum size of strings. If
% **** you want to input long numbers exceeding
% **** such a maximum size, please use READ!:LNUM.
if fixp n then make!:bf(n, 0)
else if not(numberp n or stringp n) then bflerrmsg 'read!:num
else begin integer j,m,sign; scalar ch,u,v,l,appear!.,appear!/;
j := m := 0;
sign := 1;
u := v := appear!. := appear!/ := nil;
l := explode n;
loop: ch := car l;
if digit ch then << u := ch . u; j := j + 1 >>
else if ch eq '!. then << appear!. := t; j := 0 >>
else if ch eq '!/ then << appear!/ := t; v := u; u := nil >>
else if ch eq '!- then sign := -1
else if ch memq '(!E !D !B !e !d !b) then go to jump; %JBM
if l := cdr l then goto loop else goto make;
jump: while l := cdr l do
<<if digit(ch := car l) or ch eq '!-
then v := ch . v >>;
l := reverse v;
if car l eq '!- then m := - compress cdr l
else m:= compress l;
make: u := reverse u;
v := reverse v;
if appear!/ then
return conv!:r2bf(make!:ratnum(sign*compress v,compress u),
if !:prec!: then !:prec!: else 50);
if appear!. then j := - j else j := 0;
if sign = 1 then u := compress u else u := - compress u;
return make!:bf(u, j + m)
end;
symbolic procedure abs!: nmbr;
% This function makes the absolute value of "n". N is a bigfloat
% representation of "n".
if mt!: nmbr > 0 then nmbr else make!:bf(- mt!: nmbr, ep!: nmbr);
symbolic procedure minus!: nmbr;
% This function makes the minus number of "n". N is a bigfloat
% representation of "n".
make!:bf(- mt!: nmbr, ep!: nmbr);
symbolic procedure plus!:(n1,n2);
begin scalar m1,m2,e1,e2,d; return
if (m1 := mt!: n1)=0 then n2
else if (m2 := mt!: n2)=0 then n1
else if (d := (e1 := ep!: n1)-(e2 := ep!: n2))=0
then make!:bf(m1+m2, e1)
else if d>0 then make!:bf(m1*10**d+m2,e2)
else make!:bf(m1+m2*10**(-d),e1) end;
symbolic procedure difference!:(n1,n2);
begin scalar m1,m2,e1,e2,d; return
if (m1 := mt!: n1)=0 then minus!: n2
else if (m2 := mt!: n2)=0 then n1
else if (d := (e1 := ep!: n1)-(e2 := ep!: n2))=0
then make!:bf(m1-m2, e1)
else if d>0 then make!:bf(m1*10**d-m2,e2)
else make!:bf(m1-m2*10**(-d),e1) end;
symbolic procedure times!:(n1, n2);
% This function calculates the product of "n1" and "n2".
% N1 and N2 are bigfloat representations of "n1" and "n2".
make!:bf(mt!: n1 * mt!: n2, ep!: n1 + ep!: n2);
symbolic procedure divide!:(n1,n2,k);
% This function calculates the quotient of "n1" and "n2", with the
% precision K, by rounding the ratio of "n1" and "n2" at the (K+1)th
% place. N1 and N2 are bigfloat representations of "n1" and "n2".
% K is any positive integer.
begin
n1 := conv!:mt(n1, k + preci!: n2 + 1);
n1 := make!:bf(mt!: n1 / mt!: n2, ep!: n1 - ep!: n2);
return round!:mt(n1, k)
end;
symbolic procedure min!:(n1,n2);
% This function returns the smaller of "n1" and "n2".
% N1 and N2 are bigfloat representations of "n1" and "n2".
if lessp!:(n2, n1) then n2 else n1;
symbolic procedure greaterp!:(a,b);
% This function calculates the a > b, but avoids
% generating large numbers if magnitude difference is large.
if ep!: a=ep!: b then mt!: a>mt!: b else
(((if d=0 then ma>mb else
((if d>p2 then ma>0 else if d<-p2 then mb<0
else if d>0 then ma*10**d>mb
else ma>mb*10**-d)
where p2=2*!:prec!:))
where d=ep!: a - ep!: b, ma=mt!: a, mb=mt!: b)
where a= normbf a, b=normbf b);
symbolic procedure equal!:(a,b);
%tests bfloats for a=b rapidly without generating digits. %SK
zerop mt!: a and zerop mt!: b or
ep!:(a := normbf a)=ep!:(b := normbf b) and mt!: a=mt!: b;
symbolic procedure lessp!:(n1, n2);
% This function returns T if "n1" < "n2" else returns NIL.
% N1 and N2 are bigfloat representations of "n1" and "n2".
greaterp!:(n2, n1);
symbolic procedure leq!:(n1, n2);
% This function returns T if "n1" <= "n2" else returns NIL.
% N1 and N2 are bigfloat representations of "n1" and "n2".
not greaterp!:(n1, n2);
symbolic procedure minusp!: x;
% This function returns T if "x"<0 else returns NIL.
% X is any Lisp entity.
bfp!: x and mt!: x < 0;
symbolic procedure make!:ratnum(nm,dn);
% This function constructs an internal representation
% of a rational number composed of the numerator
% NM and the denominator DN.
% NM and DN are any integers (positive or negative).
% **** Four routines in this section are temporary.
% **** That is, if your system has own routines
% **** for rational number arithmetic, you can
% **** accommodate our system to yours only by
% **** redefining these four routines.
if zerop dn then rerror(arith,1,"Zero divisor in make:ratnum")
else if dn > 0 then '!:ratnum!: . (nm . dn)
else '!:ratnum!: . (-nm . -dn);
symbolic procedure ratnump!:(x);
% This function returns T if X is a rational number
% representation, else NIL.
% X is any Lisp entity.
eqcar(x,'!:ratnum!:); %JBM Change to EQCAR.
symbolic smacro procedure numr!: rnmbr;
% This function selects the numerator of a rational number "n".
% RNMBR is a rational number representation of "n".
cadr rnmbr;
symbolic smacro procedure denm!: rnmbr;
% This function selects the denominator of a rational number "n".
% RNMBR is a rational number representation of "n".
cddr rnmbr;
symbolic procedure conv!:r2bf(rnmbr,k);
% This function converts a rational number RNMBR to a bigfloat of
% precision K, i.e., a bigfloat representation with a given
% precision. RNMBR is a rational number representation. K is a
% positive integer.
if ratnump!: rnmbr and fixp k and k > 0
then divide!:(make!:bf( numr!: rnmbr, 0),
make!:bf( denm!: rnmbr, 0),k)
else bflerrmsg 'conv!:r2bf;
endmodule;
module bfauxil; % Support for the roots package and ROUNDED domain.
% Author: Stanley L. Kameny <valley!stan@rand.org>.
% Version and Date: Mod 1.92, 20 Mar 1990.
% Copyright (c) 1988,1989,1990. Stanley L. Kameny. All Rights Reserved.
Comment support for modules allroot and isoroot, and for ROUNDED
domain logic;
fluid '(!:prec!:);
global '(bfone!* bfhalf!* bfz!*);
global '(!!nfpd !!decflt !!shbinfl vv!! !!flint);
global '(!!minflbf !!maxflbf);
bfz!* := '!:bf!: . (0 . 0);
bfone!* := '!:bf!: . (1 . 0); bfhalf!* := '!:bf!: . (5 . -1);
smacro procedure mt!: u; cadr u;
smacro procedure ep!: u; cddr u;
smacro procedure make!:bf(mt,ep); '!:bf!: . (mt . ep);
symbolic smacro procedure rndpwr j; normbf round!:mt(j,!:prec!:);
symbolic procedure retag u; if atom u then u else '!:bf!: . u;
symbolic procedure normbf x;
begin scalar mt,s;integer ep;
if (mt := mt!: x)=0 then go to ret;
if mt<0 then <<mt := -mt; s := t>>;
ep := ep!: x;
mt := mt*10;
while remainder(mt, 100000000)=0 do <<
mt := mt/100000000;
ep := ep+8 >>;
while remainder(mt, 10)=0 do <<
mt := mt/10;
ep := ep+1 >>;
ep := ep-1;
if s then mt := -mt;
ret: return '!:bf!: . (mt . ep) end;
symbolic procedure divbf(u,v); normbf divide!:(u,v,!:prec!:);
symbolic procedure bfdivide(u,v);
if atom u then u/v else divbf(u,v);
symbolic procedure timbf(u,v); rndpwr times!:(u,v);
symbolic procedure bftimes(u,v); if atom u then u*v else timbf(u,v);
symbolic procedure plubf(a,b);
% This function calculates the normalized rounded sum of a and b,
% but avoids generating large numbers if magnitude difference is large.
rndpwr
(if ma=0 then b else if mb=0 then a else
((if d=0 then make!:bf(ma+mb,ea) else
((if d>p2 then a else if d<-p2 then b
else if d>0 then make!:bf(ma*10**d+mb,eb)
else make!:bf(ma+mb*10**-d,ea))
where p2=2*!:prec!:))
where d=ea-eb where ea=ep!: a,eb=ep!: b))
where ma=mt!: a,mb=mt!: b;
symbolic procedure bfplus(u,v); if atom u then u+v else plubf(u,v);
symbolic procedure difbf(a,b);
% This function calculates the normalized rounded difference of a and b,
% but avoids generating large numbers if magnitude difference is large.
rndpwr
(if ma=0 then minus!: b else if mb=0 then a else
((if d=0 then make!:bf(ma-mb,ea) else
((if d>p2 then a else if d<-p2 then minus!: b
else if d>0 then make!:bf(ma*10**d-mb,eb)
else make!:bf(ma-mb*10**-d,ea))
where p2=2*!:prec!:))
where d=ea-eb where ea=ep!: a,eb=ep!: b))
where ma=mt!: a,mb=mt!: b;
symbolic procedure bfdiffer(u,v); if atom u then u-v else difbf(u,v);
symbolic procedure invbf u; divbf(bfone!*,u);
symbolic procedure bfinverse u; if atom u then 1.0/u else invbf u;
symbolic procedure bfminus u; if atom u then -u else minus!: u;
symbolic procedure bflessp(a,b); if atom a then a<b else grpbf(b,a);
symbolic procedure grpbf(a,b);
% This function calculates a > b, but avoids generating large numbers
% if magnitude difference is large.
<<a := mt!: a; b := mt!: b;
if d=0 then a>b else
((if d>p2 then a>0 else if d<-p2 then b<0
else if d>0 then a*10**d>b
else a>b*10**-d)
where p2=2*!:prec!:)>>
where d=ep!: a -ep!: b;
symbolic procedure bfminusp u; if atom u then minusp u else minusp!: u;
symbolic procedure bfzp u; if atom u then zerop u else mt!: u=0;
symbolic procedure bfnzp u; not bfzp u;
symbolic procedure bfloat x; if floatp x then fl2bf x else
normbf(if atom x then read!:num x else x);
symbolic procedure !!shbinflp;
begin integer n; vv!! := 9.0;
while n<300 and not errorp errorset!*('(vv!!!*1e10),nil)
do n := n+10;
return n<300 end;
symbolic procedure vv!!!*1e10; vv!! := vv!!*1.0e10;
symbolic (!!shbinfl := !!shbinflp());
symbolic procedure bfsqrt x;
% computes sqrt x by Newton's method.
if minusp!: x then terrlst(x,'bfsqrt) else
begin scalar nx,dx,dc,k2,nf;
if bfzerop!: x then return bfz!*;
k2 := !:prec!:+2;
dc := '!:bf!: . (1 . (-k2+(order!: x +3)/2));
nx := '!:bf!: . if remainder(ep!:(nx := conv!:mt(x,2)),2)=0
then ((3+2*mt!: nx/25) . (ep!: nx/2))
else ((10+2*mt!: nx/9) . ((ep!: nx-1)/2));
nf := 1;
loop: if (nf := 2*nf)>k2 then nf := k2;
dx := times!:(bfhalf!*,plus!:(divide!:(x,nx,nf),nx));
if nf>=k2 and not greaterp!:(abs!: difference!:(dx,nx),dc)
then return rndpwr nx;
nx := dx; go to loop end;
symbolic procedure realrat x;
begin scalar d,g;
if bfp!: x then go to bf;
if eqcar(x,'quotient) then
if fixp cadr x and fixp caddr x then
<<x := if (d := caddr x)<0 then -cadr x else cadr x;
d := abs d; go to ret>>
else x := cadr x/caddr x;
if zerop x then return (0 . 1);
if not floatp x then return (x . 1);
x := bfloat x;
bf: d := cddr(x := normbf x); x := cadr x;
if x=0 then return (0 . 1);
if d<0 then d := 10**(-d) else <<x := x*10**d; d := 1>>;
ret: g := gcdn(abs x,d); return (x/g) . (d/g) end;
remflag ('(fl2int),'lose);
symbolic procedure fl2int x;
<<x := fl2bf x;
(if d=0 then m else if d>0 then m*10**d else m/(10**-d))
where m=mt!: x,d=ep!: x>>;
flag ('(fl2int),'lose);
symbolic procedure fl2bf x;
(if zerop x then z else
begin scalar s,r; integer d,d1;
if x<0 then <<x := -x; s := t>>;
% convert x to an integer equivalent;
r := normbf read!:num x;
d1 := (d := ep!: r+length explode mt!: r-!:prec!:-1)/2;
x := x*10.0**-d1; x := x*10.0**(d1-d);
% int:
repeat z := plubf(z,r := normbf read!:num x)
until
mt!: cut!:ep(read!:num(abs(x := x-mt!: r*10.0**ep!: r)+0.5),0)=0;
x := mt!:(z := normbf round!:mt(z,!:prec!:-1));
return '!:bf!: . (if s then -x else x) . (d+ep!: z) end)
where !:prec!:=!!nfpd,z=bfz!*;
symbolic procedure cflot x; if floatp x then x else bf2flr bfloat x;
global '(!1!.e4); !1!.e4 := 10000.0;
symbolic procedure bf2flr u;
% Converts bfloat to float by rounding at !!nfpd digits.
% We use error1 rather than rerror, because we want to catch such an
% error in an errorset.
begin scalar x,ep,y;
x := 0.0;
u := bfloat u;
if bfzerop!: u then return x;
ep := ep!:(u := round!:mt(u,!!nfpd));
if grpbf(!!minflbf,y := abs!: u) or grpbf(y,!!maxflbf)
then error1()
else ep := 10.0**ep;
y := u := mt!: u;
lp: x := x + (y - 10000*(u := u/10000))*ep;
if u=0 then return x;
y := u; ep := ep * !1!.e4;
go to lp
end;
symbolic procedure gf2fl a; % force into float format.
if atom a then a else if eqcar(a,'!:bf!:) then bf2flr a
else (gf2fl car a) . gf2fl cdr a;
symbolic procedure gf2bf a; if a then % force into bfloat format.
if atom a then bfloat a else if eqcar(a,'!:bf!:) then a
else (gf2bf car a) . gf2bf cdr a;
symbolic procedure r2bf u;
% translate any real number object to bigfloat.
if atom u then bfloat u
else if eqcar(u,'!:bf!:) then u
else if numberp car u then divbf(i2bf!: car u,i2bf!: cdr u)
else if eqcar(u,'quotient) then
bfdivide(i2bf!: cadr u,i2bf!: caddr u)
else if eqcar(u,'!:rn!:) then r2bf cdr u
else r2bf cadr u;
symbolic procedure r2fl u;
% translate any real number object to float.
if u=0 then 0.0
else if atom u then cflot u
else if numberp car u then (cflot car u)/cdr u
else if eqcar(u,'quotient) then (cflot cadr u)/caddr u
else if eqcar(u,'!:bf!:) then bf2flr u
else if eqcar(u,'!:rn!:) then r2fl cdr u
else r2fl cadr u;
symbolic procedure gfplus(u,v);
if atom car u then gffplus(u,v) else gbfplus(u,v);
symbolic procedure gffplus(u,v); (car u+car v) . (cdr u+cdr v);
symbolic procedure gbfplus(u,v);
(plubf(car u,car v)) . plubf(cdr u,cdr v);
symbolic procedure gfdiffer(u,v);
if atom car u then gffdiff(u,v) else gbfdiff(u,v);
symbolic procedure gffdiff(u,v); (car u-car v) . (cdr u-cdr v);
symbolic procedure gbfdiff(u,v);
(difbf(car u,car v)) . difbf(cdr u,cdr v);
symbolic procedure gftimes(u,v);
if atom car u then gfftimes(u,v) else gbftimes(u,v);
symbolic procedure gfftimes(u,v);
begin scalar ru,iu,rv,iv;
ru := car u; iu := cdr u; rv := car v; iv := cdr v;
return (ru*rv-iu*iv) . (ru*iv+iu*rv) end;
symbolic procedure gbftimes(u,v);
begin scalar ru,iu,rv,iv;
ru := car u; iu := cdr u; rv := car v; iv := cdr v;
return (difbf(timbf(ru,rv),timbf(iu,iv))) .
plubf(timbf(ru,iv),timbf(iu,rv)) end;
symbolic procedure gfquotient(u,v);
if atom car u then gffquot(u,v) else gbfquot(u,v);
symbolic procedure gffquot(u,v);
begin scalar ru,iu,rv,iv,d;
ru := car u; iu := cdr u; rv := car v; iv := cdr v;
d := rv*rv+iv*iv;
return ((ru*rv+iu*iv)/d) . ((iu*rv-ru*iv)/d) end;
symbolic procedure gbfquot(u,v);
begin scalar ru,iu,rv,iv,d;
ru := car u; iu := cdr u; rv := car v; iv := cdr v;
d := plubf(timbf(rv,rv),timbf(iv,iv));
return divbf(plubf(timbf(ru,rv),timbf(iu,iv)),d) .
divbf(difbf(timbf(iu,rv),timbf(ru,iv)),d) end;
symbolic procedure gfminus u; (bfminus car u) . (bfminus cdr u);
symbolic procedure gfrotate u; (bfminus cdr u) . (car u);
symbolic procedure gfrl u; car u;
symbolic procedure gfim u; cdr u;
symbolic procedure gfzerop u;
if not atom car u then cadar u=0 and caddr u=0
else equal(u,(0.0 . 0.0));
symbolic procedure gfdot(u,v);
if atom car u then gffdot(u,v) else gbfdot(u,v);
symbolic procedure gffdot(u,v); car u*car v+cdr u*cdr v;
symbolic procedure gbfdot(u,v);
plubf(timbf(car u,car v),timbf(cdr u,cdr v));
symbolic procedure gfrsq u; gfdot(u,u);
symbolic procedure gffrsq u; car u*car u+cdr u*cdr u;
symbolic procedure gbfrsq u;
plubf(timbf(car u,car u),timbf(cdr u,cdr u));
symbolic procedure gffmult(r,u); (r*car u) . (r*cdr u);
symbolic procedure gffsqrt x;
begin scalar x0,nx,xd,xd0,rl,im; rl := gfrl x; im := gfim x;
rl := sqrt(hypot(rl,im)/2+rl/2); im := im/(2*rl); nx := rl . im;
repeat
<<x0 := nx;
nx := gffmult(0.5,gffplus(x0,gffquot(x,x0)));
xd0 := xd; xd := gffrsq gffdiff(x,gfftimes(nx,nx))>>
until xd0 and xd0-xd<=0.0; return x0 end;
symbolic procedure gbfmult(r,u);
<<r := bfloat r; (timbf(r,car u)) . (timbf(r,cdr u))>>;
symbolic procedure gbfsqrt x;
begin scalar x0,nx,xd,xd0,rl;
nx :=
<<rl := (bfsqrt timbf(bfhalf!*,plubf(bfsqrt gfrsq x,gfrl x)));
rl . timbf(bfhalf!*,divbf(gfim x,rl))>>;
repeat
<<x0 := nx;
nx := gbfmult(bfhalf!*,gbfplus(x0,gbfquot(x,x0)));
xd0 := xd; xd := gbfrsq gbfdiff(x,gbftimes(nx,nx))>>
until xd0 and mt!: difbf(xd0,xd)<=0; return x0 end;
symbolic procedure gfsqrt x;
% computes gfsqrt x by Newton's method, for both gf and gbf.
begin scalar xn,neg;
if gfzerop x then return x;
if bfminusp gfrl x then <<x := gfminus x; neg := t>>;
if bfzp gfim x then
<<x := gfrl x;
xn := rl2gfc(if atom x then sqrt x else bfsqrt x);
go to ret>>;
xn := if atom car x then gffsqrt x else gbfsqrt x;
ret: return if neg then gfrotate xn else xn end;
symbolic procedure rl2gfc x; x . if atom x then 0.0 else bfz!*;
symbolic procedure sgn x; if x>0 then 1 else if x<0 then -1 else 0;
symbolic procedure cut!:mt(nmbr,k);
% This function returns a given number "n" unchanged
% if its precision is not greater than K, else it
% cuts off its mantissa at the (K+1)th place and
% returns an equivalent number of precision K.
% **** CAUTION! No rounding is made.
% NMBR is a BIG-FLOAT representation of "n".
% K is a positive integer.
if bfp!: nmbr and fixp k and k > 0 then
if (k := length explode abs mt!: nmbr - k) <= 0 then nmbr
else decprec!:(nmbr, k)
else bflerrmsg 'cut!:mt$
symbolic procedure exptbf(x,n,a);
% Computes a*x**n in bfloat arithmetic for positive x and integer n.
begin scalar m;
lp: m := n; if m-2*(n := n/2) neq 0 then a := timbf(a,x);
if n=0 then return a;
x := timbf(x,x); go to lp end;
symbolic procedure icbrt x;
% computes icbrt by Newton's method, avoiding magnitude problems.
if abs(x := fix2 x) < 2 then x else
begin scalar s,r,g,trv,nx;
if x<0 then x := -x else s := t;
if x<100 then <<nx := 1; go to loop>>;
r := remainder(ep!:(nx := conv!:mt('!:bf!: . (x . 0),3)),3);
nx := if r=0 then (5+mt!: nx/179)*10**(ep!: nx/3)
else if r=1 or r=-2 then (10+mt!: nx/74)*10**((ep!: nx-1)/3)
else (22+mt!: nx/35)*10**((ep!: nx-2)/3);
loop: trv := nx; nx := 2*trv+x/(trv*trv); nx := nx/3;
if g and nx>=trv then go to ret;
g := t; go to loop;
ret: return if s then trv else -trv end;
symbolic procedure isqrt x;
% computes isqrt x by Newton's method, avoiding magnitude problems.
if x<0 then terrlst(x,'isqrt)
else if (x := fix2 x)=0 or x=1 then x else
begin scalar trv,nx,g;
if x<144 then <<nx := 1; go to loop>>;
nx := if remainder(ep!:(nx := conv!:mt('!:bf!: . (x . 0),2)),2)=0
then (3+2*mt!: nx/25)*10**(ep!: nx/2)
else (10+2*mt!: nx/9)*10**((ep!: nx-1)/2);
loop: trv := nx; nx := (trv + x/trv)/2;
if g and nx>=trv then return trv;
g := t; go to loop end;
symbolic procedure ilog2 x; % logb(10,2) is ~ 485/146.
if x<=0 then terrlst(x,'ilog2) else if (x := fix2 x)<3 then x-1 else
(<<y := 2**(p := p/146); while x>=y do <<y := y*2; p := p+1>>; p-1>>
where p=(length explode x-1)*485,y=1);
symbolic procedure irootn!-round(m,n);
% M,n integral, n>0, answer is nearest integer to m/n.
(m+n/2)/n;
symbolic procedure irootn!-roundup(m,n);
% M,n integral, n>0, answer is least integer >= m/n.
(if cdr quotrem<=0 then car quotrem else 1+car quotrem)
where quotrem = divide(m,n);
symbolic procedure irootn(n,r);
% n,r integral; result 's' approximates rth root of n. That is,
% if n>0 then s**r <= n < (s+1)**r else s**r >= n > (s-1)**r.
% 3363/2378 is an approximation to sqrt 2.
% This definition uses isqrt or icbrt where applicable.
if not fixp n then typerr(n,"integer")
else if not fixp r or r<=0 then typerr(r,"positive integer")
else if n<0
then if evenp r then typerr(r,"odd integer") else -irootn(-n,r)
% Modified here to use isqrt or icbrt for speed.
else if r=1 then n else if r=2 then isqrt n else if r=3 then icbrt n
else begin scalar x,newx,upb;
x:=2**(irootn!-round(1+ilog2 n,r));
newx:=x-irootn!-roundup(x-n/x**(r-1),r);
upb:=(3363*x)/2378;
if upb<newx then newx:=upb;
repeat <<x:=newx; newx:=x-irootn!-roundup(x-n/x**(r-1),r)>>
until newx>=x or newx=0;
return x
end;
put('irootn,'number!-of!-args,2); % For VALUECHK.
endmodule;
module paraset; % Parameter determining module.
% Author: Stanley L. Kameny.
% Paraset.red determines the parameters !!nfpd, !!decflt, and !!maxflexp
% for floating point numbers. !!nfpd, !!decflt and !!maxflexp are
% computed at load, but !!maxflexp may have to be fixed up at run time
% if !!flexperr is true.
symbolic;
global '(!!nfpd !!decflt !!!~xx !!yy !!maxflexp !!flexperr !!plumax);
share !!nfpd,!!decflt,!!maxflexp;
symbolic procedure find!!nfpd;
begin scalar x,y,z;integer i; x:=y:=9.0;
repeat <<x := 10.0*x+y; i := i+1>> until (z := x+1.0)=x;
if 10.0*fix(z/10)-1.0 neq x then i := i-1;
return !!nfpd:=i end;
symbolic procedure find!!decflt;
begin scalar x,y; integer i; x := y := 1.0; find!!nfpd();
repeat <<x := 10.0*x; i := i+1>> until x+y=x;
return !!decflt := i=!!nfpd end;
symbolic procedure find!!maxflexp;
begin scalar z; integer n;
!!!~xx := 9.0;
while not errorp
(z := errorset!*(
'(progn (setq !!yy (plus 9.0 (times !!!~xx 10.0)))
(read!:num !!yy)
(and (not (infinityp !!yy))
(greaterp !!yy !!!~xx))),nil))
and car z do
<<n := n+1; !!!~xx := !!yy>>;
!!flexperr := not errorp z and not car z;
return !!maxflexp := n end;
symbolic procedure infinityp u;
% Check for a representation of an IEEE floating point infinity.
not(x eq '!- or digit x) where x=car explode u;
symbolic procedure !!mfefix;
<<if !!flexperr
then begin integer n;
!!flexperr := nil;
while errorp errorset!*('(explode !!!~xx),nil)
do <<!!!~xx := !!!~xx/10.0; n := n+1>>;
!!maxflexp := !!maxflexp-n
end;
if not !!plumax then roundconstants()>>;
find!!nfpd(); find!!decflt(); find!!maxflexp();
endmodule;
module math; % Mathematical Package for REDUCE.
% Author: Stanley L. Kameny <valley!stan@rand.org>.
% Version and Date: Mod 1.61 25 May 91.
% Copyright (c) 1987, 1988, 1989, 1990, 1991 Stanley L. Kameny.
% All Rights Reserved.
%*******************************************************************
%** **
%** This math package will compute the floating point values of **
%** the usual elementary functions, namely: **
%** sin asin sind asind sinh asinh **
%** cos acos cosd acosd cosh acosh **
%** tan atan tand atand tanh atanh **
%** cot acot cotd acotd coth acoth **
%** sec asec secd asecd sech asech **
%** csc acsc cscd acscd csch acsch **
%** atan2 atan2d **
%** exp ln sqrt **
%** expt log cbrt **
%** logb hypot **
%** log10 floor **
%** ceiling **
%** round **
%** **
%** All functions are computed to the accuracy of the floating- **
%** point precision of the system set up at the time. **
%** **
%*******************************************************************
% Revisions:
% 25 May 91 atan2d added. Function list updated.
% 6 Feb 90 exp, expt, and hyperbolic functions improved.
% 6 Nov 90 find!!nfpd improved; fl2int eliminated (now in bfauxil.)
% 31 Mar 90 fl2int speeded up for very large argument.
% 25 Feb 90 expt modified.
% 15 Oct 89 atan2 and sind,asind family added.
% 8 Oct 89 !!scalsintan,sqrt,expt,and hyperbolics modified.
% 8 Oct 89 hypot,floor,ceiling,round added
% 26 Aug 89 zerop test used in expt
% 20 Jul 89 !!scalsintan revised: same limits for deg and rad
% 17 May 89 find!!nfpd revised (for binary or decimal floats)
% 18 Apr 89 global log10 -> logten (to avoid name conflict)
% 18 Apr 89 !!scalsintan revised (per: Rainer Schoepf)
% 27 Nov 88 log,expt,sqrt revised for speed
% 5 Jun 88 log -> log10; ln -> log; exptfp -> expt (revised)
% 11 Nov 87 hyperbolic fns rewritten: speeded up; improved
% 6 Nov 87 ln,atan rewritten: speeded up. sq!-1 eliminated
% 30 Oct 87 sin,(cos),tan,(cot),exp rewritten; speeded up
%*******************************************************************
%** Basic functions **
%*******************************************************************
symbolic;
global '(math!!label);
math!!label:="Math package mod 1.61, 25 May 91";
symbolic procedure terrlst (x,y);error(0,list(x," invalid for ",y));
global '(!!nfpd !!decflt !!flint !!epsqrt);
if not !!nfpd then flag('(!!nfpd !!decflt),'share);
!!epsqrt := 10.0**((-1-!!nfpd)/2);
!!flint := 10.0**!!nfpd;
symbolic procedure sqrt x;
% Computes sqrt x by Newton's method, avoiding magnitude problems.
if x<0 then terrlst(x,'sqrt) else
begin scalar trv,nx,g,l,o,c1,c2,f1; f1 := nx := o := 1.0;
if (x := float x)=0.0 or x=o then return x;
if x<o then <<x := o/x; l := t>>;
c1 := if !!decflt then 10000.0 else 8192.0; c2 := c1*c1;
while x>c2 do <<x := x/c2; f1 := f1*c1>>;
loop: trv := nx; nx := (trv + x/trv)/2;
if g and nx>=trv then go to ret;
g := t; go to loop;
ret: nx := nx*f1; return if l then o/nx else nx end;
symbolic procedure cbrt x;
begin scalar s,l,o,g,trv,nx,c1,c2,f1; f1 := nx := o := 1.0;
if (x := float x)=0.0 or abs x=o then return x
else if x<0 then x := -x else s := t;
if x<o then <<x := o/x; l := t>> else if x=o then go to ret;
c1 := if !!decflt then 1000.0 else 1024.0; c2 := c1*c1*c1;
while x>c2 do <<x := x/c2; f1 := f1*c1>>;
loop: trv := nx; nx := trv/1.5+x/(trv*trv*3);
if g and nx>=trv then go to ret;
g := t; go to loop;
ret: nx := nx*f1; if l then nx := o/nx;
return if s then nx else -nx end;
symbolic procedure hypot(p,q);
% Hypot(p,q)=sqrt(p*p+q*q) but avoids intermediate overflow.
begin scalar r;
if p<0 then p := -p;
if q<0 then q := -q;
if zerop p then return q
else if zerop q then return p
else if p<q then <<r := p; p := q; q := r>>;
if p+q=p then return p else r := q/p;
return if r<!!epsqrt then p else p*sqrt(1+r*r) end;
symbolic procedure fix2 x; if fixp x then x else fl2int x;
symbolic procedure floor x;
% Returns the largest integer less than or equal to x
% (i.e. the "greatest integer" function.)
% Note the trickiness to compensate for fact that (unlike APL's
% "FLOOR" function) FIX truncates towards zero.
% A definition of fix(x+sgn(x)*0.5) has also been suggested.
if fixp x then x
else (if x = float n then n else if x >= 0 then n else n-1)
where n = fix2 x;
symbolic procedure ceiling x;
% Returns the smallest integer greater than or equal to X.
% Note the trickiness to compensate for fact that (unlike APL's
% "FLOOR" function) FIX truncates towards zero.
if fixp x then x
else (if x = float n then n else if x >= 0 then n+1 else n)
where n = fix2 x;
symbolic procedure round x;
% Rounds to the closest integer.
% Kind of sloppy -- it's biased when the digit causing rounding is a
% five. (Changed to work properly for X<0. SLK)
if fixp x then x else if x<0 then -round(-x) else floor(x+0.5);
symbolic procedure rounddec (x,p);
% Rounds x to p decimal places, unless x must already be an integer.
if abs x>=!!flint then x
else begin scalar xl,xr,sc;
sc := 10.0**p; xl := fix(x := x*sc); xr := x-xl;
if x>0 and xr>=0.5 then xl := xl+1;
if x<0 and xr<-0.5 then xl := xl-1;
return xl/sc
end;
global '(log2 sq2 sq2!-1 logsq2 logten log1000 log1e9);
sq2 := sqrt 2.0; sq2!-1 := 1/(1+sq2);
symbolic procedure log x;
begin scalar s,lx; integer p;
if (x := float x)<=0.0 then terrlst(x,'log)
else if x-1<0 then x := 1/x else s := t;
lx := 0.0;
while x>1.0e9 do <<x := x/1.0e9; lx := lx+log1e9>>;
while x>1000 do <<x := x/1000; lx := lx+log1000>>;
while x>10 do <<x := x/10; lx := lx+logten>>;
while x>2 do <<x := x/2; lx := lx+log2>>;
if x>sq2 then <<x := x/sq2; lx := lx+logsq2>>;
lx := lx+sclogx!-1(x-1);
return if s then lx else -lx end;
symbolic procedure sclogx!-1 x;
begin scalar tx,px,lx,st,sl; integer p;
tx := px := x; p := 1; lx := 0.0;
st := x*(1-x/2);
while st+abs tx>st do
<<sl := tx . sl; tx := (px:=-px*x)/(p := p+1)>>;
for each i in sl do lx := lx+i;
return lx end;
log2 := 2*(logsq2 := sclogx!-1 sq2!-1);
log1e9 := 3*(log1000 := 3*(logten := log 10.0));
global '(pilist); remflag('(pilist),'reserved);
global(pilist := '(pii pii2 pii3 pii4 pii6 twopi !!rad2deg !!deg2rad));
flag('(pilist),'reserved);
symbolic procedure atan x;
begin scalar arg,term,termp,trv,s,g,y;
integer p;
if (x := cflot x)< 0 then x := -x else s := t;
if x > 1 then x:=1/x else g:=t;
if x < !!epsqrt then go to quad;
term := if x<0.43 then (arg := x)
else (arg := x/(1+sqrt(1+x*x)))*2;
arg := arg*arg; p := 1; trv := (termp := term)*(1-arg/3);
while trv+abs termp >trv do
<<y := termp . y; termp := (term := -term*arg)/(p := p+2)>>;
x := 0;
for each i in y do x := x+i;
quad: if not g then x := pii2-x;
if not s then x := -x;
return x end;
symbolic procedure atand x; !!rad2deg * atan x;
remflag(pilist,'reserved);
twopi := 2*(pii := 2*(pii2 := 2*(pii4:=atan 1.0)));
pii3 := 2*(pii6 := pii2/3);
!!deg2rad:=pii4/45; !!rad2deg:=45/pii4;
flag(pilist,'reserved);
fluid '(!*ddf!* !*df!* !*sf!* !*qf!*);
symbolic procedure sin x;
begin scalar !*sf!*,!*qf!*;integer p;
% test for 90 deg -> 1.0
x := !!scalsintan(x,t);
if !*qf!* then <<x := 1.0; go to ret>>;
% for x>45, compute cos of complement, else compute sin.
if x>pii4 then x := pii2-x else p := 1;
x := !!sints(x,p);
ret: return if !*sf!* then x else -x end;
symbolic procedure sind x;
begin scalar !*sf!*,!*qf!*;integer p;
% test for 90 deg -> 1.0
x := !!scalsintand(x,t);
if !*qf!* then <<x := 1.0; go to ret>>;
% for x>45, compute cos of complement, else compute sin.
if x>45.0 then x := 90.0-x else p := 1;
x := !!sints(x*!!deg2rad,p);
ret: return if !*sf!* then x else -x end;
symbolic procedure tan x;
begin scalar y,inv,!*sf!*,!*qf!*;
y:=x; x:= !!scalsintan(x,nil);
if !*qf!* then terrlst(y,'tan);
if x>pii4 then x := pii2-x else inv := t;
% For scaled x>45, compute cot else compute tan.
if x>!!epsqrt then <<x := !!sints(x,1); x := x/sqrt(1-x*x)>>;
if not inv then x := 1/x;
return if !*sf!* then -x else x end;
symbolic procedure tand x;
begin scalar y,inv,!*sf!*,!*qf!*;
y:=x; x:= !!scalsintand(x,nil);
if !*qf!* then terrlst(y,'tand);
if x>45.0 then x := 90.0-x else inv := t;
% For scaled x>45, compute cot else compute tan.
x := x*!!deg2rad;
if x>!!epsqrt then <<x := !!sints(x,1); x := x/sqrt(1-x*x)>>;
if not inv then x := 1.0/x;
return if !*sf!* then -x else x end;
global '(max!-trig!-fact); max!-trig!-fact := 10**(!!nfpd/2);
fluid '(!:prec!:);
symbolic procedure !!scalsintan(x,w);
% x is scaled to 0<=x<=90 deg, with !*sf!* = {sin>0 or tan<0}.
% w true for sin, false for tan.
%modified to avoid infinite loop for large x, after Rainer Schoepf's
%suggestion, adjusted so degree and radian input agrees- SLK.
begin scalar xf,x0;
if x<0 then x := -x else !*sf!* := t; x0 := x;
if (xf := fix2(x/twopi))>max!-trig!-fact then
terrlst(if !*sf!* then x else -x,if w then 'sin else 'tan);
x := x - float xf * twopi;
if x>pii then (if w then
<<x := twopi-x; !*sf!* := not !*sf!*>> else x := x-pii);
if x>pii2 then x:=pii-x else if not w then !*sf!*:=not !*sf!*;
!*qf!* := x>=pii2;
% the remaining tests and scaling are done separately by sin and tan
if x<x0/10.0**(!:prec!:-3) then x := 0.0;
return x end;
symbolic procedure !!scalsintand(x,w);
% x is scaled to 0<=x<=90 deg, with !*sf!* = {sin>0 or tan<0}.
% w true for sin, false for tan.
%modified to avoid infinite loop for large x, after Rainer Schoepf's
%suggestion, adjusted so degree and radian input agrees- SLK.
begin scalar xf,x0;
if x<0 then x := -x else !*sf!* := t; x0 := x;
if (xf := fix2(x/360.0))>max!-trig!-fact then
terrlst(if !*sf!* then x else -x,if w then 'sin else 'tan);
x := x - float xf * 360.0;
if x>180.0 then (if w then
<<x := 360.0-x; !*sf!* := not !*sf!*>> else x := x-180.0);
if x>90.0 then x:=180.0-x else if not w then !*sf!*:=not !*sf!*;
!*qf!* := x>=90.0;
% the remaining tests and scaling are done separately by sin and tan
if x<x0/10.0**(!:prec!:-3) then x := 0.0;
return x end;
symbolic procedure !!sints (x,p);
% Does the actual computation of the sin or cos series.
begin scalar sl,sq,st,term;
x := float x;
if x<!!epsqrt then return if p=1 then x else 1.0;
sq := x*x;
st:=if p=1 then (term:=x)*(1-sq/6) else (term:=1.0)*(1-sq/2);
while st + abs term > st do
<<sl:=term . sl;
p:=p+2.0;
term:=-term*sq/((p-1.0)*p) >>;
x:=0.0;
for each i in sl do x:=x+i;
return x end;
symbolic procedure !!sinhts x;
% Does the actual computation of the sinh for 0<x<0.91.
begin scalar p,sl,sq,st,term;
if x<2*!!epsqrt then return x;
st:= (term := x)*(1-(sq := x*x)/6); p := 1;
while st + term > st do
<<sl:=term . sl;
p:=p+2.0;
term:=term*sq/((p-1.0)*p) >>;
x := 0.0;
for each i in sl do x := x+i;
return x end;
global '(ee);
symbolic procedure exp v;
begin scalar d,nr,mr,fr,st;integer p,ip;
mr := fr := 1.0; v := float v;
if abs v>1 then <<ip := fix2 v; v := v-ip; fr := ee**ip>>;
if abs v>0.5 then v := v/2 else d := t;
if v=0.0 then go to ret;
st := mr+v;
while st+abs mr > st do
<<nr := mr . nr; mr := mr*v/(p := p+1)>>;
mr := 0.0;
for each i in nr do mr := mr+i;
ret: if not d then mr := mr*mr;
return fr*mr end;
remflag('(ee),'reserved); ee := exp 1.0; flag('(ee),'reserved);
put('expt,'number!-of!-args,2);
remflag('ee . pilist,'reserved);
symbolic procedure expt(x,y);
% Computes x**y in for argument sets that yield real values.
% Calls exp or sqrt where appropriate.
begin scalar x1,y1,s,prd;integer p,ip; x1 := x; y1 := y;
if zerop y then go to zry
else if zerop x then go to zrx;
if fixp(ip := y) or zerop(y-(ip := fix2 y)) then go to int
% Delay inversion only if exponent is integer.
else if y<0 then <<x := 1.0/x; y := -y>>;
if y-1.0=0.0 then return x
else if x=ee then return exp y
else if (prd := rounddec(y,!!nfpd-3))=0.5
then return sqrt float x
else if prd=0.25 then return sqrt sqrt float x
else if prd=1.5 then return x*sqrt float x
else if x>0 then go to ret;
x := -x;p := 1;
repeat ip := fix(prd := rounddec(y*(p := p+2),!!nfpd-3))
until p>99 or zerop(ip-prd);
if p>99 then error (0,list (x,"**",y," not real"))
else if remainder(ip,2)=1 then s := t;
ret: x := exp(y*log float x); return if s then -x else x;
zrx: if y>0 then return if floatp x or floatp y then 0.0 else 0
else error(0,"divide by zero in EXPT");
zry: if zerop x then error(0,"0**0 undefined")
else return if floatp x or floatp y then 1.0 else 1;
int: p := 1; if ip<0 then <<ip := -ip; x := 1.0/x>>;
loop: y := ip; if y-2*(ip := ip/2) neq 0 then p := p*x;
if ip=0 then return p;
x := x*x; go to loop end;
symbolic procedure float!-constants;
<<find!!decflt(); !!epsqrt := sqrt(10.0**(-1-!!nfpd));
!!flint := 10.0**!!nfpd; !!shbinflp();
sq2 := sqrt 2.0; sq2!-1 := 1/(1+sqrt 2.0);
log2 := 2*(logsq2 := sclogx!-1 sq2!-1);
log1e9 := 3*(log1000 := 3*(logten := log 10.0));
remflag(pilist,'reserved);
twopi := 2*(pii := 2*(pii2 := 2*(pii4:=atan 1.0)));
pii3 := 2*(pii6 := pii2/3);
!!deg2rad:=pii4/45; !!rad2deg:=45/pii4;
flag(pilist,'reserved);
max!-trig!-fact := 10**(!!nfpd/2);
remflag('(ee),'reserved); ee := exp 1.0; flag('(ee),'reserved);
"constants recomputed">>;
flag('ee . pilist,'reserved);
%lisp float!-constants();
%**********************************************************************
%** Functions derived from basic functions **
%**********************************************************************
symbolic procedure cos x; sin(pii2-x);
symbolic procedure cot x; tan(pii2-x);
symbolic procedure sec x; 1/cos x;
symbolic procedure csc x; 1/sin x;
symbolic procedure acot x; pii2-atan x;
symbolic procedure asin x;
if abs x<1 then
atan(if abs x<!!epsqrt then x else x/sqrt(1-x*x))
else if abs x>1 then terrlst (x,'asin)
else if x>0 then pii2 else -pii2;
symbolic procedure acos x; pii2-asin x;
symbolic procedure acsc x;
if abs x>=1 then asin(1.0/x) else terrlst(x,'acsc);
symbolic procedure asec x;
if abs x<1 then terrlst(x,'asec) else pii2-asin(1.0/x);
symbolic procedure cosd x; sind(90.0-x);
symbolic procedure cotd x; tand(90.0-x);
symbolic procedure secd x; 1/cosd x;
symbolic procedure cscd x; 1/sind x;
symbolic procedure acotd x; 90.0-atand x;
symbolic procedure asind x; !!rad2deg * asin x;
symbolic procedure acosd x; 90.0-asind x;
symbolic procedure acscd x;
if abs x>=1 then asind(1.0/x) else terrlst(x,'acscd);
symbolic procedure asecd x;
if abs x<1 then terrlst(x,'asecd) else 90.0-asind(1.0/x);
symbolic procedure sinh x;
begin scalar s;
if x<0.0 then x:=-x else s:=t;
if (x := float x)<0.91 then <<x := !!sinhts x; go to ret>>;
x := exp(-x); x := (1.0/x-x)/2;
ret: return if s then x else -x end;
symbolic procedure cosh x; <<x := exp(-abs x),(x+1.0/x)/2>>;
symbolic procedure tanh x;
if x<0.0 then -tanh(-x) else
<<x := exp(-2.0*x); (1.0-x)/(1.0+x)>>;
symbolic procedure coth x;
if x<0.0 then -coth(-x) else
<<x := exp(-2.0*x); (1.0+x)/(1.0-x)>>;
symbolic procedure asinh x; begin scalar s;
if x<0 then x:=-x else s:=t;
x:=if x<!!epsqrt then x else log (x+if x<2 then sqrt(x*x+1)
else if 1/x<!!epsqrt then x else x*sqrt(1+1/(x*x)));
return if s then x else -x end;
symbolic procedure acosh x; if x<1 then terrlst(x,'acosh)
else log (x+if 1/x<!!epsqrt then x else x*sqrt(1-1/(x*x)));
symbolic procedure atanh x; if abs x>=1 then terrlst(x,'atanh)
else if abs x<!!epsqrt then x else 0.5*log((1+x)/(1-x));
symbolic procedure acoth x;
if abs x<=1 then terrlst(x,'acoth) else atanh (1.0/x);
symbolic procedure sech x;1/cosh x;
symbolic procedure csch x;1/sinh x;
symbolic procedure asech x;
if x<=0 or x>1 then terrlst(x,'asech) else acosh (1.0/x);
symbolic procedure acsch x;
if (x:= float x)=0.0 then terrlst(x,'acsch) else asinh(1/x);
symbolic procedure ln x; log x;
symbolic procedure log10 x;
if x>0 then log x/logten else terrlst(x,'logten);
symbolic procedure logb (x,b); %log x to base b;
begin scalar a,s; a:=x>0; s:=not(b<=0 or zerop(b-1));
if a and s then return log x/log b
else terrlst((if a then list ('base,b)
else if s then list('arg,x) else list(x,b)),'logb) end;
symbolic procedure atan2(y,x);
if zerop x then pii2*sgn y else
<<(if x>0 then a else if y<0 then a-pii else a+pii)
where a=atan(y/x)>>;
symbolic procedure atan2d(y,x);
if zerop x then 90.0*sgn y else
<<(if x>0 then a else if y<0 then a-180.0 else a+180.0)
where a=!!rad2deg*atan(y/x)>>;
% A numerical factorial function.
symbolic procedure factorial n;
if not fixp n or n<0
then rerror(arith,50,list(n,"invalid factorial argument"))
else nfactorial n;
symbolic procedure nfactorial n;
% Numerical factorial function. It is assumed that n is numerical
% and non-negative.
if n>20 then fac!-part(1,n)
else begin scalar m;
m:=1;
for i:=1:n do m:=m*i;
return m;
end;
symbolic procedure fac!-part (m,n);
if m=n then m
else if m=n-1 then m*n
else (fac!-part(m,p)*fac!-part(p+1,n)) where p=(m+n)/2;
endmodule;
module rounded; % *** Support for Arbitrary Rounded Arithmetic.
% Authors: Anthony C. Hearn and Stanley L. Kameny.
% Copyright (c) 1991 RAND. All rights reserved.
Comment this module defines a rounded object as a list with two fields:
(<tag>.<structure>).
The <structure> depends on the precision. It is either a floating point
number or the stripped bfloat (mt . ep);
fluid '(!:prec!: !:print!-prec!: minprec!# rootacc!#);
global '(bfone!* epsqrt!*);
fluid '(dmode!* !*bfspace !*numval !*roundbf !*!*roundbf !*norndbf);
global '(domainlist!* !!nfpd !!flprec !!rdprec mxflbf!! mnflbf!!);
global '(!!plumax !!plumin !!timmax !!timmin !!maxflbf
!!minflbf !!fleps1 !!fleps2 !!flint !!maxflexp logten
!!maxarg);
switch rounded;
symbolic smacro procedure ep!: u; cddr u;
symbolic smacro procedure mt!: u; cadr u;
%Set value for !!flprec. It never changes.
!!flprec := !!nfpd-if !!decflt then 1 else 3;
symbolic procedure logfp x;
% floating log of x**(1/n) using bfloat logic as boost.
(log(m/10.0**p)+(p+ep!: x)*logten)
where p=length explode m-1 where m=mt!: x;
symbolic procedure roundconstants;
<<!!plumax := 0.5*10.0**!!maxflexp;
!!minflbf := invbf(!!maxflbf := '!:bf!: . 1 . !!maxflexp);
% plumin must be large enough to avoid underflow from difference.
!!plumin := 10.0**!!flprec/!!plumax;
!!timmin := 1/(!!timmax := 0.7*sqrt(10.0**!!maxflexp));
!!maxarg := logfp !!maxflbf>>;
switch bfspace,numval,roundbf,norndbf; !*bfspace := !*numval := t;
put('roundbf,'simpfg,'((t (roundbfon)) (nil (roundbfoff))));
symbolic procedure roundbfon; !*!*roundbf := t;
symbolic procedure roundbfoff; !*!*roundbf := !!rdprec > !!flprec;
global '(!!decflt rd!-tolerance!* cr!-tolerance!* yy!! bfz!*);
% put('rounded,'package!-name,'arith); % Use if ARITH autoloaded.
domainlist!* := union('(!:rd!:),domainlist!*);
put('rounded,'tag,'!:rd!:);
put('!:rd!:,'dname,'rounded);
flag('(!:rd!:),'field);
put('!:rd!:,'i2d,'i2rd!*);
put('!:rd!:,'minusp,'rd!:minusp);
put('!:rd!:,'plus,'rd!:plus);
put('!:rd!:,'times,'rd!:times);
put('!:rd!:,'difference,'rd!:difference);
put('!:rd!:,'quotient,'rd!:quotient);
put('!:rd!:,'zerop,'rd!:zerop);
put('!:rd!:,'onep,'rd!:onep);
put('!:rd!:,'prepfn,'rd!:prep);
put('!:rd!:,'prifn,'rd!:prin);
put('!:rd!:,'minus,'rd!:minus);
put('!:rd!:,'!:rn!:,'!*rd2rn);
put('!:rn!:,'!:rd!:,'!*rn2rd);
symbolic procedure round!* x;
% Returns actual number representation, as either float or bfloat.
if atom(x := cdr x) then x else '!:bf!: . x;
symbolic smacro procedure smallfloatp u;
% Returns true if <structure> is a small rounded.
atom u;
symbolic smacro procedure preci!: nmbr;
% This function counts the precision of a number "n".
% NMBR is a BIG-FLOAT representation of "n".
length explode abs mt!: nmbr$
symbolic smacro procedure make!:bf(mt,ep); '!:bf!: . (mt . ep)$
symbolic procedure mkround u; '!:rd!: . if atom u then u else cdr u;
symbolic procedure roundbfp; !*roundbf or !!rdprec > !!flprec;
symbolic procedure bftrim!: v;
normbf round!:mt(v,if !:print!-prec!:
then !:print!-prec!:
else !:prec!:-2);
%remflag('(precision print!-precision),'lose);
% probably not needed any more, since it doesn't get flagged lose.
symbolic procedure print!-precision n;
% Set the system printing precision !:print!-prec!:.
% Returns previous value.
begin scalar oldprec;
if n=0 then return !:print!-prec!:;
if n<0 then
<< oldprec := !:print!-prec!:;
!:print!-prec!: := nil;
return oldprec >>;
if n > !:prec!: then
<< msgpri(nil,"attempt to set print!-precision greater than",
"precision ignored",nil,nil);
return nil >>;
oldprec := !:print!-prec!:;
!:print!-prec!: := n;
return oldprec
end;
symbolic procedure print!_precision n;
% Alternative name.
print!-precision n;
symbolic procedure precision n;
% Set the system precision !!rdprec, bfloat precision !:prec!:,
% and rd!:onep tolerance. Returns previous value.
<<if not numberp n or n<0
then rerror(arith,2,"positive number required");
if (n := fix n)>0 then rootacc!# := max(n,6);
precision1(n,t)>>;
% share !!fleps1,!!fleps2,!:prec!:;
symbolic procedure precision1(n,bool);
begin scalar oldprec;
if n=0 then return !!rdprec;
if bool then rmsubs(); % So that old results are resimplified.
oldprec := !!rdprec;
!:prec!: := (!!rdprec := max(n,minprec!#))+2;
if !:print!-prec!: and n < !:print!-prec!:+2
then !:print!-prec!: := nil; %unset
epsqrt!* := '!:bf!: . (1 . -!:prec!:/2);
rd!-tolerance!* := '!:bf!: . 1 . (2-!:prec!:);
cr!-tolerance!* := '!:bf!: . 1 . (2*(2-!:prec!:));
if getd 'expt!* then
algebraic <<!!fleps1 := 10.0**(2-!:prec!:);
!!fleps2 := !!fleps1**2>>
else <<!!fleps1 := 10.0**(2-!:prec!:);
!!fleps2 := !!fleps1**2>>;
!*!*roundbf := !!rdprec > !!flprec or !*roundbf;
return oldprec end;
flag('(precision),'opfn); % Symbolic operator precision.
flag('(print!-precision),'opfn); % Symbolic operator print!-precision.
flag('(print!_precision),'opfn); % Symbolic operator print_precision.
symbolic procedure !*rd2rn n;
% Converts a rounded number N into a rational to the system precision.
mkrn(car x,cdr x) where x=realrat bftrim!: bfloat round!* n;
symbolic procedure !*rn2rd u;
% Converts the (tagged) rational u/v into a (tagged) rounded
% number to the system precision, after testing to number
mkround chkrn!* r2bf cdr u;
minprec!# := min(6,!!flprec-2);
precision1(!!flprec,nil); % Initial value = effective float prec.
symbolic procedure precmsg pr;
if pr>!!rdprec then
<<msgpri(nil,"precision increased to",pr,nil,nil);
precision1(pr,t)>>;
symbolic procedure rd!:simp u;
if eqcar(u,0) then nil ./ 1
else if u = '(10 . -1) then 1 ./ 1
else if null dmode!* or dmode!* eq '!:gi!:
then (if eqcar(x,'!:rn!:) then cdr x else x ./ 1)
where x = !*rd2rn('!:rd!: . u)
else if dmode!* memq '(!:rd!: !:cr!:)
then (mkround convprec!* u) ./ 1 % Must call convprec!*, since
% precision may have changed.
else (if y then !*d2q apply1(y,'!:rd!: . u)
else dmoderr('!:rd!:,dmode!*))
where y = get('!:rd!:,dmode!*);
put('!:rd!:,'simpfn,'rd!:simp);
symbolic procedure rndbfon; if not !*norndbf then
<<!*!*roundbf := t;
if !:prec!:<!!flprec+3 then
<<!*roundbf := t;
lprim "ROUNDBF turned on to increase accuracy">>>>;
symbolic procedure i2rd!* u;
% Converts integer U to tagged rounded form.
mkround chkint!* u;
symbolic procedure chkint!* u;
if !*!*roundbf then bfloat u else
((if x<65535 then float u else if x<10**!!maxflexp then cflot u
else <<rndbfon(); bfloat u>>) where x=abs u);
mnflbf!! := invbf(mxflbf!! := '!:bf!: . (1 . 800));
symbolic procedure erbdck r;
errorp r or not mt!:(r := abs!: bfloat r)=0 and
(greaterp!:(r,mxflbf!!) or greaterp!:(mnflbf!!,r));
copyd('errorp1,if !!decflt then 'erbdck else 'errorp);
symbolic procedure chkrn!* u;
if !*!*roundbf then u else bf2flck u;
symbolic procedure bf2flck u;
if !*!*roundbf then u
else if mt!: u=0 then 0.0 else
((if not grpbf(!!minflbf,r) and not grpbf(r,!!maxflbf)
then bf2flr u
else <<rndbfon(); u>>) where r := abs!: u);
symbolic procedure convchk x;
if !*!*roundbf then if atom x then bfloat x else x
else if atom x then x else bf2flck x;
symbolic procedure convprec!* u;
convchk if atom u then u else '!:bf!: . u;
symbolic procedure convprec u; convchk round!* u;
symbolic procedure rd!:minusp u; bfminusp round!* u;
symbolic procedure convprc2(u,v);
<<u := convchk round!* u; yy!! := convchk round!* v;
if !*roundbf then <<yy!! := bfloat yy!!; bfloat u>> else u>>;
symbolic procedure rdzchk(u,x,y);
if atom u then
if u=0.0 or x>0.0 and y>0.0 or x<0.0 and y<0.0 then u
else if abs u<(abs x)*!!fleps1 then 0.0 else u
else
if mt!: u=0 or mt!: x>0 and mt!: y>0 or mt!: x<0 and mt!: y<0 then u
else if lessp!:(abs!: u,times!:(abs!: x,rd!-tolerance!*)) then bfz!*
else u;
symbolic procedure rd!:plus(u,v);
(if not !*!*roundbf and atom cdr u and atom cdr v
and (z := safe!-fp!-plus(cdr u,cdr v)) then '!:rd!: . z else
begin scalar x,y;
x := convprc2(u,v); y := yy!!;
u := if not atom x then plubf(x,y) else
<<z := errorset!*(list('plus2,mkquote x,mkquote y),nil);
if errorp1 z
then <<rndbfon(); plubf(x := bfloat x,y := bfloat y)>>
else car z>>;
return mkround rdzchk(u,x,y) end) where z=nil;
symbolic procedure rd!:difference(u,v);
(if not !*!*roundbf and atom cdr u and atom cdr v
and (z := safe!-fp!-plus(cdr u,-cdr v)) then '!:rd!: . z else
begin scalar x,y;
x := convprc2(u,v); y := yy!!;
u := if not atom x then difbf(x,y) else
<<z := errorset!*(list('difference,mkquote x,mkquote y),nil);
if errorp1 z
then <<rndbfon(); difbf(x := bfloat x,y := bfloat y)>>
else car z>>;
return mkround rdzchk(u,x,if atom y then -y else minus!: y) end)
where z=nil;
symbolic procedure rd!:times(u,v);
(if not !*!*roundbf and atom cdr u and atom cdr v
and (z := safe!-fp!-times(cdr u,cdr v)) then '!:rd!: . z else
begin scalar x,y;
x := convprc2(u,v); y := yy!!;
return mkround if not atom x then timbf(x,y) else
<<z := errorset!*(list('times2,mkquote x,mkquote y),nil);
if errorp1 z then <<rndbfon(); timbf(bfloat x,bfloat y)>>
else car z>> end) where z=nil;
symbolic procedure rd!:quotient(u,v);
(if not !*!*roundbf and atom cdr u and atom cdr v
and (z := safe!-fp!-quot(cdr u,cdr v)) then '!:rd!: . z else
begin scalar x,y;
x := convprc2(u,v); y := yy!!;
if atom x and zerop y then rdqoterr();
return mkround if not atom x then
if mt!: y=0 then rdqoterr() else divbf(x,y)
else
<<z := errorset!*(list('quotient,mkquote x,mkquote y),nil);
if errorp1 z then <<rndbfon(); divbf(bfloat x,bfloat y)>>
else car z>> end) where z=nil;
symbolic procedure rdqoterr; error(0,"zero divisor in quotient");
symbolic procedure safe!-fp!-plus(x,y);
if zerop x then y else if zerop y then x else
begin scalar u;
if x>0.0 and y>0.0 then
if x<!!plumax and y<!!plumax then go to ret else return nil;
if x<0.0 and y<0.0 then
if -x<!!plumax and -y<!!plumax then go to ret else return nil;
if abs x<!!plumin and abs y<!!plumin then return nil;
ret: return
if (u := plus2(x,y))=0.0
or x>0.0 and y>0.0 or x<0.0 and y<0.0 then u
else if abs u<(abs x)*!!fleps1 then 0.0 else u end;
symbolic procedure safe!-fp!-times(x,y);
if zerop x or zerop y then 0.0
else if x=1.0 then y else if y=1.0 then x else
begin scalar u,v; u := abs x; v := abs y;
if u>=1.0 and u<=!!timmax then
if v<=!!timmax then go to ret else return nil;
if u>!!timmax then if v<=1.0 then go to ret else return nil;
if u<1.0 and u>=!!timmin then
if v>=!!timmin then go to ret else return nil;
if u<!!timmin and v<1.0 then return nil;
ret: return times2(x,y) end;
symbolic procedure safe!-fp!-quot(x,y);
if zerop y then rdqoterr()
else if zerop x then 0.0 else if y=1.0 then x else
begin scalar u,v; u := abs x; v := abs y;
if u>=1.0 and u<=!!timmax then
if v>=!!timmin then go to ret else return nil;
if u>!!timmax then if v>=1.0 then go to ret else return nil;
if u<1.0 and u>=!!timmin then
if v<=!!timmax then go to ret else return nil;
if u<!!timmin and v>1.0 then return nil;
ret: return quotient(x,y) end;
symbolic procedure rd!:zerop u; bfzp round!* u;
symbolic procedure rd!:minus u; mkround bfminus round!* u;
symbolic procedure rd!:onep u;
% We need the tolerance test since some LISPs (e.g. PSL) can print
% a number as 1.0, but it doesn't equal 1.0!
(if atom x then abs(1.0-x)<!!fleps1 else equal!:(bfone!*,bftrim!: x))
where x = round!* u;
% Since decimal input -> :rd: in all dmodes, dmode!* must be used to
% determine whether to round to current precision, but input never gets
% truncated, since precision is always increased at input time.
% to avoid inaccuracies in floating point representation, rd!:prep
% returns values in bfloat format.
symbolic procedure rd!:prep u;
if rd!:onep u then 1
else if rd!:onep rd!:minus u then -1 else rdprep1 u;
symbolic procedure rdprep1 u;
mkround(
if atom x then
if not dmode!* memq '(!:rd!: !:cr!:) or !*!*roundbf
then round!:mt(bfloat x,min(!:prec!:,!!nfpd))
else if !:prec!:>=!!nfpd then x
else bf2flr round!:mt(bfloat x,!:prec!:)
else round!:mt('!:bf!: . x,!:prec!:))
where x=cdr u;
symbolic procedure rd!:prin u;
% Printed output is rounded to 2 fewer places than internal value.
bfprin!: normbf bftrim!: bfloat round!* u;
% Support for REDUCE 3.3 switches BIGFLOAT and FLOAT.
symbolic procedure onbigfloat;
% Permit BIGFLOAT as a switch for compatibility.
begin
lprim "Please use ROUNDED instead";
!*roundbf := t; on1 'rounded
end;
symbolic procedure offbigfloat; <<!*roundbf := nil; off1 'rounded>>;
symbolic procedure onfloat;
% Permit FLOAT as a switch for compatibility.
begin
lprim "Please use ROUNDED instead";
if precision 0 neq !!flprec then
<<lprim list("Precision changed to",!!flprec);
precision1(!!flprec,t)>>;
!*roundbf := nil; on1 'rounded
end;
put('float,'simpfg,'((t (onfloat)) (nil (off1 'rounded))));
put('bigfloat,'simpfg,'((t (onbigfloat)) (nil (offbigfloat))));
switch bigfloat,float;
initdmode 'rounded;
endmodule;
module comprd; % *** Support for Complex Rounded Arithmetic.
% Authors: Anthony C. Hearn and Stanley L. Kameny.
% Copyright (c) 1989 The RAND Corporation. All rights reserved.
Comment this module defines a complex rounded as:
(<tag>. (<structure> . <structure>>) or ('!:crd!: . (rl . im))
The <structure> depends on the precision. It is either a floating point
number or the stripped bfloat (mt . ep);
fluid '(!:prec!:);
global '(bfone!* epsqrt!*);
fluid '(dmode!* !*bfspace !*numval !*prin!# !*roundbf !*!*roundbf);
global '(cr!-tolerance!* domainlist!* !!nfpd !!flprec !!rdprec bfz!*
yy!!);
domainlist!* := union('(!:cr!:),domainlist!*);
fluid '(!*complex!-rounded);
put('complex!-rounded,'tag,'!:cr!:);
put('!:cr!:,'dname,'complex!-rounded);
flag('(!:cr!:),'field);
put('!:cr!:,'i2d,'i2cr!*);
put('!:cr!:,'plus,'cr!:plus);
put('!:cr!:,'times,'cr!:times);
put('!:cr!:,'difference,'cr!:differ);
put('!:cr!:,'quotient,'cr!:quotient);
put('!:cr!:,'zerop,'cr!:zerop);
put('!:cr!:,'onep,'cr!:onep);
put('!:cr!:,'prepfn,'cr!:prep);
put('!:cr!:,'prifn,'cr!:prin);
put('!:cr!:,'minus,'cr!:minus);
put('!:cr!:,'minusp,'cr!:minusp);
put('!:cr!:,'rationalizefn,'girationalize!:);
put('!:cr!:,'!:rn!:,'!*cr2rn);
put('!:rn!:,'!:cr!:,'!*rn2cr);
put('!:rd!:,'!:cr!:,'!*rd2cr);
put('!:cr!:,'!:rd!:,'!*cr2rd);
put('!:cr!:,'!:crn!:,'!*cr2crn);
put('!:crn!:,'!:cr!:,'!*crn2cr);
put('!:gi!:,'!:cr!:,'!*gi2cr);
put('!:cr!:,'cmpxfn,'mkcr);
put('!:cr!:,'ivalue,'mkdcrn);
put('!:cr!:,'realtype,'!:rd!:);
put('!:rd!:,'cmpxtype,'!:cr!:);
symbolic procedure cr!:minusp u;
(if atom x then zerop y and x<0 else zerop car y and car x<0)
where x=cadr u,y=cddr u;
symbolic procedure striptag u; if atom u then u else cdr u;
symbolic procedure mkcr(u,v); '!:cr!: . (striptag u) . striptag v;
symbolic smacro procedure crrl x; cadr x;
symbolic smacro procedure crim x; cddr x;
symbolic smacro procedure mt!: x; cadr x;
symbolic smacro procedure ep!: x; cddr x;
symbolic smacro procedure preci!: x; length explode abs mt!: x;
symbolic procedure gf2cr!: x;
'!:cr!: . (striptag car x) . striptag cdr x;
symbolic procedure crprcd u;
(rl . im) where rl=convprec!* crrl u,im=convprec!* crim u;
symbolic procedure crprcd2(x,y);
<<x := crprcd x; yy!! := crprcd y; x>>;
% simp must call convprec!*, since precision may have changed.
symbolic procedure cr!:simp u; (gf2cr!: crprcd u) ./ 1;
put('!:cr!:,'simpfn,'cr!:simp);
symbolic procedure mkdcr u; cri!*() ./ 1;
symbolic procedure i2cr!* u;
%converts integer U to tagged cr form.
<<u := chkint!* u; mkcr(u,if atom u then 0.0 else bfz!*)>>;
symbolic procedure trimcrrl n; trimcr crrl n;
symbolic procedure trimcr n;
bftrim!: if atom n then bfloat n else '!:bf!: . n;
symbolic procedure cr2rderr;
error(0,
"complex to real type conversion requires zero imaginary part");
symbolic procedure !*cr2rn n;
% Converts a cr number n into a rational if possible.
if bfnzp retag crim n then cr2rderr() else
<<n := realrat trimcrrl n; mkrn(car n,cdr n)>>;
symbolic procedure !*rn2cr u;
% Converts the (tagged) rational u/v into a (tagged) rounded complex
% number to the system precision.
<<u := chkrn!* r2bf cdr u; mkcr(u,if atom u then 0.0 else bfz!*)>>;
symbolic procedure !*cr2crn u;
% Converts a (tagged) cr number u into a (tagged) crn.
mkcrn(realrat trimcrrl u,realrat trimcr crim u);
symbolic procedure !*crn2cr u;
% Converts a (tagged) crn number u into a (tagged) cr.
mkcr(rl,if !*roundbf then bfloat im else im)
where rl=chkrn!* r2bf cadr u where im=chkrn!* r2bf cddr u;
symbolic procedure !*cr2rd n;
if bfnzp retag crim n then cr2rderr() else '!:rd!: . crrl n;
symbolic procedure !*rd2cr u;
mkcr(x,if atom x then 0.0 else bfz!*) where x=convprec u;
symbolic procedure !*gi2cr u;
mkcr(rl,if !*roundbf then bfloat im else im)
where rl=chkint!* cadr u where im=chkint!* cddr u;
symbolic procedure errp2 r;
errorp r or erbdck caar r or erbdck cdar r;
symbolic procedure bfrsq u;
(if atom x then x*x+y*y else plubf(timbf(x,x),timbf(y,y)))
where x=car u,y=cdr u;
symbolic procedure crzchk(u,x,y);
begin
if atom car u then
if bfrsq u<(bfrsq x)*!!fleps2 then return 0.0 . 0.0 else go to ck;
if lessp!:(bfrsq u,timbf(bfrsq x,cr!-tolerance!*)) then
return bfz!* . bfz!*;
ck: return rdzchk(car u,car x,car y) . rdzchk(cdr u,cdr x,cdr y) end;
symbolic procedure cr!:plus(u,v);
begin scalar x,y; x := crprcd2(u,v); y := yy!!;
u := if !*!*roundbf then gfplus(x,y)
else if (v := safe!-crfp!-plus(x,y)) then v else
((if errp2 r then
<<rndbfon(); gfplus(x := gf2bf x,y := gf2bf y)>> else car r)
where r=errorset(list('gfplus,mkquote x,mkquote y),nil,nil));
return gf2cr!: crzchk(u,x,y) end;
symbolic procedure cr!:differ(u,v);
begin scalar x,y; x := crprcd2(u,v); y := yy!!;
u := if !*!*roundbf then gfdiffer(x,y)
else if (v := safe!-crfp!-diff(x,y)) then v else
((if errp2 r then
<<rndbfon(); gfplus(x := gf2bf x,y := gf2bf y)>> else car r)
where r=errorset(list('gfdiffer,mkquote x,mkquote y),nil,nil));
return gf2cr!: crzchk(u,x,gfminus y) end;
symbolic procedure cr!:times(u,v);
gf2cr!:
(if !*!*roundbf then gftimes(x,yy!!)
else if (u := safe!-crfp!-times(x,yy!!)) then u else
((if errp2 r then <<rndbfon(); gftimes(gf2bf x,gf2bf yy!!)>>
else car r)
where r=errorset!*(list('gftimes,mkquote x,mkquote yy!!),nil)))
where x=crprcd2(u,v);
symbolic procedure cr!:quotient(u,v);
gf2cr!:
(if gfzerop yy!! then rdqoterr()
else if !*!*roundbf then gfquotient(x,yy!!)
else if (u := safe!-crfp!-quot(x,yy!!)) then u else
((if errp2 r then
<<rndbfon(); gfquotient(gf2bf x,gf2bf yy!!)>> else car r)
where r=errorset!*(list('gfquotient,mkquote x,mkquote yy!!),nil)))
where x=crprcd2(u,v);
symbolic procedure safe!-crfp!-plus(u,v);
(if x and y then crzchk(x . y,u,v))
where x=safe!-fp!-pl(car u,car v),y=safe!-fp!-pl(cdr u,cdr v);
symbolic procedure safe!-crfp!-diff(u,v);
(if x and y then crzchk(x . y,u,gfminus v))
where x=safe!-fp!-pl(car u,-car v),y=safe!-fp!-pl(cdr u,-cdr v);
symbolic procedure safe!-crfp!-times(u,v);
begin scalar ru,iu,rv,iv,a,b;
ru := car u; iu := cdr u; rv := car v; iv := cdr v;
if not (a := safe!-fp!-times(ru,rv)) or
not (b := safe!-fp!-times(iu,iv)) then return nil;
if not(u := safe!-fp!-pl0(a,-b)) then return nil;
u := rdzchk(u,a,-b);
if not (a := safe!-fp!-times(ru,iv)) or
not (b := safe!-fp!-times(iu,rv)) then return nil;
if not(v := safe!-fp!-pl0(a,b)) then return nil;
return u . rdzchk(v,a,b) end;
symbolic procedure safe!-crfp!-quot(u,v);
% compute u * inverse v.
begin scalar ru,iu,rv,iv,a,b,dd;
ru := car u; iu := cdr u; rv := car v; iv := cdr v;
if not (a := safe!-fp!-times(rv,rv)) or
not (b := safe!-fp!-times(iv,iv)) then return nil;
rv := rv/(dd := a+b); iv := -iv/dd;
if not (a := safe!-fp!-times(ru,rv)) or
not (b := safe!-fp!-times(iu,iv)) then return nil;
if not(u := safe!-fp!-pl0(a,-b)) then return nil;
u := rdzchk(u,a,-b);
if not (a := safe!-fp!-times(ru,iv)) or
not (b := safe!-fp!-times(iu,rv)) then return nil;
if not(v := safe!-fp!-pl0(a,b)) then return nil;
return u . rdzchk(v,a,b) end;
symbolic procedure safe!-fp!-pl(x,y);
% floating plus protect from under- and over-flows but no zero
% result check.
if zerop x then y else if zerop y then x
else if x>0 and y>0 then
if x<!!plumax and y<!!plumax then plus2(x,y) else nil
else if x<0 and y<0 then
if (-x<!!plumax and -y<!!plumax) then plus2(x,y) else nil
else if abs x<!!plumin or abs y<!!plumin then nil else plus2(x,y);
symbolic procedure safe!-fp!-pl0(x,y);
% protects floating plus against under-flow only.
if zerop x then y else if zerop y then x
else if abs x<!!plumin and abs y<!!plumin then nil else plus2(x,y);
symbolic procedure crrnda(y,a);
if bfzp y then y else
<<(if b<=0 then bfz!* else normbf round!:mt(y,b))
where b=(if !*prin!# then !:prec!:-2 else !:prec!:)-max(a-1,0)>>;
symbolic procedure crrnd y;
begin scalar a,d,b; a := car (y := gf2bf y); d := cdr y;
if bfzp a or bfzp d then return y;
b := ep!: a+preci!: a-ep!: d-preci!: d;
return crrnda(a,-b) . crrnda(d,b) end;
symbolic procedure cr!:minus u; gf2cr!: gfminus crprcd u;
symbolic procedure cr!:zerop u;
bfzp retag crrl u and bfzp retag crim u;
symbolic procedure cr!:onep u;
bfzp retag crim u and rd!:onep mkround retag crrl u;
% prep works entirely in bfloat, to avoid floating point conversion
% errors.
symbolic procedure cr!:prep u;
crprep1((rdprep1 car x) . rdprep1 cdr x)
where x=crrnd((retag crrl u) . retag crim u);
symbolic procedure crprep1 u;
(if rd!:zerop d then rd!:prep a
else if rd!:zerop a then
if rd!:minusp d
then if rd!:onep rd!:minus d then list('minus,'i)
else list('times,d,'i)
else crprimp d
else if rd!:minusp d
then list('difference,rd!:prep a,crprimp rd!:minus d)
else list('plus,rd!:prep a,crprimp d)) where a=car u,d=cdr u;
symbolic procedure crprimp u;
if rd!:onep u then 'i else list('times,u,'i);
symbolic procedure cr!:prin v;
(if atom (v := cr!:prep v)
or car v eq 'times or car v memq domainlist!*
then maprin v
else <<prin2!* "("; maprin v; prin2!* ")">>)
where !*prin!#=t;
initdmode 'complex!-rounded;
endmodule;
module rdelem; % Elementary functions in rounded domain.
fluid '(!:prec!: !*!*roundbf);
global '(bfz!* bfone!* bften!* bfhalf!* !:180!* !:bf1!.5!* !:bf3!*
!:bf60!* epsqrt!* bftwo!* pii !!flprec !!rdprec !!decflt !!shbinfl
pi!/180 !180!/pi ee !!maxarg);
pi!/180 := pii/180; !180!/pi := 180/pii;
bfhalf!* := '!:bf!: . (5 . -1); !:180!* := '!:bf!: . (18 . 1);
bfone!* := '!:bf!: . (1 . 0); bften!* := '!:bf!: . (1 . 1);
!:bf1!.5!* := '!:bf!: . (15 . -1); !:bf3!* := '!:bf!: . (3 . 0);
bftwo!* := '!:bf!: . (2 . 0); !:bf60!* := '!:bf!: . (6 . 1);
fluid '(!*numval);
deflist('((exp rdexp!*) (expt rdexpt!*) (log rdlog!*) (sin rdsin!*)
(cos rdcos!*) (tan rdtan!*) (asin rdasin!*) (acos rdacos!*)
(atan rdatan!*) (sqrt rdsqrt!*) (sinh rdsinh!*) (cosh rdcosh!*)
(sec rdsec!*) (csc rdcsc!*) (cot rdcot!*) (tanh rdtanh!*)
(coth rdcoth!*) (sech rdsech!*) (csch rdcsch!*) (asinh rdasinh!*)
(acosh rdacosh!*) (acot rdacot!*) (asec rdasec!*) (acsc rdacsc!*)
(atanh rdatanh!*) (acoth rdacoth!*) (asech rdasech!*)
(acsch rdacsch!*) (logb rdlogb!*) (log10 rdlog10!*) (ln rdlog!*)
(cbrt rdcbrt!*) (atan2 rdatan2!*) (hypot rdhypot!*)
(sind rdsind!*) (cosd rdcosd!*) (asind rdasind!*) (acosd rdacosd!*)
(tand rdtand!*) (cotd rdcotd!*) (atand rdatand!*) (acotd rdacotd!*)
(secd rdsecd!*) (cscd rdcscd!*) (asecd rdasecd!*) (acscd rdacscd!*)
(atan2d rdatan2d!*) (deg2rad deg2rad!*) (rad2deg rad2deg!*)
(deg2dms deg2dms!*) (rad2dms rad2dms!*) (dms2deg dms2deg!*)
(dms2rad dms2rad!*) (norm rdnorm!*) (arg rdarg!*)
(e rde!*) (pi rdpi!*)),
'!:rd!:);
for each n in '(exp sin cos tan asin acos atan sinh cosh % log
sec csc cot tanh coth sech csch asinh acosh acot asec acsc atanh acoth
asech acsch logb log10 ln cbrt atan2 hypot sind cosd asind acosd tand
cotd atand acotd secd cscd asecd acscd atan2d deg2rad rad2deg deg2dms
rad2dms dms2deg dms2rad norm arg argd) do put(n,'simpfn,'simpiden);
deflist('((dms2deg!* simpdms) (dms2rad!* simpdms)), 'simparg);
deflist('((atan2 2) (hypot 2) (atan2d 2) (logb 2)),
'number!-of!-args);
flag('(acsc sind asind tand atand cotd acotd cscd acscd csch
acsch deg2rad rad2deg sgn),'odd);
flag('(cosd secd),'even);
flag('(sech),'nonzero);
symbolic smacro procedure mt!: x; cadr x;
symbolic smacro procedure ep!: x; cddr x;
symbolic procedure rdexp!* u; mkround
(if not atom x then exp!* x
else if x>!!maxarg then <<rndbfon(); exp!* bfloat x>>
else if x<-!!maxarg then 0.0 else exp x)
where x=convprec u;
symbolic procedure rdsqrt!* u;
mkround(if atom x then sqrt x else bfsqrt x)
where x=convprec u;
symbolic procedure rdexpt!*(u,v); mkround
(if not atom x then texpt!:any(x,y) else
((if z>!!maxarg then <<rndbfon(); texpt!:any(bfloat x,bfloat y)>>
else if z<-!!maxarg then 0.0 else expt(x,y))
where z=y*logfp bfloat x))
where x=convprec u,y=convprec v;
symbolic procedure rdlog!* u;
mkround(if atom x then log x else log!* x)
where x=convprec u;
symbolic procedure rdsgn!* u;
(if atom x then sgn x else sgn mt!: x) where x=round!* u;
symbolic procedure rdatan2!*(u,v);
mkround(if atom x then atan2(x,y) else atan2!*(x,y))
where x=convprec u,y=convprec v;
symbolic procedure rdatan2d!*(u,v);
mkround(if atom x then atan2d(x,y) else rad2deg!: atan2!*(x,y))
where x=convprec u,y=convprec v;
% symbolic procedure atan2!*(y,x);
% if mt!: x=0 then if (y := mt!: y)=0 then bfz!* else
% <<x := timbf(pi!/2!*(),bfhalf!*); if y<0 then minus!: x else x>>
% else <<(if mt!: x>0 then a
% else if mt!: y<0 then difbf(a,pi!/2!*())
% else plubf(a,pi!/2!*()))
% where a=atan!* divbf(y,x)>>;
symbolic procedure atan2!*(y,x);
if mt!: x=0 then if (y := mt!: y)=0 then bfz!* else
<<x := pi!/2!*(); if y<0 then minus!: x else x>>
else <<(if mt!: x>0 then a
else if mt!: y<0 then difbf(a,pi!*())
else plubf(a,pi!*()))
where a=atan!* divbf(y,x)>>;
symbolic procedure atan2d!*(y,x);
if mt!: x=0 then if (y := mt!: y)=0 then bfz!* else
<<x := timbf(!:180!*,bfhalf!*); if y<0 then minus!: x else x>>
else <<(if mt!: x>0 then a
else if mt!: y<0 then difbf(a,!:180!*) else plubf(a,!:180!*))
where a=rad2deg!: atan!* divbf(y,x)>>;
symbolic procedure rde!*; mkround if !*!*roundbf then e!*() else ee;
symbolic procedure rdpi!*; mkround if !*!*roundbf then pi!*() else pii;
symbolic procedure pi!/2!*; timbf(bfhalf!*,pi!*());
symbolic procedure deg2rad!* u;
mkround(if atom x then deg2rad x else deg2rad!: x)
where x=convprec u;
symbolic procedure rad2deg!* u;
mkround(if atom x then rad2deg x else rad2deg!: x)
where x=convprec u;
symbolic procedure deg2rad x; x*pi!/180;
symbolic procedure rad2deg x; x*!180!/pi;
symbolic procedure deg2rad!: x; divbf(timbf(x,pi!*()),!:180!*);
symbolic procedure rad2deg!: x; divbf(timbf(x,!:180!*),pi!*());
symbolic procedure rdsin!* u;
mkround (if atom x then sin x else sin!* x)
where x=convprec u;
symbolic procedure rdsind!* u;
mkround (if atom x then sind x else sin!* deg2rad!: x)
where x=convprec u;
symbolic procedure rdcos!* u;
mkround(if atom x then cos x else cos!* x)
where x=convprec u;
symbolic procedure rdcosd!* u;
mkround(if atom x then cosd x else cos!* deg2rad!: x)
where x=convprec u;
symbolic procedure rdtan!* u;
mkround(if atom x then tan x else tan!* x)
where x=convprec u;
symbolic procedure rdtand!* u;
mkround(if atom x then tand x else tan!* deg2rad!: x)
where x=convprec u;
symbolic procedure rdasin!* u;
mkround(if atom x then asin x else asin!* x)
where x=convprec u;
symbolic procedure rdasind!* u;
mkround(if atom x then asind x else rad2deg!: asin!* x)
where x=convprec u;
symbolic procedure rdacos!* u;
mkround(if atom x then acos x else acos!* x)
where x=convprec u;
symbolic procedure rdacosd!* u;
mkround(if atom x then acosd x else rad2deg!: acos!* x)
where x=convprec u;
symbolic procedure rdatan!* u;
mkround(if atom x then atan x else atan!* x)
where x=convprec u;
symbolic procedure rdatand!* u;
mkround(if atom x then atand x else rad2deg!: atan!* x)
where x=convprec u;
symbolic procedure rdsinh!* u;
mkround(if atom x then sinh x else sinh!* x)
where x=convprec u;
symbolic procedure rdcosh!* u;
mkround(if atom x then cosh x else cosh!* x)
where x=convprec u;
% these redefine functions that are in bfelem, and are faster.
symbolic procedure sinh!* x;
timbf(bfhalf!*,difbf(y,invbf y)) where y=exp!* x;
symbolic procedure cosh!* x;
timbf(bfhalf!*,plubf(y,invbf y)) where y=exp!* x;
% no bfelem functions after this point.
symbolic procedure rdsec!* u;
mkround(if atom x then sec x else invbf cos!* x)
where x=convprec u;
symbolic procedure rdsecd!* u;
mkround(if atom x then secd x else invbf cos!* deg2rad!: x)
where x=convprec u;
symbolic procedure rdcsc!* u;
mkround(if atom x then csc x else invbf sin!* x)
where x=convprec u;
symbolic procedure rdcscd!* u;
mkround(if atom x then cscd x else invbf sin!* deg2rad!: x)
where x=convprec u;
symbolic procedure rdcot!* u;
mkround(if atom x then cot x else tan!* difbf(pi!/2!*(),x))
where x=convprec u;
symbolic procedure rdcotd!* u;
mkround(if atom x then cotd x else tan!* difbf(pi!/2!*(),deg2rad!: x))
where x=convprec u;
symbolic procedure rdtanh!* u;
mkround(if atom x then tanh x else divbf(sinh!* x,cosh!* x))
where x=convprec u;
symbolic procedure rdcoth!* u;
mkround(if atom x then coth x else divbf(cosh!* x,sinh!* x))
where x=convprec u;
symbolic procedure rdsech!* u;
mkround(if atom x then sech x else invbf cosh!* x)
where x=convprec u;
symbolic procedure rdcsch!* u;
mkround(if atom x then csch x else invbf sinh!* x)
where x=convprec u;
symbolic procedure rdasinh!* u;
mkround(if atom x then asinh x else asinh!* x)
where x=convprec u;
symbolic procedure rdacosh!* u;
mkround(if atom x then acosh x else acosh!* x)
where x=convprec u;
symbolic procedure asinh!* x; begin scalar s;
if minusp!: x then x := minus!: x else s := t;
x := if leq!:(x,epsqrt!*) then x
else log!* plubf(x,
if lessp!:(x,bftwo!*) then bfsqrt plubf(timbf(x,x),bfone!*)
else if lessp!:(invbf x,epsqrt!*) then x
else timbf(x,bfsqrt plubf(bfone!*,divbf(bfone!*,timbf(x,x)))));
return if s then x else minus!: x end;
symbolic procedure acosh!* x;
if lessp!:(x,bfone!*) then terrlst(x,'acosh)
else log!* plubf(x,if leq!:(invbf x,epsqrt!*) then x
else timbf(x,bfsqrt difbf(bfone!*,divbf(bfone!*,timbf(x,x)))));
symbolic procedure rdacot!* u;
mkround(if atom x then acot x
else difbf(pi!/2!*(),atan!* x))
where x=convprec u;
symbolic procedure rdacotd!* u;
mkround(if atom x then acotd x
else rad2deg!: difbf(pi!/2!*(),atan!* x))
where x=convprec u;
symbolic procedure rdasec!* u; % not yet
mkround(if atom x then asec x else
difbf(pi!/2!*(),asin!* invbf x))
where x=convprec u;
symbolic procedure rdasecd!* u; % not yet
mkround(if atom x then asecd x else
rad2deg!: difbf(pi!/2!*(),asin!* invbf x))
where x=convprec u;
symbolic procedure rdacsc!* u;
mkround(if atom x then acsc x else asin!* invbf x)
where x=convprec u;
symbolic procedure rdacscd!* u;
mkround(if atom x then acscd x else rad2deg!: asin!* invbf x)
where x=convprec u;
symbolic procedure rdatanh!* u;
mkround(if atom x then atanh x else atanh!* x)
where x=convprec u;
symbolic procedure atanh!* x;
if not greaterp!:(bfone!*,abs!: x) then terrlst(x,'atanh)
else if leq!:(abs!: x,epsqrt!*) then x
else timbf(bfhalf!*,
log!* divbf(plubf(bfone!*,x),difbf(bfone!*,x)));
symbolic procedure rdacoth!* u;
mkround(if atom x then acoth x else atanh!* invbf x)
where x=convprec u;
symbolic procedure rdasech!* u; % not from here down
mkround(if atom x then asech x
else if leq!:(x,bfz!*) or greaterp!:(x,bfone!*)
then terrlst(x,'asech) else acosh!* invbf x)
where x=convprec u;
symbolic procedure rdacsch!* u;
mkround(if atom x then acsch x
else if mt!: x=0 then terrlst(x,'acsh) else asinh!* invbf x)
where x=convprec u;
symbolic procedure rdlogb!*(u,v);
mkround(if atom x then logb(x,b) else logb!*(x,b))
where x=convprec u,b=convprec v;
symbolic procedure rdlog10!* u;
mkround(if atom x then log10 x else logb!*(x,bften!*))
where x=convprec u;
symbolic procedure logb!* (x,b); %log x to base b;
begin scalar a,s;
a := greaterp!:(x,bfz!*);
s := not(leq!:(b,bfz!*) or equal!:(b,bfone!*));
if a and s then return divbf(log!* x,log!* b)
else terrlst((if a then list ('base,b)
else if s then list('arg,x) else list(x,b)),'logb) end;
symbolic procedure rdcbrt!* u;
mkround(if atom x then cbrt x else cbrt!* x)
where x=convprec u;
symbolic procedure cbrt!* x;
begin scalar s,l,g,u,nx,r; u := bfone!*;
if mt!: x=0 or equal!:(abs!: x,u) then return x
else if minusp!: x then x := minus!: x else s := t;
if lessp!:(x,u) then <<x := divbf(u,x); l := t>>
else if equal!:(x,u) then go to ret;
nx := '!:bf!: .
<<r := remainder(ep!:(nx := conv!:mt(x,3)),3);
if r=0 then (5+mt!: nx/179) . (ep!: nx/3)
else if r=1 or r=-2 then (10+mt!: nx/74) . ((ep!: nx-1)/3)
else (22+mt!: nx/35) . ((ep!: nx-2)/3)>>;
loop: r := nx;
nx := plubf(divbf(r,!:bf1!.5!*),
divbf(x,timbf(r,timbf(r,!:bf3!*))));
if g and leq!:(r,nx) then go to ret;
g := t; go to loop;
ret: if l then nx := divbf(u,nx);
return if s then nx else minus!: nx end;
symbolic procedure rdhypot!*(u,v);
mkround(if atom p then hypot(p,q) else hypot!*(p,q))
where p=convprec u,q=convprec v;
symbolic procedure hypot!*(p,q);
% Hypot(p,q)=sqrt(p*p+q*q) but avoids intermediate swell.
begin scalar r;
if minusp!: p then p := minus!: p; if minusp!: q then q := minus!: q;
if mt!: p=0 then return q else if mt!: q=0 then return p
else if lessp!:(p,q) then <<r := p; p := q; q := r>>;
r := divbf(q,p);
return if lessp!:(r,epsqrt!*) then p
else timbf(p,bfsqrt plubf(bfone!*,timbf(r,r))) end;
symbolic procedure simpdms l;
% Converts argument of form ({d,m,s}) to rd ((d m s)) if possible.
if cdr l or atom (l := car l) or not eqcar(l,'list)
or length l neq 4 then nil else
begin scalar fl; integer n;
l := for each a in cdr l collect
<<n := n+1;
if not (null(a := simprd list a) and (fl := t))
then a := car a;
if not fl and n>2 then i2rd!* round2a!* a else a>>;
if not fl then return list l end;
symbolic procedure round2a!* a; if atom a then a else round!* a;
symbolic procedure dms2rad!* u; deg2rad!* dms2deg!* u;
symbolic procedure dms2deg!* u;
mkround(if atom caddr l then dms2deg l else dms2deg!: l)
where l=list(round2a!* car u,round2a!* cadr u,round!* caddr u);
symbolic procedure dms2deg l; ((caddr l/60.0+cadr l)/60.0+car l);
symbolic procedure dms2deg!: l;
plubf(bfloat car l,divbf(plubf(bfloat cadr l,
divbf(bfloat caddr l,!:bf60!*)),!:bf60!*));
symbolic procedure rad2dms x; deg2dms rad2deg x;
symbolic procedure rad2dms!* u; deg2dms!* rad2deg!* u;
symbolic procedure deg2dms!* u;
mklist3!*(if atom x then deg2dms x else deg2dms!: x)
where x=round2a!* u;
symbolic procedure mklist3!* x; % floats seconds if not integer.
'list . list(car x,cadr x,
<<(if atom s and zerop(s-fix s) then fix s
else if not atom s and integerp!: s then conv!:bf2i s
else mk!*sq !*f2q mkround s) where s=caddr x>>);
symbolic procedure deg2dms x; % dms output in form list(d,m,s);
begin integer d,m;
m := fix(x := 60.0*(x-(d := fix2 x)));
return list(d,m,60.0*(x-m)) end;
symbolic procedure deg2dms!: x; % dms output in form list(d,m,s).
begin integer d,m;
d := conv!:bf2i x;
m := conv!:bf2i(x := timbf(!:bf60!*,difbf(x,bfloat d)));
return list(d,m,timbf(!:bf60!*,difbf(x,bfloat m))) end;
symbolic procedure rdnorm!* u; if rd!:minusp u then rd!:minus u else u;
symbolic procedure rdarg!* u;
if rd!:minusp u then rdpi!*() else rdzero!*();
% the following bfloat definitions are needed in addition to files
% smbflot and bfelem.red to support rdelem.
global '(!:bf!-1 !:bf!-2 !:bf!-0!.5 !:bf!-0!.25);
!:bf!-1 := '!:bf!: . (1 . 0); !:bf!-2 := '!:bf!: . (2 . 0);
symbolic smacro procedure preci!: nmbr;
% Returns precision of the bigfloat nmbr.
length explode abs mt!: nmbr$
symbolic smacro procedure make!:bf(mt,ep); '!:bf!: . (mt . ep);
symbolic procedure texpt!:(nmbr, k);
% This function calculates the Kth power of "n" up to the precision
% specified by !:PREC!:. %SK
% NMBR is a BIG-FLOAT representation of "n" and K an integer.
if not fixp k then bflerrmsg 'texpt!: % use texpt!:any in this case.
else if k=0 then !:bf!-1
else if k=1 then nmbr
else if k<0 then invbf texpt!:(nmbr,-k) %SK
else exptbf(nmbr,k,!:bf!-1); %SK
symbolic procedure quotient!:(n1, n2);
% This function calculates the integer quotient of "n1"
% and "n2", just as the "QUOTIENT" for integers does.
% **** For calculating the quotient up to a necessary
% **** precision, please use DIVIDE!:.
% N1 and N2 are BIG-FLOAT representations of "n1" and "n2".
begin integer e1, e2;
if (e1 := ep!: n1) = (e2 := ep!: n2) then return
make!:bf(mt!: n1 / mt!: n2, 0)
else if e1 > e2 then return
quotient!:(incprec!:(n1, e1 - e2) , n2)
else return
quotient!:(n1, incprec!:(n2, e2 - e1));
end$
symbolic procedure texpt!:any(x, y);
%modified by SK to use bfsqrt and exp!*, invbf and timbf.
% This function calculates the power x**y, where "x"
% and "y" are any numbers. The precision of
% the result is specified by !:PREC!:. % SK
% **** For a negative "x", this function returns
% **** -(-x)**y unless "y" is an integer.
% X is a BIG-FLOAT representation of "x".
% Y is either an integer, a floating-point number,
% or a BIG-FLOAT number, i.e., a BIG-FLOAT
% representation of "y".
if equal!:(x,e!*()) then exp!* bfloat y
else if fixp y then texpt!:(x, y)
else if integerp!: y then texpt!:(x,conv!:bf2i y)
else if not(bfp!: y or bfp!:(y := read!:num y))
then bflerrmsg 'texpt!:any % read!:num probably not necessary.
else if minusp!: y then invbf texpt!:any(x,minus!: y) %SK
else if equal!:(y,!:bf!-0!.5) then bfsqrt x %SK
else if equal!:(y,!:bf!-0!.25) then bfsqrt bfsqrt x %SK
else begin integer n; scalar xp, yp;
n := (if !:prec!: then !:prec!:
else max(preci!: x, preci!: y));
if minusp!: x then xp:=minus!: x else xp := x;
if integerp!: times!:(y,!:bf!-2) then
<< xp := incprec!:(xp, 1);
yp := texpt!:(xp, conv!:bf2i y);
yp := times!:(yp, sqrt!:(xp, n + 1));
yp := round!:mt(yp, n) >>
else
<< yp := timbf(y, log!:(xp, n + 1)); %SK
yp := exp!:(yp, n) >>;
return (if minusp!: x then minus!: yp else yp);
end$
symbolic procedure cut!:ep(nmbr, k);
% This function returns a given number "n" unchanged
% if its exponent is not less than K, else it
% cuts off its mantissa and returns an equivalent
% number of exponent K.
% **** CAUTION! No rounding is made.
% NMBR is a BIG-FLOAT representation of "n".
% K is an integer (positive or negative).
if bfp!: nmbr and fixp k then
if (k := k - ep!: nmbr) <= 0 then nmbr
else decprec!:(nmbr, k)
else bflerrmsg 'cut!:ep$
symbolic procedure integerp!: x;
% This function returns T if X is a BIG-FLOAT
% representing an integer, else it returns NIL.
% X is any LISP entity.
%JBM Critique: this is pretty slow. Couldn't we just check the
%JBM Critique: exponent in relation to the precision?
bfp!: x and
(ep!: x >= 0 or
equal!:(x,i2bf!: conv!:bf2i x));
symbolic procedure conv!:bf2i nmbr;
% This function converts a <BIG-FLOAT>, i.e., a BIG-FLOAT
% representation of "n", to an integer. The result
% is the integer part of "n".
% **** For getting the nearest integer to "n", please use
% **** the combination MT!:( CONV!:EP(NMBR,0)).
% NMBR is a BIG-FLOAT representation of the number "n".
if ep!:(nmbr := cut!:ep(nmbr, 0)) = 0 then mt!: nmbr
else mt!: nmbr * 10 ** ep!: nmbr;
symbolic procedure round!:ep(nmbr, k);
% This function rounds a number "n" and returns an
% equivalent number having the exponent K if
% the exponent of "n" is less than K, else
% it returns the given number unchanged.
% NMBR is a BIG-FLOAT representation of "n".
% K is an integer (positive or negative).
if bfp!: nmbr and fixp k then
if (k := k - 1 - ep!: nmbr) < 0 then nmbr
else if k = 0 then round!:last nmbr
else round!:last decprec!:(nmbr, k)
else bflerrmsg 'round!:ep$
endmodule;
module crelem; % Complex elementary functions for complex rounded.
fluid '(!*!*roundbf);
global '(!!flprec !!rdprec bfz!* bftwo!* bfone!* bfhalf!*);
symbolic smacro procedure tagrl x; '!:rd!: . cadr x;
symbolic smacro procedure tagim x; '!:rd!: . cddr x;
symbolic procedure rdone!*; if !*!*roundbf then bfone!* else 1.0;
symbolic procedure rdtwo!*; if !*!*roundbf then bftwo!* else 2.0;
symbolic procedure rdhalf!*; if !*!*roundbf then bfhalf!* else 0.5;
symbolic procedure rdzero!*; if !*!*roundbf then bfz!* else 0.0;
symbolic procedure crhalf!*; mkcr(rdhalf!*(),rdzero!*());
symbolic procedure cri!*; mkcr(rdzero!*(),rdone!*());
symbolic procedure cri!/2; mkcr(rdzero!*(),rdhalf!*());
symbolic procedure cr2i!*; mkcr(rdzero!*(),rdtwo!*());
deflist('((expt crexpt!*) (sin crsin!*) (cos crcos!*) (tan crtan!*)
(asin crasin!*) (acos cracos!*) (atan cratan!*)
(cot crcot!*) (acot cracot!*) (sec crsec!*) (asec crasec!*)
(sind crsind!*) (cosd crcosd!*) (tand crtand!*)
(asind crasind!*) (acosd cracosd!*) (atand cratand!*)
(cotd crcotd!*) (acotd cracotd!*) (secd crsecd!*)
(asecd crasecd!*) (arg crarg!*) (argd crargd!*)
(csc crcsc!*) (acsc cracsc!*) (sinh crsinh!*) (cosh crcosh!*)
(cscd crcscd!*) (acscd cracscd!*)
(asinh crasinh!*) (acosh cracosh!*) (tanh crtanh!*)
(coth crcoth!*) (atanh cratanh!*) (acoth cracoth!*)
(sech crsech!*) (csch crcsch!*) (asech crasech!*)
(acsch cracsch!*)
(sqrt crsqrt!*) (norm crnorm!*) (arg crarg!*) (log crlog!*)
(exp crexp!*) (logb crlogb!*) (e cre!*) (pi crpi!*)),'!:cr!:);
symbolic procedure cre!*; mkcr(rde!*(),rdzero!*());
symbolic procedure crpi!*; mkcr(rdpi!*(),rdzero!*());
symbolic procedure crexpt!*(u,v);
if cr!:zerop(cr!:differ(v,crhalf!*())) then crsqrt!* u
else crexp!* cr!:times(v,crlog!* u);
symbolic procedure crnorm!* u; rdhypot!*(tagrl u,tagim u);
symbolic procedure crarg!* u; rdatan2!*(tagim u,tagrl u);
symbolic procedure crargd!* u; rdatan2d!*(tagim u,tagrl u);
symbolic procedure crsqrt!* u; gf2cr!: gfsqrt crprcd u;
symbolic procedure crr2d!* u; mkcr(rad2deg!* tagrl u,rad2deg!* tagim u);
symbolic procedure crd2r!* u; mkcr(deg2rad!* tagrl u,deg2rad!* tagim u);
symbolic procedure crsin!* u;
mkcr(rd!:times(rdsin!* rl,rdcosh!* im),
rd!:times(rdcos!* rl,rdsinh!* im))
where rl=tagrl u,im=tagim u;
symbolic procedure crsind!* u;
mkcr(rd!:times(rdsind!* rl,rdcosh!* deg2rad!* im),
rd!:times(rdcos!* rl,rdsinh!* deg2rad!* im))
where rl=tagrl u,im=tagim u;
symbolic procedure crcos!* u;
mkcr(rd!:times(rdcos!* rl,rdcosh!* im),
rd!:minus rd!:times(rdsin!* rl,rdsinh!* im))
where rl=tagrl u,im=tagim u;
symbolic procedure crcosd!* u;
mkcr(rd!:times(rdcosd!* rl,rdcosh!* deg2rad!* im),
rd!:minus rd!:times(rdsind!* rl,rdsinh!* deg2rad!* im))
where rl=tagrl u,im=tagim u;
symbolic procedure crtan!* u;
cr!:times(cri!*(),cr!:quotient(cr!:differ(y,x),cr!:plus(y,x)))
where x=crexp!*(cr!:times(cr2i!*(),u)),y=i2cr!* 1;
symbolic procedure crtand!* u;
cr!:times(cri!*(),cr!:quotient(cr!:differ(y,x),cr!:plus(y,x)))
where x=crexp!*(cr!:times(cr2i!*(),crd2r!* u)),y=i2cr!* 1;
symbolic procedure crcot!* u;
cr!:times(cri!*(),cr!:quotient(cr!:plus(x,y),cr!:differ(x,y)))
where x=crexp!*(cr!:times(cr2i!*(),u)),y=i2cr!* 1;
symbolic procedure crcotd!* u;
cr!:times(cri!*(),cr!:quotient(cr!:plus(x,y),cr!:differ(x,y)))
where x=crexp!*(cr!:times(cr2i!*(),crd2r!* u)),y=i2cr!* 1;
symbolic procedure crlog!* u;
mkcr(rdlog!* crnorm!* u,crarg!* u);
symbolic procedure crlogb!*(u,b); cr!:quotient(crlog!* u,crlog!* b);
symbolic procedure timesi!* u; cr!:times(cri!*(),u);
symbolic procedure crasin!* u; cr!:minus timesi!* crasinh!* timesi!* u;
symbolic procedure crasind!* u;
crr2d!* cr!:minus timesi!* crasinh!* timesi!* u;
symbolic procedure cracos!* u;
cr!:plus(cr!:times(i2cr!* 0.5,crpi!*()),
timesi!* crasinh!* timesi!* u);
symbolic procedure cracosd!* u;
crr2d!* cr!:plus(cr!:times(i2cr!* 0.5,crpi!*()),
timesi!* crasinh!* timesi!* u);
symbolic procedure cratan!* u;
cr!:times(cri!/2(),crlog!* cr!:quotient(
cr!:plus(cri!*(),u),cr!:differ(cri!*(),u)));
symbolic procedure cratand!* u;
crr2d!* cr!:times(cri!/2(),crlog!* cr!:quotient(
cr!:plus(cri!*(),u),cr!:differ(cri!*(),u)));
symbolic procedure cracot!* u;
cr!:times(cri!/2(),crlog!* cr!:quotient(
cr!:differ(u,cri!*()),cr!:plus(cri!*(),u)));
symbolic procedure cracotd!* u;
crr2d!* cr!:times(cri!/2(),crlog!* cr!:quotient(
cr!:differ(u,cri!*()),cr!:plus(cri!*(),u)));
symbolic procedure crsec!* u; cr!:quotient(i2cr!* 1,crcos!* u);
symbolic procedure crsecd!* u;
cr!:quotient(i2cr!* 1,crcos!* crd2r!* u);
symbolic procedure crcsc!* u; cr!:quotient(i2cr!* 1,crsin!* u);
symbolic procedure crcscd!* u;
cr!:quotient(i2cr!* 1,crsin!* crd2r!* u);
symbolic procedure crasec!* u; cracos!* cr!:quotient(i2cr!* 1,u);
symbolic procedure crasecd!* u;
crr2d!* cracos!* cr!:quotient(i2cr!* 1,u);
symbolic procedure cracsc!* u; crasin!* cr!:quotient(i2cr!* 1,u);
symbolic procedure cracscd!* u;
crr2d!* crasin!* cr!:quotient(i2cr!* 1,u);
symbolic procedure crsinh!* u;
cr!:times(i2cr!* 0.5,cr!:differ(y,cr!:quotient(i2cr!* 1,y)))
where y=crexp!* u;
symbolic procedure crcosh!* u;
cr!:times(i2cr!* 0.5,cr!:plus(y,cr!:quotient(i2cr!* 1,y)))
where y=crexp!* u;
symbolic procedure crtanh!* u;
cr!:quotient(cr!:differ(x,y),cr!:plus(x,y))
where x=crexp!*(cr!:times(i2cr!* 2,u)),y=i2cr!* 1;
symbolic procedure crcoth!* u;
cr!:quotient(cr!:plus(x,y),cr!:differ(x,y))
where x=crexp!*(cr!:times(i2cr!* 2,u)),y=i2cr!* 1;
symbolic procedure crsech!* u;
cr!:quotient(i2cr!* 2,cr!:plus(y,cr!:quotient(i2cr!* 1,y)))
where y=crexp!* u;
symbolic procedure crcsch!* u;
cr!:quotient(i2cr!* 2,cr!:differ(y,cr!:quotient(i2cr!* 1,y)))
where y=crexp!* u;
symbolic procedure crasinh!* u;
crlog!* cr!:plus(u,
if bflessp(round!* crnorm!* u,rdtwo!*())
then crsqrt!* cr!:plus(i2cr!* 1,s)
else cr!:times(u,
crsqrt!* cr!:plus(i2cr!* 1,cr!:quotient(i2cr!* 1,s))))
where s=cr!:times(u,u);
symbolic procedure cracosh!* u;
crlog!* cr!:plus(u,crsqrt!* cr!:differ(cr!:times(u,u),i2cr!* 1));
symbolic procedure cratanh!* u;
cr!:times(i2cr!* 0.5,crlog!* cr!:quotient(cr!:plus(i2cr!* 1,u),
cr!:differ(i2cr!* 1,u)));
symbolic procedure cracoth!* u;
cr!:times(i2cr!* 0.5,crlog!* cr!:quotient(cr!:plus(i2cr!* 1,u),
cr!:differ(u,i2cr!* 1)));
symbolic procedure crasech!* u; cracosh!* cr!:quotient(i2cr!* 1,u);
symbolic procedure cracsch!* u; crasinh!* cr!:quotient(i2cr!* 1,u);
symbolic procedure crexp!* u;
<<u := tagim u; mkcr(rd!:times(r,rdcos!* u),rd!:times(r,rdsin!* u))>>
where r=rdexp!* tagrl u;
endmodule;
module bfelem; % Bigfloat elementary constants and functions.
% Author: T. Sasaki, 1979.
% Modifications by: Anthony C. Hearn, Jed B. Marti, Stanley L. Kameny.
symbolic;
fluid '(!:prec!:);
global '(bfsaveprec!*);
global '(!:bf!-pi %PI to 20 digits.
!:bf!-0 %0.0
!:bf!-1 %1.0
!:bf!-e %E to 20 digits
!:bf!-0!.5 %0.5
!:bf!-0!.25 %0.25
!:bf!-0!.1 %0.1
!:bf!-1!.72 %1.72
!:bf!-0!.42 %0.42
!:bf!-0!.72 %0.72
);
% *** Tables for Elementary Function and Constant Values ***
symbolic smacro procedure make!:bf(mt,ep); '!:bf!: . (mt . ep)$
symbolic smacro procedure mt!: u; cadr u;
symbolic smacro procedure preci!: nmbr;
% Returns precision of the bigfloat nmbr.
length explode abs mt!: nmbr$
symbolic smacro procedure ep!: u; cddr u;
symbolic procedure allfixp l;
% Returns T if all of L are FIXP.
null l or fixp car l and allfixp cdr l;
symbolic procedure read!:lnum(l);
% This function reads a long number "n" represented by a list in a way
% described below, and constructs a BIG-FLOAT representation of "n".
% L is a list of integers, the first element of which gives the order of
% "n" and all the next elements when concatenated give the mantissa of
% "n".
% **** ORDER(n)=k if 10**k <= ABS(n) < 10**(k+1).
% **** Except for the first element, all integers in L
% **** should not begin with "0" because some
% **** systems suppress leading zeros.
% JBM: Fix some kludgy coding here.
% JBM: Add BFSAVEPREC!* precision saver.
if not allfixp l then bflerrmsg 'read!:lnum
else begin scalar mt, ep, k, sign, u, v, dcnt;
mt := dcnt := 0; %JBM
% ep := car(u := l) + 1; %JBM
u := l;
ep := add1 car u;
sign := if minusp cadr l then -1 else 1; %JBM
while u:=cdr u do
<< k := length explode(v := abs car u); %JBM
% k := 0; %JBM
% while v do << k := k + 1; v := cdr v >>; %JBM
mt := mt * 10**k + v; %JBM
ep := ep - k;
dcnt := dcnt + k; % JBM
if bfsaveprec!* and dcnt > bfsaveprec!* then %JBM
u := '(nil) >>; %JBM
return make!:bf(sign * mt, ep);
end$
deflist('((exp exp!*) (expt bfexpt!:) (log log!*) (sin sin!*)
(cos cos!*) (tan tan!*) (asin asin!*) (acos acos!*)
(atan atan!*) (sqrt sqrt!*) (sinh sinh!*) (cosh cosh!*)
(e e!*) (pi pi!*)),
'!:bf!:);
symbolic procedure bfexpt!:(u,v);
% Calculates u**v, including case u<0.
if minusp!: u
then multd(texpt!:any(minus!: u,v),
!*q2f if null numr simp list('difference,v,
'(quotient 1 2))
then simp 'i
else mksq(list('expt,'(minus 1),v),1))
else texpt!:any(u,v);
symbolic procedure exp!* u; exp!:(u,!:prec!:);
symbolic procedure log!* u; log!:(u,!:prec!:);
symbolic procedure sin!* u; sin!:(u,!:prec!:);
symbolic procedure cos!* u; cos!:(u,!:prec!:);
symbolic procedure tan!* u; tan!:(u,!:prec!:);
symbolic procedure asin!* u; asin!:(u,!:prec!:);
symbolic procedure acos!* u; acos!:(u,!:prec!:);
symbolic procedure atan!* u; atan!:(u,!:prec!:);
symbolic procedure sqrt!* u; sqrt!:(u,!:prec!:);
symbolic procedure pi!*;
if !:prec!:>1000 then !:bigpi !:prec!: else !:pi !:prec!:;
symbolic procedure e!*; !:e !:prec!:;
%*************************************************************
%** **
%** 3-1. Elementary CONSTANTS. **
%** **
%*************************************************************
symbolic procedure !:pi k;
% This function calculates the value of the circular
% constant "PI", with the precision K, by
% using Machin's well known identity:
% PI = 16*atan(1/5) - 4*atan(1/239).
% Calculation is performed mainly on integers.
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg '!:pi
else if k <= 20 then
% round!:mt(make!:bf(314159265358979323846, -20), k)
round!:mt(!:bf!-pi, k) %JBM
else
begin integer k3,s,ss,m,n,x; scalar u;
u := get!:const('!:pi, k);
if u neq "NOT FOUND" then return u;
ss := n := 10 ** (k3 := k + 3) / 5;
x := -5 ** 2;
m := 1;
while n neq 0 do <<n := n/x; ss := ss + n/(m := m + 2)>>;
s := n := 10 ** k3 / 239;
x := -239 ** 2;
m := 1;
while n neq 0 do << n := n / x; s := s + n / (m := m + 2) >>;
u := round!:mt(make!:bf(16 * ss - 4 * s, - k3), k);
save!:const('!:pi, u);
return u;
end$
symbolic procedure !:bigpi k;
% This function calculates the value of the circular
% constant "PI", with the precision K, by the
% arithmetic-geometric mean method. (See,
% R. Brent, JACM Vol.23, #2, pp.242-251(1976).)
% K is a positive integer.
% **** This function should be used only when you
% **** need "PI" of precision higher than 1000.
if not fixp k or k <= 0 then bflerrmsg '!:bigpi
else begin integer k2, n; scalar dcut, half, x, y, u, v;
u := get!:const('!:pi, k);
if u neq "NOT FOUND" then return u;
k2 := k + 2;
half := !:bf!-0!.5; %JBM
dcut := make!:bf(10, - k2);
x := i2bf!:(n := 1);
y := divide!:(x, !:sqrt2 k2, k2);
u := !:bf!-0!.25; %JBM
while greaterp!:(abs!: difference!:(x, y), dcut) do
<< v := x;
x := times!:(plus!:(x, y), half);
y := sqrt!:(cut!:ep(times!:(y, v), - k2), k2);
v := difference!:(x, v);
v := times!:(times!:(v, v),i2bf!: n);
u := difference!:(u, cut!:ep(v, - k2));
n := 2*n>>;
v := cut!:mt(texpt!:(plus!:(x, y), 2), k2);
u := divide!:(v, times!:(i2bf!: 4, u), k);
save!:const('!:pi, u);
return u;
end$
symbolic procedure !:e k;
% This function calculates the value of "e", the base
% of the natural logarithm, with the precision K,
% by summing the Taylor series for exp(x=1).
% Calculation is performed mainly on integers.
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg '!:e
else if k <= 20 then
% round!:mt(make!:bf(271828182845904523536, -20), k)
round!:mt(!:bf!-e, k) %JBM
else begin integer k2, ans, m, n; scalar u;
u := get!:const('!:e, k);
if u neq "NOT FOUND" then return u;
k2 := k + 2;
m := 1;
n := 10 ** k2;
ans := 0;
while n neq 0 do ans := ans + (n := n / (m := m + 1));
ans := ans + 2 * 10 ** k2;
u := round!:mt(make!:bf(ans, - k2), k);
save!:const('!:e2, u);
return u;
end$
symbolic procedure !:e01(k);
% This function calculates exp(0.1), the value of the
% exponential function at the point 0.1, with
% the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:e01, k);
if u neq "NOT FOUND" then return u;
u := exp!:(!:bf!-0!.1, k); %JBM
save!:const('!:e01, u);
return u;
end$
symbolic procedure !:log2 k;
% This function calculates log(2), the natural
% logarithm of 2, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:log2, k);
if u neq "NOT FOUND" then return u;
u := log!:(i2bf!: 2, k);
save!:const('!:log2, u);
return u;
end$
symbolic procedure !:log3 k;
% This function calculates log(3), the natural
% logarithm of 3, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:log3, k);
if u neq "NOT FOUND" then return u;
u := log!:(i2bf!: 3, k);
save!:const('!:log3, u);
return u;
end$
symbolic procedure !:log5 k;
% This function calculates log(5), the natural
% logarithm of 5, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:log5, k);
if u neq "NOT FOUND" then return u;
u := log!:(i2bf!: 5, k);
save!:const('!:log5, u);
return u;
end$
symbolic procedure !:log10 k;
% This function calculates log(10), the natural
% logarithm of 10, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:log10, k);
if u neq "NOT FOUND" then return u;
u := log!:(i2bf!: 10, k);
save!:const('!:log10, u);
return u;
end$
symbolic procedure !:logpi k;
% This function calculates log(PI), the natural
% logarithm of "PI", with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:logpi, k);
if u neq "NOT FOUND" then return u;
u := log!:(!:pi(k + 2), k);
save!:const('!:logpi, u);
return u
end$
symbolic procedure !:sqrt2(k);
% This function calculates SQRT(2), the square root
% of 2, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:sqrt2, k);
if u neq "NOT FOUND" then return u;
u := sqrt!:(i2bf!: 2, k);
save!:const('!:sqrt2, u);
return u;
end$
symbolic procedure !:sqrt3(k);
% This function calculates SQRT(3), the square root
% of 3, with the precision K.
% K is a positive integer.
begin scalar u;
u:=get!:const('!:sqrt3, k);
if u neq "NOT FOUND" then return u;
u := sqrt!:(i2bf!: 3, k);
save!:const('!:sqrt3, u);
return u;
end$
symbolic procedure !:sqrt5 k;
% This function calculates SQRT(5), the square root
% of 5, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:sqrt5, k);
if u neq "NOT FOUND" then return u;
u := sqrt!:(i2bf!: 5, k);
save!:const('!:sqrt5, u);
return u;
end$
symbolic procedure !:sqrt10 k;
% This function calculates SQRT(10), the square root
% of 10, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:sqrt10, k);
if u neq "NOT FOUND" then return u;
u := sqrt!:(i2bf!: 10, k);
save!:const('!:sqrt10, u);
return u;
end$
symbolic procedure !:sqrtpi k;
% This function calculates SQRT(PI), the square root
% of "PI", with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:sqrtpi, k);
if u neq "NOT FOUND" then return u;
u := sqrt!:(!:pi(k + 2), k);
save!:const('!:sqrtpi, u);
return u;
end$
symbolic procedure !:sqrte k;
% This function calculates SQRT(e), the square root
% of "e", with the precision K.
% K is a positive integer.
begin scalar u;
u:=get!:const('!:sqrte, k);
if u neq "NOT FOUND" then return u;
u := sqrt!:(!:e(k + 2), k);
save!:const('!:sqrte, u);
return u;
end$
symbolic procedure !:cbrt2 k;
% This function calculates CBRT(2), the cube root
% of 2, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:cbrt2, k);
if u neq "NOT FOUND" then return u;
u := cbrt!:(i2bf!: 2, k);
save!:const('!:cbrt2, u);
return u;
end$
symbolic procedure !:cbrt3 k;
% This function calculates CBRT(3), the cube root
% of 3, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:cbrt3, k);
if u neq "NOT FOUND" then return u;
u := cbrt!:(i2bf!: 3, k);
save!:const('!:cbrt3, u);
return u;
end$
symbolic procedure !:cbrt5 k;
% This function calculates CBRT(5), the cube root
% of 5, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:cbrt5, k);
if u neq "NOT FOUND" then return u;
u := cbrt!:(i2bf!: 5, k);
save!:const('!:cbrt5, u);
return u;
end$
symbolic procedure !:cbrt10 k;
% This function calculates CBRT(10), the cube root
% of 10, with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:cbrt10, k);
if u neq "NOT FOUND" then return u;
u := cbrt!:(i2bf!: 10, k);
save!:const('!:cbrt10, u);
return u;
end$
symbolic procedure !:cbrtpi k;
% This function calculates CBRT(PI), the cube root
% of "PI", with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:cbrtpi, k);
if u neq "NOT FOUND" then return u;
u := cbrt!:(!:pi(k + 2), k);
save!:const('!:cbrtpi, u);
return u;
end$
symbolic procedure !:cbrte k;
% This function calculates CBRT(e), the cube root
% of "e", with the precision K.
% K is a positive integer.
begin scalar u;
u := get!:const('!:cbrte, k);
if u neq "NOT FOUND" then return u;
u := cbrt!:(!:e(k + 2), k);
save!:const('!:cbrte, u);
return u;
end$
%*************************************************************
%** **
%** 3-2. Routines for saving CONSTANTS. **
%** **
%*************************************************************
symbolic procedure get!:const(cnst, k);
% This function returns the value of constant CNST
% of the precision K, if it was calculated
% previously with, at least, the precision K,
% else it returns "NOT FOUND".
% CNST is the name of the constant (to be quoted).
% K is a positive integer.
if atom cnst and fixp k and k > 0 then
begin scalar u;
u := get(cnst, 'save!:c);
if null u or car u < k then return "NOT FOUND"
else if car u = k then return cdr u
else return round!:mt(cdr u, k);
end
else bflerrmsg 'get!:const$
symbolic procedure save!:const(cnst, nmbr);
% This function saves the value of constant CNST
% for the later use.
% CNST is the name of the constant (to be quoted).
% NMBR is a BIG-FLOAT representation of the value.
if atom cnst and bfp!: nmbr then
put(cnst, 'save!:c, preci!: nmbr . nmbr)
else bflerrmsg 'save!:const$
symbolic procedure set!:const(cnst, l);
% This function sets the value of constant CNST.
% CNST is the name of the constant (to be quoted).
% L is a list of integers, which represents the
% value of the constant in the way described
% in the function READ!:LNUM.
save!:const(cnst, read!:lnum l)$
% Setting the constants.
set!:const( '!:pi , '( 0 3141 59265 35897 93238 46264
33832 79502 88419 71693 99375 105820 9749 44592 30781
64062 86208 99862 80348 25342 11706 79821 48086 51328
23066 47093 84460 95505 82231 72535 94081 28481 1174
5028410 2701 93852 11055 59644 62294 89549 30381 96442
88109 8) )$
set!:const( '!:e , '( 0 2718 28182 84590 45235 36028
74713 52662 49775 72470 93699 95957 49669 67627 72407
66303 53547 59457 13821 78525 16642 74274 66391 93200
30599 21817 41359 66290 43572 90033 42952 60595 63073
81323 28627 943490 7632 33829 88075 31952 510190 1157
38341 9) )$
set!:const( '!:e01 , '( 0 1105 17091 80756 47624 81170
78264 90246 66822 45471 94737 51871 87928 63289 44096
79667 47654 30298 91433 18970 74865 36329 2) )$
set!:const( '!:log2 , '(-1 6931 47180 55994 53094 17232
12145 81765 68075 50013 43602 55254 1206 800094 93393
62196 96947 15605 86332 69964 18687 54200 2) )$
set!:const( '!:log3 , '( 0 1098 61228 866810 9691 39524
52369 22525 70464 74905 57822 74945 17346 94333 63749
42932 18608 96687 36157 54813 73208 87879 7) )$
set!:const( '!:log5 , '( 0 1609 43791 2434100 374 60075
93332 26187 63952 56013 54268 51772 19126 47891 47417
898770 7657 764630 1338 78093 179610 7999 7) )$
set!:const( '!:log10 , '( 0 2302 58509 29940 456840 1799
14546 84364 20760 11014 88628 77297 60333 27900 96757
26096 77352 48023 599720 5089 59829 83419 7) )$
set!:const( '!:logpi , '( 0 1144 72988 5849400 174 14342
73513 53058 71164 72948 12915 31157 15136 23071 47213
77698 848260 7978 36232 70275 48970 77020 1) )$
set!:const( '!:sqrt2 , '( 0 1414 21356 23730 95048 80168
872420 96980 7856 96718 75376 94807 31766 79737 99073
24784 621070 38850 3875 34327 64157 27350 1) )$
set!:const( '!:sqrt3 , '( 0 17320 5080 75688 77293 52744
634150 5872 36694 28052 53810 38062 805580 6979 45193
301690 88000 3708 11461 86757 24857 56756 3) )$
set!:const( '!:sqrt5 , '( 0 22360 6797 74997 89696 40917
36687 31276 235440 6183 59611 52572 42708 97245 4105
209256 37804 89941 441440 8378 78227 49695 1) )$
set!:const( '!:sqrt10, '( 0 3162 277660 1683 79331 99889
35444 32718 53371 95551 39325 21682 685750 4852 79259
44386 39238 22134 424810 8379 30029 51873 47))$
set!:const( '!:sqrtpi, '( 0 1772 453850 9055 16027 29816
74833 41145 18279 75494 56122 38712 821380 7789 85291
12845 91032 18137 49506 56738 54466 54162 3) )$
set!:const( '!:sqrte , '( 0 1648 721270 7001 28146 8486
507878 14163 57165 3776100 710 14801 15750 79311 64066
10211 94215 60863 27765 20056 36664 30028 7) )$
set!:const( '!:cbrt2 , '( 0 1259 92104 98948 73164 7672
106072 78228 350570 2514 64701 5079800 819 75112 15529
96765 13959 48372 93965 62436 25509 41543 1) )$
set!:const( '!:cbrt3 , '( 0 1442 249570 30740 8382 32163
83107 80109 58839 18692 53499 35057 75464 16194 54168
75968 29997 33985 47554 79705 64525 66868 4) )$
set!:const( '!:cbrt5 , '( 0 1709 97594 66766 96989 35310
88725 43860 10986 80551 105430 5492 43828 61707 44429
592050 4173 21625 71870 10020 18900 220450 ) )$
set!:const( '!:cbrt10, '( 0 2154 4346900 318 83721 75929
35665 19350 49525 93449 42192 10858 24892 35506 34641
11066 48340 80018 544150 3543 24327 61012 6) )$
set!:const( '!:cbrtpi, '( 0 1464 59188 75615 232630 2014
25272 63790 39173 85968 55627 93717 43572 55937 13839
36497 98286 26614 56820 67820 353820 89750 ) )$
set!:const( '!:cbrte , '( 0 1395 61242 50860 89528 62812
531960 2586 83759 79065 15199 40698 26175 167060 3173
90156 45951 84696 97888 17295 83022 41352 1) )$
%*************************************************************
%** **
%** 4-1. Elementary FUNCTIONS. **
%** **
%*************************************************************
symbolic procedure sqrt!:(x, k);
% This function calculates SQRT(x), the square root
% of "x", with the precision K, by Newton's
% iteration method.
% X is a BIG-FLOAT representation of "x", x >= 0.
% K is a positive integer.
if minusp!: x or not fixp k or k <= 0 then bflerrmsg 'sqrt!:
else if bfzerop!: x then i2bf!: 0
else begin integer k2,ncut,nfig; scalar dcut,half,dy,y,y0,u;
k2 := k + 2;
ncut := k2 - (order!: x + 1) / 2;
half := !:bf!-0!.5; %JBM
dcut := make!:bf(10, - ncut);
dy := make!:bf(20, - ncut);
y0 := conv!:mt(x, 2);
if remainder(ep!: y0, 2) = 0 then
y0 := make!:bf(3 + 2 * mt!: y0 / 25, ep!: y0 / 2)
else y0 := make!:bf(10 + 2 * mt!: y0 / 9, (ep!: y0 - 1) / 2);
nfig := 1;
while nfig < k2 or greaterp!:(abs!: dy, dcut) do
<< if (nfig := 2 * nfig) > k2 then nfig := k2;
u := divide!:(x, y0, nfig);
y := times!:(plus!:(y0, u), half);
dy := difference!:(y, y0);
y0 := y >>;
return round!:mt(y, k);
end$
symbolic procedure cbrt!:(x, k);
% This function calculates CBRT(x), the cube root
% of "x", with the precision K, by Newton's
% iteration method.
% X is a BIG-FLOAT representation of any real "x".
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg 'cbrt!:
else if bfzerop!: x then i2bf!: 0
else if minusp!: x then minus!: cbrt!:(minus!: x, k)
else begin integer k2, ncut, nfig, j; scalar dcut, thre, dy, y, u;
k2 := k + 2;
ncut := k2 - (order!: x + 2) / 3;
thre := i2bf!: 3;
dcut := make!:bf(10, - ncut);
dy := make!:bf(20, - ncut);
y := conv!:mt(x, 3);
if (j := remainder(ep!: y, 3)) = 0 then
y := make!:bf(5 + mt!: y / 167, ep!: y / 3)
else if j = 1 or j = -2 then
y := make!:bf(10 + mt!: y / 75, (ep!: y - 1) / 3)
else y := make!:bf(22 + 2 * mt!: y / 75, (ep!: y - 2) / 3);
nfig := 1;
while nfig < k2 or greaterp!:(abs!: dy, dcut) do
<< if (nfig := 2 * nfig) > k2 then nfig := k2;
u := cut!:mt(times!:(y, y), nfig);
u := divide!:(x, u, nfig);
j :=order!:(u := difference!:(u, y)) + ncut - k2;
dy := divide!:(u, thre, max(1, nfig + j));
y := plus!:(y, dy) >>;
return round!:mt(y, k);
end$
symbolic procedure exp!:(x, k);
% This function calculates exp(x), the value of
% the exponential function at the point "x",
% with the precision K, by summing terms of
% the Taylor series for exp(z), 0 < z < 1.
% X is a BIG-FLOAT representation of any real "x".
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg 'exp!:
else if bfzerop!: x then i2bf!: 1
else begin integer k2, m; scalar one, q, r, y, yq, yr, save!:p;
k2 := k + 2;
one := i2bf!: 1;
q := i2bf!:(m := conv!:bf2i(y := abs!: x));
r := difference!:(y, q);
if bfzerop!: q then yq := one
else << save!:p := !:prec!:;
!:prec!: := k2;
yq := texpt!:(!:e k2, m);
!:prec!: := save!:p >>;
if bfzerop!: r then yr:=one
else begin integer j, n; scalar dcut, fctrial, ri, tm;
dcut := make!:bf(10, - k2);
yr := ri := tm := one;
m := 1;
j := 0;
while greaterp!:(tm, dcut) do
<< fctrial := i2bf!:(m := m * (j := j + 1));
ri := cut!:ep(times!:(ri, r), - k2);
n := max(1, k2 - order!: fctrial + order!: ri);
tm := divide!:(ri, fctrial, n);
yr := plus!:(yr,tm);
if remainder(j,10)=0 then yr := cut!:ep(yr, - k2) >>;
end;
y := cut!:mt(times!:(yq, yr), k + 1);
return (if minusp!: x then divide!:(one, y, k)
else round!:last y);
end$
symbolic procedure log!:(x, k);
% This function calculates log(x), the value of the
% logarithmic function at the point "x", with
% the precision K, by summing terms of the
% Taylor series for log(1+z), 0 < z < 0.10518.
% X is a BIG-FLOAT representation of "x", x > 0.
% K is a positive integer.
if minusp!: x or bfzerop!: x or
not fixp k or k <= 0 then bflerrmsg 'log!:
else if equal!:(x,i2bf!: 1) then i2bf!: 0
else begin integer k2,m; scalar eee,es,one,sign,l,y,z,save!:p;
k2 := k + 2;
one := i2bf!: 1;
eee := !:e k2;
es := !:e01 k2;
if greaterp!:(x, one) then << sign := one; y := x >>
else << sign := minus!: one; y := divide!:(one, x, k2) >>;
if lessp!:(y, eee) then << m := 0; z := y >>
else << if (m := (order!: y * 23) / 10) = 0 then z := y
else << save!:p := !:prec!:;
!:prec!: := k2;
z := divide!:(y, texpt!:(eee, m), k2);
!:prec!: := save!:p >>;
while greaterp!:(z, eee) do
<< m := m+1; z := divide!:(z, eee, k2) >> >>;
l := i2bf!: m;
y := !:bf!-0!.1;
while greaterp!:(z, es) do
<< l := plus!:(l, y); z := divide!:(z, es, k2) >>;
z := difference!:(z, one);
begin integer n; scalar dcut, tm, zi;
y := tm := zi := z;
z := minus!: z;
dcut := make!:bf(10, - k2);
m := 1;
while greaterp!:(abs!: tm, dcut) do
<< zi := cut!:ep(times!:(zi, z), - k2);
n := max(1, k2 + order!: zi);
tm := divide!:(zi,i2bf!:(m := m + 1), n);
y := plus!:(y, tm);
if zerop remainder(m,10) then y := cut!:ep(y,-k2)>>;
end;
y := plus!:(y, l);
return round!:mt(times!:(sign, y), k);
end$
symbolic procedure sin!:(x, k);
% This function calculates sin(x), the value of
% the sine function at the point "x", with
% the precision K, by summing terms of the
% Taylor series for sin(z), 0 < z < PI/4.
% X is a BIG-FLOAT representation of any real "x".
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg 'sin!:
else if bfzerop!: x then i2bf!: 0
else if minusp!: x then minus!: sin!:(minus!: x, k)
else begin integer k2, m; scalar pi4, sign, q, r, y;
k2 := k + 2;
m := preci!: x;
pi4 := times!:(!:pi(k2 + m), !:bf!-0!.25);
if lessp!:(x, pi4) then << m := 0; r := x >>
else << m := conv!:bf2i(q := quotient!:(x, pi4));
r := difference!:(x, times!:(q, pi4)) >>;
sign := i2bf!: 1;
if m >= 8 then m := remainder(m, 8);
if m >= 4 then << sign := minus!: sign; m := m - 4>>;
if m = 0 then go to sn
else if onep m then go to m1
else if m = 2 then go to m2
else go to m3;
m1: r := cut!:mt(difference!:(pi4, r), k2);
return times!:(sign, cos!:(r, k));
m2: r := cut!:mt(r, k2);
return times!:(sign, cos!:(r, k));
m3: r := cut!:mt(difference!:(pi4, r), k2);
sn: begin integer j, n, ncut; scalar dcut, fctrial, ri, tm;
ncut := k2 - min(0, order!: r + 1);
dcut := make!:bf(10, - ncut);
y := ri := tm := r;
r := minus!: cut!:ep(times!:(r, r), - ncut);
m := j := 1;
while greaterp!:(abs!: tm, dcut) do
<< j := j + 2;
fctrial := i2bf!:(m := m*j*(j - 1));
ri := cut!:ep(times!:(ri, r), - ncut);
n := max(1,k2 - order!: fctrial + order!: ri);
tm := divide!:(ri, fctrial, n);
y := plus!:(y, tm);
if zerop remainder(j,20) then y := cut!:ep(y,-ncut)>>;
end;
return round!:mt(times!:(sign, y), k);
end$
symbolic procedure cos!:(x, k);
% This function calculates cos(x), the value of
% the cosine function at the point "x", with
% the precision K, by summing terms of the
% Taylor series for cos(z), 0 < z < PI/4.
% X is a BIG-FLOAT representation of any real "x".
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg 'cos!:
else if bfzerop!: x then i2bf!: 1
else if minusp!: x then cos!:(minus!: x, k)
else begin integer k2, m; scalar pi4, sign, q, r, y;
k2 := k + 2;
m := preci!: x;
pi4 := times!:(!:pi(k2 + m), !:bf!-0!.25);
if lessp!:(x, pi4) then << m := 0; r := x >>
else << m := conv!:bf2i(q := quotient!:(x, pi4));
r := difference!:(x, times!:(q, pi4)) >>;
sign := i2bf!: 1;
if m >= 8 then m := remainder(m, 8);
if m >= 4 then << sign := minus!: sign; m := m - 4 >>;
if m >= 2 then sign := minus!: sign;
if m = 0 then go to cs
else if m = 1 then go to m1
else if m = 2 then go to m2
else go to m3;
m1: r := cut!:mt(difference!:(pi4, r), k2);
return times!:(sign, sin!:(r, k));
m2: r := cut!:mt(r, k2);
return times!:(sign, sin!:(r, k));
m3: r := cut!:mt(difference!:(pi4, r), k2);
cs: begin integer j, n; scalar dcut, fctrial, ri, tm;
dcut := make!:bf(10, - k2);
y := ri := tm := i2bf!: 1;
r := minus!: cut!:ep(times!:(r, r), - k2);
m := 1;
j := 0;
while greaterp!:(abs!: tm, dcut) do
<< j := j + 2;
fctrial := i2bf!:(m := m * j * (j - 1));
ri := cut!:ep(times!:(ri, r), - k2);
n := max(1, k2 - order!: fctrial + order!: ri);
tm := divide!:(ri, fctrial, n);
y := plus!:(y, tm);
if zerop remainder(j,20) then y := cut!:ep(y,-k2)>>;
end;
return round!:mt(times!:(sign, y), k);
end$
symbolic procedure tan!:(x, k);
% This function calculates tan(x), the value of
% the tangent function at the point "x",
% with the precision K, by calculating
% sin(x) or cos(x) = sin(PI/2-x).
% X is a BIG-FLOAT representation of any real "x",
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg 'tan!:
else if bfzerop!: x then i2bf!: 0
else if minusp!: x then minus!: tan!:(minus!: x, k)
else begin integer k2, m; scalar one, pi4, sign, q, r;
k2 := k + 2;
one := i2bf!: 1;
m := preci!: x;
pi4 := times!:(!:pi(k2 + m), !:bf!-0!.25);
if lessp!:(x, pi4) then << m := 0; r := x >>
else << m := conv!:bf2i(q := quotient!:(x, pi4));
r := difference!:(x, times!:(q, pi4)) >>;
if m >= 4 then m := remainder(m, 4);
if m >= 2 then sign := minus!: one else sign := one;
if m = 1 or m = 3 then r := difference!:(pi4, r);
r := cut!:mt(r, k2);
if m = 0 or m = 3 then go to m03 else go to m12;
m03: r := sin!:(r, k2);
q := difference!:(one, times!:(r, r));
q := sqrt!:(cut!:mt(q, k2), k2);
return times!:(sign, divide!:(r, q, k));
m12: r := sin!:(r, k2);
q := difference!:(one, times!:(r, r));
q := sqrt!:(cut!:mt(q, k2), k2);
return times!:(sign, divide!:(q, r, k));
end$
symbolic procedure asin!:(x, k);
% This function calculates asin(x), the value of
% the arcsine function at the point "x",
% with the precision K, by calculating
% atan(x/SQRT(1-x**2)) by ATAN!:.
% The answer is in the range [-PI/2 , PI/2].
% X is a BIG-FLOAT representation of "x", IxI <= 1;
% K is a positive integer.
if greaterp!:(abs!: x, i2bf!: 1) or
not fixp k or k <= 0 then bflerrmsg 'asin!:
else if minusp!: x then minus!: asin!:(minus!: x, k)
else begin integer k2; scalar one, y;
k2 := k + 2;
one := i2bf!: 1;
if lessp!:(difference!:(one, x), make!:bf(10, - k2))
then return round!:mt(times!:(!:pi add1 k,!:bf!-0!.5),k);
%JBM
y := cut!:mt(difference!:(one, times!:(x, x)), k2);
y := divide!:(x, sqrt!:(y, k2), k2);
return atan!:(y, k);
end$
symbolic procedure acos!:(x, k);
% This function calculates acos(x), the value of
% the arccosine function at the point "x",
% with the precision K, by calculating
% atan(SQRT(1-x**2)/x) if x > 0 or
% atan(SQRT(1-x**2)/x) + PI if x < 0.
% The answer is in the range [0 , PI].
% X is a BIG-FLOAT representation of "x", IxI <= 1.
% K is a positive integer.
if greaterp!:(abs!: x, i2bf!: 1) or
not fixp k or k <= 0 then bflerrmsg 'acos!:
else begin integer k2; scalar y;
k2 := k + 2;
if lessp!:(abs!: x, make!:bf(50, - k2))
then return round!:mt(times!:(!:pi add1 k,!:bf!-0!.5),k);
%JBM
y := difference!:(i2bf!: 1, times!:(x, x));
y := cut!:mt(y, k2);
y := divide!:(sqrt!:(y, k2), abs!: x, k2);
return (if minusp!: x then
round!:mt(difference!:(!:pi(k + 1), atan!:(y, k)), k)
else atan!:(y, k) );
end$
symbolic procedure atan!:(x, k);
% This function calculates atan(x), the value of the
% arctangent function at the point "x", with
% the precision K, by summing terms of the
% Taylor series for atan(z) if 0 < z < 0.42.
% Otherwise the following identities are used:
% atan(x) = PI/2 - atan(1/x) if 1 < x and
% atan(x) = 2*atan(x/(1+SQRT(1+x**2)))
% if 0.42 <= x <= 1.
% The answer is in the range [-PI/2 , PI/2].
% X is a BIG-FLOAT representation of any real "x".
% K is a positive integer.
if not fixp k or k <= 0 then bflerrmsg 'atan!:
else if bfzerop!: x then i2bf!: 0
else if minusp!: x then minus!: atan!:(minus!: x, k)
else begin integer k2; scalar one, pi4, y, z;
k2 := k + 2;
one := i2bf!: 1;
pi4 := times!:(!:pi k2, !:bf!-0!.25); %JBM
if equal!:(x, one) then return round!:mt(pi4, k);
if greaterp!:(x, one) then return
round!:mt(difference!:(plus!:(pi4, pi4),
atan!:(divide!:(one,x,k2),k + 1)),k);
if lessp!:(x, !:bf!-0!.42) then go to at;
y := plus!:(one, cut!:mt(times!:(x, x), k2));
y := plus!:(one, sqrt!:(y, k2));
y := atan!:(divide!:(x, y, k2), k + 1);
return round!:mt(times!:(y, i2bf!: 2), k);
at: begin integer m, n, ncut; scalar dcut, tm, zi;
ncut := k2 - min(0, order!: x + 1);
y := tm := zi := x;
z := minus!: cut!:ep(times!:(x, x), - ncut);
dcut := make!:bf(10, - ncut);
m := 1;
while greaterp!:(abs!: tm, dcut) do
<< zi := cut!:ep(times!:(zi, z), - ncut);
n := max(1, k2 + order!: zi);
tm := divide!:(zi, i2bf!:(m := m + 2), n);
y := plus!:(y, tm);
if zerop remainder(m,20) then y := cut!:ep(y,-ncut)>>;
end;
return round!:mt(y, k)
end;
%Miscellaneous constants (added by JBM).
!:bf!-pi := make!:bf(314159265358979323846, -20);
!:bf!-e := make!:bf(271828182845904523536, -20);
!:bf!-0!.25 := make!:bf(25,-2);
!:bf!-0!.5 := make!:bf(5,-1);
!:bf!-0!.1 := make!:bf(1,-1);
!:bf!-1!.72 := make!:bf(172,-2);
!:bf!-0!.42 := make!:bf(42,-2);
endmodule;
end;