module rounded; % *** Support for Arbitrary Rounded Arithmetic.
% Authors: Anthony C. Hearn and Stanley L. Kameny.
% Last updated: 23 June 1993.
% Copyright (c) 2000, Anthony C. Hearn. 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);
exports chkint!*, chkrn!*, convprec, convprec!*, deg2rad!*,
i2rd!*, logfp, mkround, rd!:difference, rd!:minus, rd!:minusp,
rd!:onep, rd!:plus, rd!:prep, rd!:prin, rd!:quotient,
rd!:simp, rd!:times, rd!:zerop, rdprep1, rdqoterr, rdzchk,
rndbfon, round!*, roundbfoff, roundbfon, roundconstants,
safe!-fp!-times;
imports !*d2q, !:difference, !:minus, !:minusp, !:zerop, abs!:, aeval,
apply1, bf2flr, bfdiffer, bfexplode0, bfinverse, bflessp,
bfloat, bfminus, bfminusp, bfprin!:, bftrim!:, bfzerop!:,
bfzp, ceiling, copyd, deg2rad!*, difbf, divbf, dmoderr,
ep!:, eqcar, equal!:, errorp, errorset!*, fl2int,
fl2rd, float!-bfp, floor, ft2rn1, geq, greaterp!:, grpbf,
i2bf!:, initdmode, invbf, leq, lessp!:, log, lprim, lshift,
make!:ibf, make!:rd, minus!:, minusp!:, mkquote, msgpri, mt!:,
neq, normbf, off1, on1, over, plubf, preci!:, r2bf, rd2fl,
rd!:forcebf, realrat, rerror, retag, rmsubs, round!:mt, setk,
sqrt, timbf, times!:, union;
fluid '(!:prec!: !:bprec!: !:print!-prec!: minprec!# rootacc!#!#);
fluid '(dmode!* !*bfspace !*numval !*roundbf !*!*roundbf !*norndbf);
fluid '(!*noconvert);
global '(bfone!* epsqrt!* log2of10);
global '(domainlist!* !!nfpd !!nbfpd !!flprec !!rdprec mxflbf!!
mnflbf!!);
global '(!!plumax !!plumin !!timmax !!timmin !!maxflbf !!minflbf
!!fleps1 !!fleps2 !!flint !!maxbflexp log2 !!maxarg);
global '(rd!-tolerance!* cr!-tolerance!* yy!! bfz!* !!smlsin);
switch rounded;
%Set value for !!flprec. It never changes.
!!flprec := !!nfpd - 3;
!!smlsin := 10.0^-(2+!!flprec);
symbolic procedure logfp x;
% floating log of x**(1/n) using bfloat logic as boost.
(log(m/float lshift(1,p))+(p+ep!: x)*log2)
where p=(preci!: x - 1) where m=mt!: x;
symbolic procedure roundconstants;
<<!!plumax := float(2**(!!maxbflexp -1));
!!minflbf := invbf(!!maxflbf := make!:ibf (1,!!maxbflexp));
% plumin must be large enough to avoid underflow from difference.
!!plumin := 10.0**!!flprec/!!plumax;
!!timmin := 1/(!!timmax := sqrt(!!plumax));
!!maxarg := logfp !!maxflbf>>;
switch bfspace,numval,roundbf; % norndbf.
!*bfspace := nil;
!*numval := t;
put('roundbf,'simpfg,'((t (roundbfon)) (nil (roundbfoff))));
symbolic procedure roundbfon; !*!*roundbf := t;
symbolic procedure roundbfoff; !*!*roundbf := !!rdprec > !!flprec;
% 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!:,'rootfn,'rd!:root);
put('!:rd!:,'!:rn!:,'!*rd2rn);
put('!:rn!:,'!:rd!:,'!*rn2rd);
symbolic procedure round!* x;
% Returns actual number representation, as either float or bfloat.
% retag cdr x;
if float!-bfp x then rd2fl x else x;
symbolic procedure mkround u;
% inverse operation to round!*, i.e. tags a naked float
if atom u then make!:rd u else u;
%symbolic procedure roundbfp; !*roundbf or !!rdprec > !!flprec;
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 precision0 n;
% called from algebraic call of precision.
if n member '((nil) () (reset))
then <<rootacc!#!# := nil; precision !!flprec>>
else if cdr n
or not numberp(n := prepsq simp!* aeval {'fix,prepsq simp!* car n})
or n<0 then
rerror(arith,5,"positive numeric value or `RESET' required")
else <<if n>0 then rootacc!#!# := max(n,6); precision n>>;
put('precision,'psopfn,'precision0);
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,6,"positive number required");
precision1(n,t)>>;
log2of10 := log 10 / log 2;
symbolic procedure decprec2internal p;
ceiling(p * log2of10) + 3;
% symbolic procedure internal2decprec p;
% floor ((p - 3) / log2of10);
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 := if !*roundbf then n else max(n,minprec!#))+2;
if !:print!-prec!: and n < !:print!-prec!:+2
then !:print!-prec!: := nil; %unset
!:bprec!: := decprec2internal !:prec!:;
epsqrt!* := make!:ibf(1, -!:bprec!:/2);
rd!-tolerance!* := make!:ibf(1, 6-!:bprec!:);
cr!-tolerance!* := make!:ibf(1, 2*(6-!:bprec!:));
% if !!rdprec <= !!flprec then
% <<!!fleps1 := 1.0/float(2**(!:bprec!: - 2));
% !!fleps2 := !!fleps1**2>>;
!*!*roundbf := !!rdprec > !!flprec or !*roundbf;
return oldprec end;
flag('(print!-precision),'opfn); % Symbolic operator print!-precision.
flag('(print_precision),'opfn); % Symbolic operator print_precision.
symbolic procedure !*rd2rn x;
% Converts a rounded number N into a rational to the system precision.
% Elegant form: uses both rd2rn1 and realrat... and choses the best,
% but uses a heuristic to avoid the extra work when not needed.
begin scalar n,p,r,r1,r2,d1,d2,ov;
if rd!:zerop x then return '!:rn!: . (0 . 1);
p := precision 0;
r := rd2rn1 x;
r1 := '!:rn!: . r;
if abs car r<10 or cdr r<10
or 2*max(length explode cdr r,length explode abs car r)<p+1
then go to ret;
r2 := '!:rn!: . realrat bftrim!: rd!:forcebf x;
precision(2+p);
d1 := !:difference(x,r1); if !:minusp d1 then d1 := !:minus d1;
d2 := !:difference(x,r2); if !:minusp d2 then d2 := !:minus d2;
if !:zerop d2 or !:minusp !:difference(d2,d1) then ov := t;
precision p;
ret: return if ov then r2 else r1 end;
symbolic procedure rd2rn1 n;
if float!-bfp n then ft2rn1 rd2fl n else bf2rn1 n;
symbolic procedure bf2rn1 n;
% Here, the nonzero input n is always a binary bigfloat
begin scalar negp,a,p0,p1,q0,q1,w,flagg,nn,r0,r1;
if mt!: n<0 then <<negp := t; n := minus!: n>>;
nn := n;
top: a := ((if d=0 then m else lshift(m,d))
where m=mt!: n,d=ep!: n);
n := difbf(n,normbf i2bf!: a);
if not flagg
then <<flagg := t; p0 := 1; p1 := a; q0 := 0; q1 := 1>>
else <<w := p0 + a*p1; p0 := p1; p1 := w; r0 := r1;
w := q0 + a*q1; q0 := q1; q1 := w>>;
r1 := abs!: difbf(nn,divbf(i2bf!: p1,i2bf!: q1));
% temporary write statement here
% if !*trrd2rn1 then << write p1 . q1," -> ",r1; terpri()>>;
if bfzerop!: n or bfzerop!: r1
then return if negp then (-p1) . q1 else p1 . q1
else if r0 and not greaterp!:(r0,r1)
then return if negp then (-p0) . q0 else p0 . q0;
n := invbf n;
go to top
end;
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.
%!!fleps1 := 1.0/float(2**(!:bprec!: - 6));
!!fleps1 := 2.0**(6 - !:bprec!:);
!!fleps2 := !!fleps1**2;
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 null atom u and car u=0 then nil ./ 1
else if null dmode!* or dmode!* eq '!:gi!:
then (if eqcar(x,'!:rn!:) then cdr x else x ./ 1)
where x = !*rd2rn make!: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,make!: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 msd!: x <= !!maxbflexp then float u
else <<rndbfon(); bfloat u>>) where x=abs u);
mnflbf!! := invbf(mxflbf!! := make!:ibf (1, 800));
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 retag u;
symbolic procedure convprec u; convchk round!* u;
symbolic procedure rd!:minusp u;
% bfminusp round!* u;
if float!-bfp u then minusp rd2fl u else minusp!: u;
symbolic procedure convprc2(u,v);
<<u := convprec u; yy!! := convprec 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 make!: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 errorp 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 make!: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 errorp 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 make!: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 errorp z then <<rndbfon(); timbf(bfloat x,bfloat y)>>
else car z>> end) where z=nil;
symbolic procedure rd!:quotient(u,v);
if !:zerop v then rerror(arith,7,"division by zero") else
(if not !*!*roundbf and atom cdr u and atom cdr v
and (z := safe!-fp!-quot(cdr u,cdr v)) then make!: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 errorp 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;
if float!-bfp u then zerop rd2fl u else mt!: u = 0;
symbolic procedure rd!:minus u;
% mkround bfminus round!* u;
if float!-bfp u then fl2rd (- rd2fl u) else minus!: 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 float!-bfp u then abs(1.0 - rd2fl u)<!!fleps1
else equal!:(bfone!*,bftrim!: u);
symbolic procedure rd!:root(u,n);
if float!-bfp u then fl2rd expt(rd2fl u,recip float n)
else texpt!:any(u,quotient!:(bfone!*,i2bf!: n));
% 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 !*noconvert then rdprep1 u
else if rd!:onep u then 1
else if rd!:onep rd!:minus u then -1 else rdprep1 u;
%symbolic procedure rdprep1 u;
% if float!-bfp u then
% if not dmode!* memq '(!:rd!: !:cr!:) or !*!*roundbf
% then round!:mt(bfloat rd2fl u,min(!:bprec!:,!!nbfpd))
% else if !:bprec!:>!!nbfpd then u
% else fl2rd bf2flr round!:mt(bfloat rd2fl u,!:bprec!:)
% else round!:mt(u,!:bprec!:);
symbolic procedure rdprep1 u;
% Using cdr u to get actual float leads to various glitches.
if float!-bfp u then u else round!:mt(u,!:bprec!:);
symbolic procedure rd!:prin u;
% Printed output is rounded to 2 fewer places than internal value.
bfprin!: bftrim!: rd!:forcebf u;
symbolic procedure rd!:explode u;
bfexplode0 bftrim!: rd!:forcebf u;
initdmode 'rounded;
endmodule;
end;