1module smlbflot; % Basic support for bigfloat arithmetic. 2 3% Authors: S.L. Kameny and T. Sasaki. 4% Modified for binary bigfloat arithmetic by Iain Beckingham and Rainer 5% Schoepf. 6% Modified for double precision printing by Herbert Melenk. 7% Modified to allow *very* large numbers to be compressed (for PSL) by 8% Winfried Neun. 9 10% Redistribution and use in source and binary forms, with or without 11% modification, are permitted provided that the following conditions are met: 12% 13% * Redistributions of source code must retain the relevant copyright 14% notice, this list of conditions and the following disclaimer. 15% * Redistributions in binary form must reproduce the above copyright 16% notice, this list of conditions and the following disclaimer in the 17% documentation and/or other materials provided with the distribution. 18% 19% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 20% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 21% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 22% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 23% CONTRIBUTORS 24% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 25% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 26% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 27% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 28% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 29% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 30% POSSIBILITY OF SUCH DAMAGE. 31% 32 33 34% Last change made Oct 6, 1999. 35 36exports abs!:, bfexplode0, bflerrmsg, bfprin!:, bftrim!:, bfzerop!:, 37 conv!:mt, cut!:ep, cut!:mt, decimal2internal, decprec!:, 38 difference!:, divide!:, equal!:, fl2bf, greaterp!:, incprec!:, 39 leq!:, lessp!:, max!:, max!:, max2!:, min!:, min2!:, minus!:, 40 minusp!:, order!:, plus!:, read!:num, round!:mt, round!:last, 41 times!:; 42 43imports aconc, ashift, bfp!:, ceiling, conv!:bf2i, ep!:, eqcar, floor, 44 geq, i2bf!:, leq, lshift, make!:ibf, msd!:, mt!:, neq, normbf, 45 oddintp, preci!:, precision, prin2!*, rerror, retag, reversip; 46 47fluid '(!*bfspace !*fullprec !*nat !:prec!: !:bprec!: !:print!-prec!: 48 !:upper!-sci!: !:lower!-sci!:); 49 50global '(!!nfpd !!nbfpd bften!* bfz!* fort_exponent); 51 52switch bfspace,fullprec; 53 54flag('(fort_exponent),'share); 55 56!*bfspace := nil; 57 58% !*fullprec := t; 59 60!:lower!-sci!: := 10; 61 62!:upper!-sci!: := 5; 63 64symbolic procedure bflerrmsg u; 65 % Revised error message for BFLOAT module, using error, not rederr. 66 error(0,{"Invalid argument to",u}); 67 68symbolic procedure bfzerop!: u; 69 % This is possibly too restricted a definition. 70 mt!: u = 0; 71 72symbolic procedure fl2bf x; 73 (if zerop x then bfz!* else 74 if not fp!-finite x then rederr "Floating point infinity or NaN" else 75 begin scalar s,r; integer d; 76 if x<0 then <<x := -x; s := t>>; 77 % convert x to an integer equivalent; 78 r := normbf read!:num x; 79 d := ep!: r+msd!: mt!: r; 80 x := x*2.0**-d; x := x + 0.5/2.0**!:bprec!:; 81 x := fix(x*2.0**!:bprec!:); 82 return make!:ibf (if s then -x else x, d - !:bprec!:) end) 83 where !:bprec!:=!!nbfpd; 84 85symbolic procedure bfprin!: u; 86% if preci!: u>!!nbfpd then bfprin0 u 87% else (bfprin0 u where !*bfspace=nil); 88 bfprin0 u; 89 90symbolic procedure divide!-by!-power!-of!-ten (x, n); 91 if n < 0 then bflerrmsg 'divide!-by!-power!-of!-ten 92 else << 93 while n > 0 do << 94 if oddintp n then x := normbf divide!: (x, f, !:bprec!:); 95 n := lshift (n, -1); 96 f := normbf cut!:mt (times!: (f, f), !:bprec!:) >>; 97 x >> where f := bften!*; 98 99symbolic procedure multiply!-by!-power!-of!-ten (x, n); 100 if n < 0 then bflerrmsg 'multiply!-by!-power!-of!-ten 101 else << 102 while n > 0 do << 103 if oddintp n then x := normbf times!: (x, f); 104 n := lshift (n, -1); 105 f := normbf cut!:mt (times!: (f, f), !:bprec!:) >>; 106 normbf cut!:mt (x, !:bprec!:) >> where f := bften!*; 107 108global '(!!log2of10); 109 110symbolic procedure round!:dec (x, p); 111 % 112 % rounds bigfloat x to p decimal places 113 % 114 begin scalar setpr; integer m, ex; 115 if null p then p := if !:print!-prec!: then !:print!-prec!: 116 else !:prec!: - 2 117 else if p > precision 0 then setpr := precision p; 118 x := round!:dec1 (x,p); 119 m := car x; ex := cdr x; 120 x := i2bf!: m; 121 if ex < 0 then x := divide!-by!-power!-of!-ten (x, -ex) 122 else if ex > 0 then x := multiply!-by!-power!-of!-ten (x, ex); 123 if setpr then precision setpr; 124 return round!:mt (x, ceiling (p * !!log2of10)) 125 end; 126 127symbolic procedure round!:dec1 (x, p); 128 % 129 % rounds bigfloat x to p decimal places 130 % returns pair (m . ex) of mantissa and exponent to base 10, 131 % m having exactly p digits 132 % performs all calculations at at least current precision, 133 % but increases the precision of the calculations to log10(x) 134 % if this is larger 135 % 136 if bfzerop!: x then cdr x 137 else (begin scalar exact, lo, sign; integer ex, k, m, n, l; 138 % 139 % We need to calculate the number k so that 10^(k+1) > |x| >= 10^k 140 % k = floor (log10 |x|) = floor (log2 |x| / log2of10); 141 % We can easily compute n so that 2^(n+1) > |x| >= 2^n, 142 % i.e., n = floor (log2 |x|), since this is just order!:(x). 143 % Since n+1 > log2 |x| >= n, it follows that 144 % floor ((n+1) / log2of10) >= k >= floor (n / log2of10) 145 % I.e., if both bounds agree, we know k, otherwise we have to check. 146 % 147 if mt!: x < 0 then <<sign := t; x := minus!: x>>; 148 n := order!: x; 149 % 150 % The division by log2of10 has to be done with precision larger than 151 % the precision of n. In particular, log2of10 has to be calculated 152 % to a larger precision. Instead of dividing by log2of10, we 153 % multiply by log10of2. 154 % 155 l := msd!: abs n; 156 <<lo := divide!: (!:log2 !:bprec!:, !:log10 !:bprec!:, !:bprec!:); 157 k := conv!:bf2i times!: (i2bf!: n, lo); 158 if k = conv!:bf2i times!: (i2bf!: (n+1), lo) then exact := t>> 159 where !:bprec!: := max (!!nbfpd, l + 7); 160 % 161 % For the following calculation the precision must be increased by 162 % the precision of n. The is necessary to ensure that the mantissa 163 % is calculated correctly for large values of the exponent. This is 164 % due to the fact that if we multiply the number x by 10^n its 165 % precision will be decreased by n. 166 % 167 !:bprec!: := !:bprec!: + l; 168 % 169 % since conv!:bf2i rounds always towards 0, we must correct for n<0 170 % 171 if n < 0 then k := k - 1; 172 ex := k - p + 1; 173 if ex < 0 then x := multiply!-by!-power!-of!-ten (x, -ex) 174 else if ex > 0 then x := divide!-by!-power!-of!-ten (x, ex); 175 if exact then nil 176 else <<l := length explode conv!:bf2i x; 177 if l < p then <<x := times!: (x, bften!*); ex := ex - 1>> 178 else if l > p then <<x := divide!: (x, bften!*, !:bprec!:); 179 ex := ex + 1>>>>; 180 % 181 % do rounding 182 % 183 x := plus!:(x, bfhalf!*); 184 % Add an "epsilon" just to be sure (e.g., for on complex,rounded; 185 % print_precision 15; 3.23456789012345+7). 186 x := plus!:(x,fl2bf(0.1^18)); 187 m := conv!:bf2i x; 188 if length explode m > p then <<m := (m + 5) / 10; ex := ex + 1>>; 189 if sign then m := -m; 190 return (m . ex); 191 end) where !:bprec!: := !:bprec!:; 192 193% symbolic procedure internal2decimal (x, p); 194 % 195 % converts bigfloat x to decimal format, with precision p 196 % Result is a pair (m . e), so that x = m*10^e, with 197 % m having exactly p decimal digits. 198 % Calculation is done with the current precision, 199 % but at least with p + 2. 200 % 201% begin scalar setpr; 202% if null p then p := if !:print!-prec!: then !:print!-prec!: 203% else !:prec!: - 2 204% else if p > precision 0 then setpr := precision p; 205% x := round!:dec1 (x,p); 206% if setpr then precision setpr; 207% return x 208% end; 209 210 211symbolic procedure bfprin0 u; 212 begin scalar r; integer m, ex; 213 r := round!:dec1 (u, if !:print!-prec!: then !:print!-prec!: 214 else !:prec!: - 2); 215 m := car r; ex := cdr r; 216 bfprin0x (m, ex) 217 end; 218 219symbolic procedure bfprin0x(m,ex); 220 begin scalar lst; integer dotpos; 221 lst := bfexplode0x(m,ex); 222 ex := cadr lst; 223 dotpos := caddr lst; 224 lst := car lst; 225 return bfprin!:lst (lst,ex,dotpos) 226 end; 227 228symbolic procedure bfexplode0 u; 229 % returns a list (lst ex dotpos) where 230 % lst = list of characters in mantissa 231 % (ie optional sign and digits) 232 % ex = decimal exponent 233 % dotpos = position of decimal point in lst 234 % (note that the sign is counted) 235 begin scalar r; integer m, ex; 236 r := round!:dec1 (u,if !:print!-prec!: then !:print!-prec!: 237 else !:prec!: - 2); 238 m := car r; ex := cdr r; 239 return bfexplode0x (m, ex) 240 end; 241 242 243symbolic procedure bfexplode0x (m, ex); 244 begin scalar lst, s; integer dotpos, l; 245 if m<0 then <<s := t; m := -m>>; 246 lst := explode m; 247 l := length lst; 248 if ex neq 0 and (l+ex < -!:lower!-sci!: or l+ex > !:upper!-sci!:) 249 then <<dotpos := 1; 250 ex := ex + l - 1>> 251 else <<dotpos := l + ex; 252 ex := 0; 253 if dotpos > l - 1 254 % 255 % add dotpos - l + 1 zeroes at the end 256 % 257 then lst := nconc!*(lst,nlist('!0,dotpos - l + 1)) 258 else while dotpos < 1 do <<lst := '!0 . lst; 259 dotpos := dotpos + 1>>>>; 260 if s then <<lst := '!- . lst; dotpos := dotpos + 1>>; 261 if null !*fullprec 262 then <<lst := reversip lst; 263 while lst and car lst eq '!0 and length lst > dotpos + 1 do 264 lst := cdr lst; 265 lst := reversip lst>>; 266 return {lst, ex, dotpos} 267 end; 268 269symbolic procedure bfprin!:lst (lst, ex, dotpos); 270 begin scalar result,ee,w; integer j; 271 ee:='e; 272 if !*fort and liter(w:=reval fort_exponent) then ee:=w else w:=nil; 273 if car lst eq '!- and dotpos = 1 then <<dotpos := 2; ex := ex - 1>>; 274 if ex neq 0 then if car lst eq '!- then <<ex := ex + dotpos - 2; 275 dotpos := 2>> 276 else <<ex := ex + dotpos - 1; dotpos := 1>> 277 else if dotpos = length lst then dotpos := -1; 278 for each char in lst do << 279 result := char . result; j := j + 1; dotpos := dotpos - 1; 280 if j=5 then <<if !*nat and !*bfspace then result := '! . result; 281 j := 0>>; 282 if dotpos = 0 then <<result := '!. . result; j := j + 1>>; 283 if j=5 then <<if !*nat and !*bfspace then result := '! . result; 284 j := 0>>>>; 285 if ex neq 0 or w then << 286 if not (!*nat and !*bfspace) then result := ee . result 287 else if j=0 then <<result := '! . ee . result; j := 2>> 288 else if j=1 then <<result := '! . ee . '! . result; j := 4>> 289 else if j=2 290 then <<result := '! . '! . ee . '! . result; j := 0>> 291 else if j=3 then <<result := '! . ee . '! . result; j := 0>> 292 else if j=4 then <<result := '! . ee . '! . result; j := 2>>; 293 lst := if ex > 0 then '!+ . explode ex else explode ex; 294 for each char in lst do << 295 result := char . result; j := j + 1; 296 if j=5 then <<if !*nat and !*bfspace then result := '! . result; 297 j := 0>>>>>>; 298 % if !*nat then for each char in reversip result do prin2!* char else 299 prin2!* smallcompress reversip result 300 end; 301 302symbolic procedure smallcompress (li); 303 begin scalar ll; 304 if (ll := length li)>1000 305 then <<li := smallsplit(li,ll/2); 306 return concat(smallcompress car li,smallcompress cdr li)>> 307 else return compress ('!" . reversip('!" . reversip li)) 308 end; 309 310symbolic procedure smallsplit (li,ln); 311 begin scalar aa; 312 for i:=1:ln do <<aa := car li . aa; li := rest li>>; 313 return (reversip aa) . li 314 end; 315 316symbolic procedure scientific_notation n; 317 begin scalar oldu,oldl; 318 oldu := !:upper!-sci!:; oldl := !:lower!-sci!: + 1; 319 if fixp n 320 then << 321 if n<0 322 then rerror(arith,1, 323 {"Invalid argument to scientific_notation:",n}); 324 !:lower!-sci!: := n - 1; !:upper!-sci!: := n; 325 >> 326 else if eqcar(n,'list) and length n=3 327 then if not (fixp cadr n and fixp caddr n) 328 then rerror(arith,2, 329 {"Invalid argument to scientific_notation:",n}) 330 else <<!:upper!-sci!: := cadr n; 331 !:lower!-sci!: := caddr n - 1 >>; 332 return {'list,oldu,oldl} % Return previous range. 333 end; 334 335flag('(scientific_notation), 'opfn); 336 337symbolic procedure order!: nmbr; 338 % This function counts the order of a number "n". NMBR is a bigfloat 339 % representation of "n". 340 % **** ORDER(n)=k if 2**k <= ABS(n) < 2**(k+1) 341 % **** when n is not 0, and ORDER(0)=0. 342 % 343 if mt!: nmbr = 0 then 0 else preci!: nmbr + ep!: nmbr - 1; 344 345symbolic inline procedure decprec!:(nmbr, k); 346 make!:ibf(ashift(mt!: nmbr,-k), ep!: nmbr + k); 347 348symbolic inline procedure incprec!:(nmbr, k); 349 make!:ibf(ashift(mt!: nmbr,k), ep!: nmbr - k); 350 351symbolic procedure conv!:mt(nmbr, k); 352 % This function converts a number "n" to an equivalent number of 353 % binary precision K by rounding "n" or adding "0"s to "n". 354 % NMBR is a binary bigfloat representation of "n". 355 % K is a positive integer. 356 if bfp!: nmbr and fixp k and k > 0 357 then if (k := preci!: nmbr - k) = 0 then nmbr 358 else if k < 0 then incprec!:(nmbr, -k) 359 else round!:last(decprec!:(nmbr, k - 1)) 360 else bflerrmsg 'conv!:mt; 361 362symbolic procedure round!:mt(nmbr, k); 363 % This function rounds a number "n" at the (K+1)th place and returns 364 % an equivalent number of binary precision K if the precision of "n" 365 % is greater than K, else it returns the given number unchanged. 366 % NMBR is a bigfloat representation of "n". K is a positive integer. 367 if bfp!: nmbr and fixp k and k > 0 then 368 if (k := preci!: nmbr - k - 1) < 0 then nmbr 369 else if k = 0 then round!:last nmbr 370 else round!:last decprec!:(nmbr, k) 371 else bflerrmsg 'round!:mt; 372 373symbolic procedure round!:ep(nmbr, k); 374% This function rounds a number "n" and returns an 375% equivalent number having the exponent K if 376% the exponent of "n" is less than K, else 377% it returns the given number unchanged. 378% NMBR is a BINARY BIG-FLOAT representation of "n". 379% K is an integer (positive or negative). 380 if bfp!: nmbr and fixp k then 381 if (k := k - 1 - ep!: nmbr) < 0 then nmbr 382 else if k = 0 then round!:last nmbr 383 else round!:last decprec!:(nmbr, k) 384 else bflerrmsg 'round!:ep$ 385 386symbolic procedure round!:last nmbr; 387 % This function rounds a number "n" at its last place. 388 % NMBR is a binary bigfloat representation of "n". 389 << if m < 0 then << m := -m; s := t >>; 390 m := if oddintp m then lshift (m, -1) + 1 391 else lshift (m, -1); 392 if s then m := -m; 393 make!:ibf (m, e) >> 394 where m := mt!: nmbr, e := ep!: nmbr + 1, s := nil; 395 396symbolic procedure cut!:mt(nmbr,k); 397% This function returns a given number "n" unchanged 398% if its binary precision is not greater than K, else it 399% cuts off its mantissa at the (K+1)th place and 400% returns an equivalent number of precision K. 401% **** CAUTION! No rounding is made. 402% NMBR is a BINARY BIG-FLOAT representation of "n". 403% K is a positive integer. 404 if bfp!: nmbr and fixp k and k > 0 then 405 if (k := preci!: nmbr - k) <= 0 then nmbr 406 else decprec!:(nmbr, k) 407 else bflerrmsg 'cut!:mt$ 408 409symbolic procedure cut!:ep(nmbr, k); 410% This function returns a given number "n" unchanged 411% if its exponent is not less than K, else it 412% cuts off its mantissa and returns an equivalent 413% number of exponent K. 414% **** CAUTION! No rounding is made. 415% NMBR is a BINARY BIG-FLOAT representation of "n". 416% K is an integer (positive or negative). 417 if bfp!: nmbr and fixp k then 418 if (k := k - ep!: nmbr) <= 0 then nmbr 419 else decprec!:(nmbr, k) 420 else bflerrmsg 'cut!:ep$ 421 422symbolic procedure bftrim!: v; 423 normbf round!:mt(v,!:bprec!: - 3); 424 425symbolic procedure decimal2internal (base10,exp10); 426 if exp10 >= 0 then i2bf!: (base10 * 10**exp10) 427 else divide!-by!-power!-of!-ten (i2bf!: base10, -exp10); 428 429symbolic procedure read!:num(n); 430 % This function reads a number or a number-like entity N 431 % and constructs a bigfloat representation of it. 432 % N is an integer, a floating-point number, or a string 433 % representing a number. 434 % **** If the system does not accept or may incorrectly 435 % **** accept the floating-point numbers, you can 436 % **** input them as strings such as "1.234E-56", 437 % **** "-78.90 D+12" , "+3456 B -78", or "901/234". 438 % **** A rational number in a string form is converted 439 % **** to a bigfloat of precision !:PREC!: if 440 % **** !:PREC!: is not NIL, else the precision of 441 % **** the result is set 170. 442 % **** Some systems set the maximum size of strings. If 443 % **** you want to input long numbers exceeding 444 % **** such a maximum size, please use READ!:LNUM. 445 if fixp n then make!:ibf(n, 0) 446 else if not(numberp n or stringp n) then bflerrmsg 'read!:num 447 else begin integer j,m,sign; scalar ch,u,v,l,appear!.,appear!/; 448 j := m := 0; 449 sign := 1; 450 u := v := appear!. := appear!/ := nil; 451 l := explode n; 452 loop: ch := car l; 453 if digit ch then << u := ch . u; j := j + 1 >> 454 else if ch eq '!. then << appear!. := t; j := 0 >> 455 else if ch eq '!/ then << appear!/ := t; v := u; u := nil >> 456 else if ch eq '!- then sign := -1 457 else if ch memq '(!E !D !B !e !d !b) then go to jump; %JBM 458 if l := cdr l then goto loop else goto make; 459 jump: while l := cdr l do 460 <<if digit(ch := car l) or ch eq '!- 461 then v := ch . v >>; 462 l := reverse v; 463 % Was erroneously smallcompress. 464 if car l eq '!- then m := - compress cdr l 465 else m:= compress l; 466 make: u := reverse u; 467 v := reverse v; 468 if appear!/ then 469 return conv!:r2bf(make!:ratnum(sign*compress v,compress u), 470 if !:bprec!: then !:bprec!: else 170); 471 if appear!. then j := - j else j := 0; 472 if sign = 1 then u := compress u else u := - compress u; 473 return round!:mt (decimal2internal (u, j + m), !:bprec!:) 474 where !:bprec!: := if !:bprec!: then !:bprec!: 475 else msd!: abs u 476 end; 477 478symbolic procedure abs!: nmbr; 479 % This function makes the absolute value of "n". N is a binary 480 % bigfloat representation of "n". 481 if mt!: nmbr > 0 then nmbr else make!:ibf(- mt!: nmbr, ep!: nmbr); 482 483symbolic procedure minus!: nmbr; 484 % This function makes the minus number of "n". N is a binary 485 % bigfloat representation of "n". 486 make!:ibf(- mt!: nmbr, ep!: nmbr); 487 488symbolic procedure plus!:(n1,n2); 489 begin scalar m1,m2,e1,e2,d; return 490 if (m1 := mt!: n1)=0 then n2 491 else if (m2 := mt!: n2)=0 then n1 492 else if (d := (e1 := ep!: n1)-(e2 := ep!: n2))=0 493 then make!:ibf(m1+m2, e1) 494 else if d>0 then make!:ibf(ashift(m1,d)+m2,e2) 495 else make!:ibf(m1+ashift(m2,-d),e1) end; 496 497symbolic procedure difference!:(n1,n2); 498 begin scalar m1,m2,e1,e2,d; return 499 if (m1 := mt!: n1)=0 then minus!: n2 500 else if (m2 := mt!: n2)=0 then n1 501 else if (d := (e1 := ep!: n1)-(e2 := ep!: n2))=0 502 then make!:ibf(m1 - m2, e1) 503 else if d>0 then make!:ibf(ashift(m1,d) - m2,e2) 504 else make!:ibf(m1 - ashift(m2,-d),e1) end; 505 506symbolic procedure times!:(n1, n2); 507 % This function calculates the product of "n1" and "n2". 508 % N1 and N2 are bigfloat representations of "n1" and "n2". 509 make!:ibf(mt!: n1 * mt!: n2, ep!: n1 + ep!: n2); 510 511symbolic procedure divide!:(n1,n2,k); 512 % This function calculates the quotient of "n1" and "n2", with the 513 % precision K, by rounding the ratio of "n1" and "n2" at the (K+1)th 514 % place. N1 and N2 are bigfloat representations of "n1" and "n2". 515 % K is any positive integer. 516 begin 517 n1 := conv!:mt(n1, k + preci!: n2 + 1); 518 n1 := make!:ibf(mt!: n1 / mt!: n2, ep!: n1 - ep!: n2); 519 return round!:mt(n1, k) 520 end; 521 522symbolic procedure max2!:(a,b); 523 % This function returns the larger of "n1" and "n2". 524 % N1 and N2 are bigfloat representations of "n1" and "n2". 525 if greaterp!:(a,b) then a else b; 526 527macro procedure max!: x; expand(cdr x,'max2!:); 528 529symbolic procedure min2!:(a,b); 530 % This function returns the smaller of "n1" and "n2". 531 % N1 and N2 are binary bigfloat representations of "n1" and "n2". 532 if greaterp!:(a,b) then b else a; 533 534macro procedure min!: x; expand(cdr x,'min2!:); 535 536symbolic procedure greaterp!:(a,b); 537% this function calculates the a > b, but avoids 538% generating large numbers if magnitude difference is large. 539 if mt!: a=0 then mt!: b<0 540 else if mt!: b=0 then mt!: a>0 541 else if ep!: a=ep!: b then mt!: a>mt!: b else 542 (((if d=0 then ma>mb else 543 ((if d>p2 then ma>0 else if d<-p2 then mb<0 544 else if d>0 then ashift(ma,d)>mb 545 else ma>ashift(mb,-d)) 546 where p2=2*!:bprec!:)) 547 where d=ep!: a - ep!: b, ma=mt!: a, mb=mt!: b) 548 where a= normbf a, b=normbf b); 549 550symbolic procedure equal!:(a,b); 551 %tests bfloats for a=b rapidly without generating digits. %SK 552 zerop mt!: a and zerop mt!: b or 553 ep!:(a := normbf a)=ep!:(b := normbf b) and mt!: a=mt!: b; 554 555symbolic procedure lessp!:(n1, n2); 556 % This function returns T if "n1" < "n2" else returns NIL. 557 % N1 and N2 are bigfloat representations of "n1" and "n2". 558 greaterp!:(n2, n1); 559 560symbolic procedure leq!:(n1, n2); 561 % This function returns T if "n1" <= "n2" else returns NIL. 562 % N1 and N2 are bigfloat representations of "n1" and "n2". 563 not greaterp!:(n1, n2); 564 565symbolic procedure minusp!: x; 566 % This function returns T if "x"<0 else returns NIL. 567 % X is any Lisp entity. 568 bfp!: x and mt!: x < 0; 569 570symbolic procedure make!:ratnum(nm,dn); 571 % This function constructs an internal representation 572 % of a rational number composed of the numerator 573 % NM and the denominator DN. 574 % NM and DN are any integers (positive or negative). 575 % **** Four routines in this section are temporary. 576 % **** That is, if your system has own routines 577 % **** for rational number arithmetic, you can 578 % **** accommodate our system to yours only by 579 % **** redefining these four routines. 580 if zerop dn then rerror(arith,3,"Zero divisor in make:ratnum") 581 else if dn > 0 then '!:ratnum!: . (nm . dn) 582 else '!:ratnum!: . (-nm . -dn); 583 584symbolic procedure ratnump!:(x); 585 % This function returns T if X is a rational number 586 % representation, else NIL. 587 % X is any Lisp entity. 588 eqcar(x,'!:ratnum!:); %JBM Change to EQCAR. 589 590symbolic inline procedure numr!: rnmbr; 591 % This function selects the numerator of a rational number "n". 592 % RNMBR is a rational number representation of "n". 593 cadr rnmbr; 594 595symbolic inline procedure denm!: rnmbr; 596 % This function selects the denominator of a rational number "n". 597 % RNMBR is a rational number representation of "n". 598 cddr rnmbr; 599 600symbolic procedure conv!:r2bf(rnmbr,k); 601 % This function converts a rational number RNMBR to a bigfloat of 602 % precision K, i.e., a bigfloat representation with a given 603 % precision. RNMBR is a rational number representation. K is a 604 % positive integer. 605 if ratnump!: rnmbr and fixp k and k > 0 606 then divide!:(make!:ibf( numr!: rnmbr, 0), 607 make!:ibf( denm!: rnmbr, 0),k) 608 else bflerrmsg 'conv!:r2bf; 609 610endmodule; 611 612end; 613