module smlbflot; % Basic support for bigfloat arithmetic.
% Authors: S.L. Kameny and T. Sasaki.
% Modified for binary bigfloat arithmetic by Iain Beckingham and Rainer
% Schoepf.
% Modified for double precision printing by Herbert Melenk.
% Modified to allow *very* large numbers to be compressed (for PSL) by
% Winfried Neun.
% Last change made Oct 6, 1999.
exports abs!:, bfexplode0, bflerrmsg, bfprin!:, bftrim!:, bfzerop!:,
conv!:mt, cut!:ep, cut!:mt, decimal2internal, decprec!:,
difference!:, divide!:, equal!:, fl2bf, greaterp!:, incprec!:,
leq!:, lessp!:, max!:, max!:, max2!:, min!:, min2!:, minus!:,
minusp!:, order!:, plus!:, read!:num, round!:mt, round!:last,
times!:;
imports aconc, ashift, bfp!:, ceiling, conv!:bf2i, ep!:, eqcar, floor,
geq, i2bf!:, leq, lshift, make!:ibf, msd!:, mt!:, neq, normbf,
oddintp, preci!:, precision, prin2!*, rerror, retag, reversip;
fluid '(!*bfspace !*fullprec !*nat !:prec!: !:bprec!: !:print!-prec!:
!:upper!-sci!: !:lower!-sci!:);
global '(!!nfpd !!nbfpd bften!* bfz!* fort_exponent);
switch bfspace,fullprec;
flag('(fort_exponent),'share);
!*bfspace := nil;
% !*fullprec := t;
!:lower!-sci!: := 10;
!:upper!-sci!: := 5;
symbolic procedure bflerrmsg u;
% Revised error message for BFLOAT module, using error, not rederr.
error(0,{"Invalid argument to",u});
symbolic procedure bfzerop!: u;
% This is possibly too restricted a definition.
mt!: u = 0;
symbolic procedure fl2bf x;
(if zerop x then bfz!* else
begin scalar s,r; integer d;
if x<0 then <<x := -x; s := t>>;
% convert x to an integer equivalent;
r := normbf read!:num x;
d := ep!: r+msd!: mt!: r;
x := x*2.0**-d; x := x + 0.5/2**!:bprec!:;
x := fix(x*2**!:bprec!:);
return make!:ibf (if s then -x else x, d - !:bprec!:) end)
where !:bprec!:=!!nbfpd;
symbolic procedure bfprin!: u;
% if preci!: u>!!nbfpd then bfprin0 u
% else (bfprin0 u where !*bfspace=nil);
bfprin0 u;
symbolic procedure divide!-by!-power!-of!-ten (x, n);
if n < 0 then bflerrmsg 'divide!-by!-power!-of!-ten
else <<
while n > 0 do <<
if oddintp n then x := normbf divide!: (x, f, !:bprec!:);
n := lshift (n, -1);
f := normbf cut!:mt (times!: (f, f), !:bprec!:) >>;
x >> where f := bften!*;
symbolic procedure multiply!-by!-power!-of!-ten (x, n);
if n < 0 then bflerrmsg 'multiply!-by!-power!-of!-ten
else <<
while n > 0 do <<
if oddintp n then x := normbf times!: (x, f);
n := lshift (n, -1);
f := normbf cut!:mt (times!: (f, f), !:bprec!:) >>;
normbf cut!:mt (x, !:bprec!:) >> where f := bften!*;
global '(log2of10);
symbolic procedure round!:dec (x, p);
%
% rounds bigfloat x to p decimal places
%
begin scalar setpr; integer m, ex;
if null p then p := if !:print!-prec!: then !:print!-prec!:
else !:prec!: - 2
else if p > precision 0 then setpr := precision p;
x := round!:dec1 (x,p);
m := car x; ex := cdr x;
x := i2bf!: m;
if ex < 0 then x := divide!-by!-power!-of!-ten (x, -ex)
else if ex > 0 then x := multiply!-by!-power!-of!-ten (x, ex);
if setpr then precision setpr;
return round!:mt (x, ceiling (p * log2of10))
end;
symbolic procedure round!:dec1 (x, p);
%
% rounds bigfloat x to p decimal places
% returns pair (m . ex) of mantissa and exponent to base 10,
% m having exactly p digits
% performs all calculations at at least current precision,
% but increases the precision of the calculations to log10(x)
% if this is larger
%
if bfzerop!: x then cdr x
else (begin scalar exact, lo, sign; integer ex, k, m, n, l;
%
% We need to calculate the number k so that 10^(k+1) > |x| >= 10^k
% k = floor (log10 |x|) = floor (log2 |x| / log2of10);
% We can easily compute n so that 2^(n+1) > |x| >= 2^n,
% i.e., n = floor (log2 |x|), since this is just order!:(x).
% Since n+1 > log2 |x| >= n, it follows that
% floor ((n+1) / log2of10) >= k >= floor (n / log2of10)
% I.e., if both bounds agree, we know k, otherwise we have to check.
%
if mt!: x < 0 then <<sign := t; x := minus!: x>>;
n := order!: x;
%
% The division by log2of10 has to be done with precision larger than
% the precision of n. In particular, log2of10 has to be calculated
% to a larger precision. Instead of dividing by log2of10, we
% multiply by log10of2.
%
l := msd!: abs n;
<<lo := divide!: (!:log2 !:bprec!:, !:log10 !:bprec!:, !:bprec!:);
k := conv!:bf2i times!: (i2bf!: n, lo);
if k = conv!:bf2i times!: (i2bf!: (n+1), lo) then exact := t>>
where !:bprec!: := max (!!nbfpd, l + 7);
%
% For the following calculation the precision must be increased by
% the precision of n. The is necessary to ensure that the mantissa
% is calculated correctly for large values of the exponent. This is
% due to the fact that if we multiply the number x by 10^n its
% precision will be decreased by n.
%
!:bprec!: := !:bprec!: + l;
%
% since conv!:bf2i rounds always towards 0, we must correct for n<0
%
if n < 0 then k := k - 1;
ex := k - p + 1;
if ex < 0 then x := multiply!-by!-power!-of!-ten (x, -ex)
else if ex > 0 then x := divide!-by!-power!-of!-ten (x, ex);
if exact then nil
else <<l := length explode conv!:bf2i x;
if l < p then <<x := times!: (x, bften!*); ex := ex - 1>>
else if l > p then <<x := divide!: (x, bften!*, !:bprec!:);
ex := ex + 1>>>>;
%
% do rounding
%
x := plus!:(x, bfhalf!*);
% Add an "epsilon" just to be sure (e.g., for on complex,rounded;
% print_precision 15; 3.23456789012345+7).
x := plus!:(x,fl2bf(0.1^18));
m := conv!:bf2i x;
if length explode m > p then <<m := (m + 5) / 10; ex := ex + 1>>;
if sign then m := -m;
return (m . ex);
end) where !:bprec!: := !:bprec!:;
% symbolic procedure internal2decimal (x, p);
%
% converts bigfloat x to decimal format, with precision p
% Result is a pair (m . e), so that x = m*10^e, with
% m having exactly p decimal digits.
% Calculation is done with the current precision,
% but at least with p + 2.
%
% begin scalar setpr;
% if null p then p := if !:print!-prec!: then !:print!-prec!:
% else !:prec!: - 2
% else if p > precision 0 then setpr := precision p;
% x := round!:dec1 (x,p);
% if setpr then precision setpr;
% return x
% end;
symbolic procedure bfprin0 u;
begin scalar r; integer m, ex;
r := round!:dec1 (u, if !:print!-prec!: then !:print!-prec!:
else !:prec!: - 2);
m := car r; ex := cdr r;
bfprin0x (m, ex)
end;
symbolic procedure bfprin0x(m,ex);
begin scalar lst; integer dotpos;
lst := bfexplode0x(m,ex);
ex := cadr lst;
dotpos := caddr lst;
lst := car lst;
return bfprin!:lst (lst,ex,dotpos)
end;
symbolic procedure bfexplode0 u;
% returns a list (lst ex dotpos) where
% lst = list of characters in mantissa
% (ie optional sign and digits)
% ex = decimal exponent
% dotpos = position of decimal point in lst
% (note that the sign is counted)
begin scalar r; integer m, ex;
r := round!:dec1 (u,if !:print!-prec!: then !:print!-prec!:
else !:prec!: - 2);
m := car r; ex := cdr r;
return bfexplode0x (m, ex)
end;
symbolic procedure bfexplode0x (m, ex);
begin scalar lst, s; integer dotpos, l;
if m<0 then <<s := t; m := -m>>;
lst := explode m;
l := length lst;
if ex neq 0 and (l+ex < -!:lower!-sci!: or l+ex > !:upper!-sci!:)
then <<dotpos := 1;
ex := ex + l - 1>>
else <<dotpos := l + ex;
ex := 0;
if dotpos > l - 1
%
% add dotpos - l + 1 zeroes at the end
%
then lst := nconc!*(lst,nlist('!0,dotpos - l + 1))
else while dotpos < 1 do <<lst := '!0 . lst;
dotpos := dotpos + 1>>>>;
if s then <<lst := '!- . lst; dotpos := dotpos + 1>>;
if null !*fullprec
then <<lst := reversip lst;
while lst and car lst eq '!0 and length lst > dotpos + 1 do
lst := cdr lst;
lst := reversip lst>>;
return {lst, ex, dotpos}
end;
symbolic procedure bfprin!:lst (lst, ex, dotpos);
begin scalar result,ee,w; integer j;
ee:='E;
if !*fort and liter(w:=reval fort_exponent) then ee:=w else w:=nil;
if car lst eq '!- and dotpos = 1 then <<dotpos := 2; ex := ex - 1>>;
if ex neq 0 then if car lst eq '!- then <<ex := ex + dotpos - 2;
dotpos := 2>>
else <<ex := ex + dotpos - 1; dotpos := 1>>
else if dotpos = length lst then dotpos := -1;
for each char in lst do <<
result := char . result; j := j + 1; dotpos := dotpos - 1;
if j=5 then <<if !*nat and !*bfspace then result := '! . result;
j := 0>>;
if dotpos = 0 then <<result := '!. . result; j := j + 1>>;
if j=5 then <<if !*nat and !*bfspace then result := '! . result;
j := 0>>>>;
if ex neq 0 or w then <<
if not (!*nat and !*bfspace) then result := ee . result
else if j=0 then <<result := '! . ee . result; j := 2>>
else if j=1 then <<result := '! . ee . '! . result; j := 4>>
else if j=2
then <<result := '! . '! . ee . '! . result; j := 0>>
else if j=3 then <<result := '! . ee . '! . result; j := 0>>
else if j=4 then <<result := '! . ee . '! . result; j := 2>>;
lst := if ex > 0 then '!+ . explode ex else explode ex;
for each char in lst do <<
result := char . result; j := j + 1;
if j=5 then <<if !*nat and !*bfspace then result := '! . result;
j := 0>>>>>>;
% if !*nat then for each char in reversip result do prin2!* char else
prin2!* smallcompress reversip result
end;
symbolic procedure smallcompress (li);
begin scalar ll;
if (ll := length li)>1000
then <<li := smallsplit(li,ll/2);
return concat(smallcompress car li,smallcompress cdr li)>>
else return compress ('!" . reversip('!" . reversip li))
end;
symbolic procedure smallsplit (li,ln);
begin scalar aa;
for i:=1:ln do <<aa := car li . aa; li := rest li>>;
return (reversip aa) . li
end;
symbolic procedure scientific_notation n;
begin scalar oldu,oldl;
oldu := !:upper!-sci!:; oldl := !:lower!-sci!: + 1;
if fixp n
then <<
if n<0
then rerror(arith,1,
{"Invalid argument to scientific_notation:",n});
!:lower!-sci!: := n - 1; !:upper!-sci!: := n;
>>
else if eqcar(n,'list) and length n=3
then if not (fixp cadr n and fixp caddr n)
then rerror(arith,2,
{"Invalid argument to scientific_notation:",n})
else <<!:upper!-sci!: := cadr n;
!:lower!-sci!: := caddr n - 1 >>;
return {'list,oldu,oldl} % Return previous range.
end;
flag('(scientific_notation), 'opfn);
symbolic procedure order!: nmbr;
% This function counts the order of a number "n". NMBR is a bigfloat
% representation of "n".
% **** ORDER(n)=k if 2**k <= ABS(n) < 2**(k+1)
% **** when n is not 0, and ORDER(0)=0.
%
if mt!: nmbr = 0 then 0 else preci!: nmbr + ep!: nmbr - 1;
symbolic smacro procedure decprec!:(nmbr, k);
make!:ibf(ashift(mt!: nmbr,-k), ep!: nmbr + k);
symbolic smacro procedure incprec!:(nmbr, k);
make!:ibf(ashift(mt!: nmbr,k), ep!: nmbr - k);
symbolic procedure conv!:mt(nmbr, k);
% This function converts a number "n" to an equivalent number of
% binary precision K by rounding "n" or adding "0"s to "n".
% NMBR is a binary 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 binary 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!: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 BINARY 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$
symbolic procedure round!:last nmbr;
% This function rounds a number "n" at its last place.
% NMBR is a binary bigfloat representation of "n".
<< if m < 0 then << m := -m; s := t >>;
m := if oddintp m then lshift (m, -1) + 1
else lshift (m, -1);
if s then m := -m;
make!:ibf (m, e) >>
where m := mt!: nmbr, e := ep!: nmbr + 1, s := nil;
symbolic procedure cut!:mt(nmbr,k);
% This function returns a given number "n" unchanged
% if its binary 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 BINARY BIG-FLOAT 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 decprec!:(nmbr, k)
else bflerrmsg 'cut!:mt$
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 BINARY 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 bftrim!: v;
normbf round!:mt(v,!:bprec!: - 3);
symbolic procedure decimal2internal (base10,exp10);
if exp10 >= 0 then i2bf!: (base10 * 10**exp10)
else divide!-by!-power!-of!-ten (i2bf!: base10, -exp10);
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 170.
% **** 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!:ibf(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 := - smallcompress cdr l
else m:= smallcompress l;
make: u := reverse u;
v := reverse v;
if appear!/ then
return conv!:r2bf(make!:ratnum(sign*compress v,compress u),
if !:bprec!: then !:bprec!: else 170);
if appear!. then j := - j else j := 0;
if sign = 1 then u := compress u else u := - compress u;
return round!:mt (decimal2internal (u, j + m), !:bprec!:)
where !:bprec!: := if !:bprec!: then !:bprec!:
else msd!: abs u
end;
symbolic procedure abs!: nmbr;
% This function makes the absolute value of "n". N is a binary
% bigfloat representation of "n".
if mt!: nmbr > 0 then nmbr else make!:ibf(- mt!: nmbr, ep!: nmbr);
symbolic procedure minus!: nmbr;
% This function makes the minus number of "n". N is a binary
% bigfloat representation of "n".
make!:ibf(- 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!:ibf(m1+m2, e1)
else if d>0 then make!:ibf(ashift(m1,d)+m2,e2)
else make!:ibf(m1+ashift(m2,-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!:ibf(m1 - m2, e1)
else if d>0 then make!:ibf(ashift(m1,d) - m2,e2)
else make!:ibf(m1 - ashift(m2,-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!:ibf(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!:ibf(mt!: n1 / mt!: n2, ep!: n1 - ep!: n2);
return round!:mt(n1, k)
end;
symbolic procedure max2!:(a,b);
% This function returns the larger of "n1" and "n2".
% N1 and N2 are bigfloat representations of "n1" and "n2".
if greaterp!:(a,b) then a else b;
macro procedure max!: x; expand(cdr x,'max2!:);
symbolic procedure min2!:(a,b);
% This function returns the smaller of "n1" and "n2".
% N1 and N2 are binary bigfloat representations of "n1" and "n2".
if greaterp!:(a,b) then b else a;
macro procedure min!: x; expand(cdr x,'min2!:);
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 ashift(ma,d)>mb
else ma>ashift(mb,-d))
where p2=2*!:bprec!:))
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,3,"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!:ibf( numr!: rnmbr, 0),
make!:ibf( denm!: rnmbr, 0),k)
else bflerrmsg 'conv!:r2bf;
endmodule;
end;