Artifact c8ef1a864101b26aa26644bc3b78611589c5c251581f67d16be5d7bf7f94d22f:
- File
r34.1/src/arith.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 132889) [annotate] [blame] [check-ins using] [more...]
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; % u is always bigfloat. % 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,m; x := 0.0; 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(); if ep<0 then <<ep := ep+!!nfpd; m := t>>; if (ep := 10.0**ep)=0.0 then error1(); y := u := mt!: u; lp: x := x + (y - 10000*(u := u/10000))*ep; if u=0 then return if not m then x else x/!!flint; 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.0/cos x; symbolic procedure csc x; 1.0/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!-plus(x,y); if zerop x then y else if zerop y then x else if x>0.0 and y>0.0 then if x<!!plumax and y<!!plumax then plus2(x,y) else nil else if x<0.0 and y<0.0 then if -x<!!plumax and -y<!!plumax then plus2(x,y) else nil else if abs x<!!plumin and abs y<!!plumin then nil else (if u=0.0 then u else if abs u<!!fleps1*abs x then 0.0 else u) where u = plus2(x,y); 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('(cotd 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 abs 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 := 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!*) (atan2 cratan2!*) (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 cratan2!*(y,x); begin scalar q,p; q := crsqrt!* cr!:plus(cr!:times(y,y),cr!:times(x,x)); if cr!:zerop q then error(0,list("invalid arguments to ",'atan2)); y := cr!:quotient(y,q); x := cr!:quotient(x,q); p := rdpi!*(); if cr!:zerop x then <<q := rd!:quotient(p,i2rd!* 2); return if bfminusp retag crrl y then rd!:minus q else q>>; q := cratan!* cr!:quotient(y,x); if bfminusp retag crrl x then <<p := !*rd2cr p; q := if bfminusp retag crrl y then cr!:differ(q,p) else cr!:plus(q,p)>>; % bfzp x is probably impossible? return q end; 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;