1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2% 3% File: PXNK:arith387.sl for Intel 386 /387 4% Description: Generic arithmetic: 387 supplementum 5% Author: Haerbert Melenk 6% Created: 14-May 1990 7% Modified: 8% Mode: Lisp 9% Package: Utilities 10% Status: Open Source: BSD License 11% 12% Redistribution and use in source and binary forms, with or without 13% modification, are permitted provided that the following conditions are met: 14% 15% * Redistributions of source code must retain the relevant copyright 16% notice, this list of conditions and the following disclaimer. 17% * Redistributions in binary form must reproduce the above copyright 18% notice, this list of conditions and the following disclaimer in the 19% documentation and/or other materials provided with the distribution. 20% 21% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 22% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 23% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 24% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 25% CONTRIBUTORS 26% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 27% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 28% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 29% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 30% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 31% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 32% POSSIBILITY OF SUCH DAMAGE. 33% 34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 35% 36% 05-03-1992 (HM) 37% version with using round-to-nearest mode for arithmetic operations 38% 39% 27-nov-1991 (HM) 40% moved to masked coprocessor mode: exceptions are tested by investigating 41% the result exponent (IEEE convention). 42% 43% the functions are mapped to the corresponding Cmacros, 44% which then are expanded to coprocessor instructions 45% 46% This module can be loaded to a PSL arithmetic which is 47% compiled avoiding the coprocessor instructions. 48 49(fluid '(fcontrolword** fcontrolwordadr** fcontrolwordn** fcontrolwordadrn**)) 50 51(setq fcontrolword** 16#0f7f) % 387 mask: truncate, all exceptions masked 52(setq fcontrolwordadr** % address of mask 53 (plus2 symval (times2 4 (id2int 'fcontrolword**)))) 54 55(setq fcontrolwordn** 16#037f) % 387 mask: round nearest, all exceptions masked 56(setq fcontrolwordadrn** % address of mask 57 (plus2 symval (times2 4 (id2int 'fcontrolwordn**)))) 58 59 60(compiletime (progn 61 (put 'clear387 'opencode 62 '((*move ($fluid fcontrolwordadr**) (reg 5)) 63 (fldcw (displacement (reg 5) 0)) 64 )) 65 (put 'clear387n 'opencode 66 '((*move ($fluid fcontrolwordadrn**) (reg 5)) 67 (fldcw (displacement (reg 5) 0)) 68 )) 69 70 (define-constant IEEEbias 1023) 71 (define-constant IEEEmask 2047) 72 (ds floathighword (x) (floathighorder (wdifference x 4))) 73 (ds floatlowword (x) (floatloworder (wdifference x 4))) 74 (ds IEEEexpt(u) (wdifference (wand IEEEmask 75 (wshift (floatHighword u) -20)) 76 IEEEbias)) 77 (ds *testNan(x y) 78 (when (weq IEEEmask (wand IEEEmask 79 (wshift (getmem (wplus2 4 x)) -20))) 80 (ftsignal y))) 81 (define-constant ZEROexp (minus IEEEbias)) 82 83)) 84 85% (fluid '(*arith387*)) 86% (setq *arith387* t) 87% (*freset) 88 89(de ftsignal(u) 90 (stderror (bldmsg "Floating point exception in %w" u))) 91 92(off r2i) 93 94(de *wfloat (x y) 95 (clear387n) 96 (*wfloat x y) 97 (*testNan x "int-float-conversion")) 98 99(de *fplus2 (x y z) 100 (clear387n) 101 (*fplus2 x y z) 102 (*testNan x "float-plus")) 103 104(de *fdifference (x y z) 105 (clear387n) 106 (*fdifference x y z) 107 (*testNan x "float-difference")) 108 109(de *ftimes2 (x y z) 110 (clear387n) 111 (*ftimes2 x y z) 112 (*testnan x "float-times")) 113 114(de *fquotient (x y z) 115 (clear387n) 116 (*fquotient x y z) 117 (*testNan x "float-quotient")) 118 119% as the Cmacro is intelligent enough, we omit the T and NIL parameters 120 121(de *fgreaterp (x y) 122 (clear387n) 123 (*fgreaterp x y )) 124 125(de *flessp (x y) 126 (clear387n) 127 (*flessp x y )) 128 129(de *wfix (f) 130 (clear387) 131 (*wfix f)) 132 133% (re)set for the 387 coprocessor status: 134% clear exception flag 135% set rounding mode to "towards zero" 136% set interrupt mask to zerodiv + overflow + invalid 137 138(on r2i) 139 140(lap '((*entry uxsin expr 1) 141 (*move (displacement (reg 1) 0) (reg 3)) 142 (*move (displacement (reg 2) 0) (reg 3)) 143 (fld (displacement (reg 2) 0)) 144 (fsin) 145 (fstp (displacement (reg 1) 0)) 146 (*linke 0 ftestresult expr 1))) 147 148 149(lap '((*entry uxcos expr 1) 150 (*move (displacement (reg 1) 0) (reg 3)) 151 (*move (displacement (reg 2) 0) (reg 3)) 152 (fld (displacement (reg 2) 0)) 153 (fcos) 154 (fstp (displacement (reg 1) 0)) 155 (*linke 0 ftestresult expr 1))) 156 157 158(lap '((*entry uxtan expr 1) 159 (*move (displacement (reg 1) 0) (reg 3)) 160 (*move (displacement (reg 2) 0) (reg 3)) 161 (fld (displacement (reg 2) 0)) 162 (fptan) 163 (fstp (displacement (reg 1) 0)) 164 (*linke 0 ftestresult expr 1))) 165 166(lap '((*entry uxatan expr 1) 167 (fld (displacement (reg 2) 0)) % the atan argument 168 (fld (displacement (reg 3) 0)) % should be a one 169 (fpatan) 170 (fstp (displacement (reg 1) 0)) 171 (*linke 0 ftestresult expr 1))) 172 173 174(lap '((*entry uuxexp expr 1) 175 (*move (displacement (reg 1) 0) (reg 3)) 176 (*move (displacement (reg 2) 0) (reg 3)) 177 (fld (displacement (reg 2) 0)) 178 (fldl2e) 179 (fmulp) 180 (f2xm1) 181 (fld1) 182 (faddp) 183 (fstp (displacement (reg 1) 0)) 184 (*linke 0 ftestresult expr 1))) 185 186(lap '((*entry uxlog expr 1) 187 (*move (displacement (reg 1) 0) (reg 3)) 188 (*move (displacement (reg 2) 0) (reg 3)) 189 (Fldln2) 190 (fld (displacement (reg 2) 0)) 191 (fyl2x) 192 (fstp (displacement (reg 1) 0)) 193 (*linke 0 ftestresult expr 1))) 194 195(flag '(factorial) 'lose) 196 197 198(lap '((*entry uxsqrt expr 1) 199 (*move (displacement (reg 1) 0) (reg 3)) 200 (*move (displacement (reg 2) 0) (reg 3)) 201 (fld (displacement (reg 2) 0)) 202 (fsqrt) 203 (fstp (displacement (reg 1) 0)) 204 (*linke 0 ftestresult expr 1))) 205 206(de ftestresult(x) (*testNan x "elem. function")) 207 208%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 209% 210% File: $pu/fast-math.sl 211% Description: Some useful mathematical functions for PSL 212% Author: Contributions from Galway, Griss, Irish, Morrison, and others 213% Created: 214% Modified: Wed Jan 7 08:16:45 1987 (Russ Fish) 215% Mode: Lisp 216% Package: Utilities 217% Status: Experimental 218% Notes: 219% Compiletime: 220% Runtime: 221% 222% (c) Copyright 1987, University of Utah, all rights reserved. 223% 224%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 225% 226% Revisions: 227% 228% Tue Oct 20 18:23:35 1987 (Jim Cobb) 229% Added mathlib-style argument range checking. 230% Wed Jan 7 08:17:05 1987 (Russ Fish) 231% Converted mathlib.sl to fast-math.sl, using C math external functions. 232% Also added if_system for VAX to copyd external- routines to the names 233% actually used in the functions. 234% Mar 11 1986, Russ Fish 235% More precise definitions of floating point constants now work. 236% 9 Jan 1985 1716-PST (Cris Perdue) 237% Installed new definition of MOD with compiletime load of 238% useful to support it. Old definition was buggy. 239% 02-Dec-83 18:07:29 (Nancy Kendzierski) 240% Translated from Rlisp to Lisp. 241% 02-Dec-83 15:54:03 Nancy Kendzierski 242% Changed 'newnam stuff to setq's of fluids, since 'newnam is an rlisp-only 243% hack. 244% 01-Dec-83 19:21:50 Nancy Kendzierski 245% Changed numbers to correct round-off problem w/un-rlisp. Since this was 246% done by hand, the process was prone to error. I don't know if I got them 247% all correct. 248% 1 July 1983, Martin Griss 249% LPRIM replaced by ERRORPRINTF so works from BARE-PSL 250% MATHLIB.RED, 16-Dec-82 21:56:52, Edit by GALWAY 251% Various fixes and enhancements too numerous for me to remember. 252% Includes fixes in SQRT function, modifications of RANDOM and other 253% functions to bring them more in line with Common Lisp, addition of MOD 254% and FLOOR. 255% <PSL.UTIL>MATHLIB.RED.13, 13-Sep-82 08:49:52, Edit by BENSON 256% Bug in EXP, changed 2**N to 2.0**N 257% <PSL.UTIL>MATHLIB.RED.12, 2-Sep-82 09:22:19, Edit by BENSON 258% Changed all calls in REDERR to calls on STDERROR 259% <PSL.UTIL>MATHLIB.RED.2, 17-Jan-82 15:48:21, Edit by GRISS 260% changed for PSL 261% 262%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 263 264%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 265% Most of these routines not very heavily tested. 266% 267% Should these names be changed so that they all begin with an F or some 268% other distinguishing mark? Are they in conflict with anything? Or should 269% we wait until we have packages? 270% Consider using Sasaki's BigFloat package -- it has all this and more, to 271% arbitrary precision. The only drawback is speed. 272% 273%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 274 275(compiletime (load useful if-system)) 276 277%***************** Constants previously declared as NewNam's ***************** 278(fluid '(number2pi numberpi numberpi!/2 numberpi!/4 number3pi!/4 279 number!-2pi number!-pi number!-pi!/2 number!-pi!/4 number!-3pi!/4 280 numbere numberinversee naturallog2 naturallog10)) 281 282(setq number2pi 6.2831853071795864770) 283(setq numberpi 3.1415926535897932385) 284(setq numberpi!/2 1.57079632679489661925) 285(setq numberpi!/4 0.785398163397448309625) 286(setq number3pi!/4 2.356194490192344928875) 287 288% It's probably a mistake to have two definitions, with and without dashes... 289(setq number!-2pi 6.2831853071795864770) 290(setq number!-pi 3.1415926535897932385) 291(setq number!-pi!/2 1.57079632679489661925) 292(setq number!-pi!/4 0.785398163397448309625) 293(setq number!-3pi!/4 2.356194490192344928875) 294 295(setq numbere 2.71828182845904523536) 296(setq numberinversee 0.367879441171442321596) % 1/e 297(setq naturallog2 0.6931471805599453094) 298(setq naturallog10 2.302585092994045684) 299 300%% Do this becuase the Vax external functions have "external-" tacked 301%% on the front. 302 303%********************* Basic functions *************************************** 304 305% Mathematically, mod(x,m) should be x-m*floor(x/m) 306(de mod (x m) 307 (let* ((pair (divide x m)) 308 (q (car pair)) 309 (r (cdr pair))) 310 % Quotient in PSL division is always truncated toward 0. We 311 % adjust Q to be the "floor" of the true quotient. 312 (if (and (lessp q 0) (not (equal r 0))) 313 (setq q (sub1 q))) 314 (difference x (times m q)))) 315 316(de floor (x) 317 % Returns the largest integer less than or equal to X. (I.e. the "greatest 318 % integer" function.) 319 (if (fixp x) 320 x 321 (prog (n) 322 (setq n (fix x)) 323 % Note the trickiness to compensate for fact that (unlike APL's 324 % "FLOOR" function) FIX truncates towards zero. 325 (return (cond ((equal x (float n)) n) 326 ((geq x 0) n) 327 (t (difference n 1))))))) 328 329(de ceiling (x) 330 % Returns the smallest integer greater than or equal to X. 331 (if (fixp x) 332 x 333 (prog (n) 334 (setq n (fix x)) 335 % Note the trickiness to compensate for fact that (unlike APL's 336 % "FLOOR" function) FIX truncates towards zero. 337 (return (cond ((equal x (float n)) n) 338 ((greaterp x 0) (plus n 1)) 339 (t n)))))) 340 341(de round (x) 342 % Rounds to the closest integer. 343 % Kind of sloppy -- it's biased when the digit causing rounding is a five, 344 345 % it's a bit weird with negative arguments, round(-2.5)= -2. 346 (if (fixp x) 347 x 348 (floor (plus x 0.5)))) 349 350%***************** Trigonometric Functions *********************************** 351 352% Trig functions are all in radians. The following few functions may be used 353 354% to convert to/from degrees, or degrees/minutes/seconds. 355(de degreestoradians (x) 356 (times x 0.0174532925199432957694)) % 2*pi/360 357 358(de radianstodegrees (x) 359 (times x 57.2957795130823208761)) % 360/(2*pi) 360 361% 360/(2*pi) 362(de radianstodms (x) 363 % Converts radians to a list of degrees, minutes, and seconds (rounded, not 364 365 % truncated, to the nearest integer). 366 (prog (degs mins) 367 (setq x (radianstodegrees x)) 368 (setq degs (fix x)) 369 (setq x (times 60 (difference x degs))) 370 (setq mins (fix x)) 371 (return (list degs mins (round (times 60 (difference x mins))))))) 372 373(de dmstoradians (degs mins sex) 374 % Converts degrees, minutes, seconds to radians. 375 (degreestoradians 376 (plus degs (quotient mins 60.0) (quotient sex 3600.0)))) 377 378(de sin (x) 379 (setq x (float x)) 380 (prog (result) 381 (setq result (gtfltn)) 382 (uxsin (floatbase result) (floatbase (fltinf x))) 383 (return (mkfltn result)))) 384 385(de cos (x) 386 (setq x (float x)) 387 (prog (result) 388 (setq result (gtfltn)) 389 (uxcos (floatbase result) (floatbase (fltinf x))) 390 (return (mkfltn result)))) 391 392(de tan (x) 393 (quotient (sin x) (cos x))) 394 395(de cot (x) 396 (quotient (cos x) (sin x))) 397 398(de sec (x) 399 (quotient 1.0 (cos x))) 400 401(de csc (x) 402 (quotient 1.0 (sin x))) 403 404(de sind (x) 405 (sin (degreestoradians x))) 406 407(de cosd (x) 408 (cos (degreestoradians x))) 409 410(de tand (x) 411 (tan (degreestoradians x))) 412 413(de cotd (x) 414 (cot (degreestoradians x))) 415 416(de secd (x) 417 (sec (degreestoradians x))) 418 419(de cscd (x) 420 (csc (degreestoradians x))) 421 422(de atan (x) 423 (setq x (float x)) 424 (prog (result) 425 (setq result (gtfltn)) 426 (uxatan (floatbase result) (floatbase (fltinf x)) (floatbase(fltinf 1.0))) 427 (return (mkfltn result)))) 428 429(de acot (x) 430 (cond ((minusp x) (if (lessp x -1.0) 431 (minus (atan (quotient -1.0 x))) 432 (plus number!-pi!/2 (atan (minus x))))) 433 ((greaterp x 1.0) (atan (quotient 1.0 x))) 434 (t (difference numberpi!/2 (atan x))))) 435 436(de asec (x) 437 (acos (quotient 1.0 x))) 438 439(de acsc (x) 440 (asin (quotient 1.0 x))) 441 442(de asind (x) 443 (radianstodegrees (asin x))) 444 445(de acosd (x) 446 (radianstodegrees (acos x))) 447 448(de atand (x) 449 (radianstodegrees (atan x))) 450 451(de acotd (x) 452 (radianstodegrees (acot x))) 453 454(de asecd (x) 455 (radianstodegrees (asec x))) 456 457(de acscd (x) 458 (radianstodegrees (acsc x))) 459 460%****************** Roots and such ******************************************* 461 462(de sqrt (x) 463 (prog (result) 464 (when (lessp x 0.00000E+000) 465 (stderror (list "SQRT given negative argument:" x))) 466 (setq result (gtfltn)) 467 (uxsqrt (floatbase result) (floatbase (fltinf (float x)))) 468 (return (mkfltn result)))) 469 470%******************** Logs and Exponentials ********************************** 471 472(de exp (x) 473 % Returns the exponential (ie, e**x) of its floatnum argument as 474 % a flonum. The argument is scaled to the interval -ln 2 to 0 475 (prog (n result) 476 (setq n (ceiling (quotient x naturallog2))) 477 (setq x (difference x (times n naturallog2))) 478 (setq result (gtfltn)) 479 (uuxexp (floatbase result) (floatbase (fltinf (float x)))) 480 (return 481 (times2 (expt 2.0 n) (mkfltn result))))) 482 483(de log (x) 484 (prog (result) 485 (when (leq x 0.00000E+000) 486 (stderror (list "LOG given non-positive argument:" x))) 487 (setq result (gtfltn)) 488 (uxlog (floatbase result) (floatbase (fltinf (float x)))) 489 (return (mkfltn result)))) 490 491(de log2 (x) 492 (quotient (log x) naturallog2)) 493 494(de log10 (x) 495 (quotient (log x) naturallog10)) 496 497%********************* Random Number Generator ******************************* 498 499% The declarations below constitute a linear, congruential 500% random number generator (see Knuth, "The Art of Computer 501% Programming: Volume 2: Seminumerical Algorithms", pp9-24). 502% With the given constants it has a period of 392931 and 503% potency 6. To have deterministic behaviour, set 504% RANDOMSEED. 505% 506% Constants are: 6 2 507% modulus: 392931 = 3 * 7 * 11 508% multiplier: 232 = 3 * 7 * 11 + 1 509% increment: 65537 is prime 510% 511% Would benefit from being recoded in SysLisp, when full word integers should 512% be used with "automatic" modular arithmetic (see Knuth). Perhaps we should 513% have a longer period version? 514% By E. Benson, W. Galway and M. Griss 515(fluid '(randomseed randommodulus)) 516 517(setq randommodulus 392931) 518 519(setq randomseed (remainder (time) randommodulus)) 520 521(de next!-random!-number () 522 % Returns a pseudo-random number between 0 and RandomModulus-1 (inclusive). 523 524 (setq randomseed 525 (remainder (plus (times 232 randomseed) 65537) randommodulus))) 526 527(de random (n) 528 % Return a pseudo-random number uniformly selected from the range 0..N-1. 529 % NOTE that this used to be called RandomMod(N). Needs to be made more 530 % compatible with Common LISP's random? 531 (fix (quotient (times (float n) (next!-random!-number)) randommodulus))) 532 533(de transfersign (s val) 534 % Transfers the sign of S to Val by returning abs(Val) if S >= 0, 535 % otherwise -abs(Val). 536 (if (geq s 0) 537 (abs val) 538 (minus (abs val)))) 539 540(de dmstodegrees (degs mins sex) 541 % Converts degrees, minutes, seconds to degrees 542 (plus degs (quotient mins 60.0) (quotient sex 3600.0))) 543 544(de degreestodms (x) 545 % Converts degrees to a list of degrees, minutes, and seconds (all integers, 546 % rounded, not truncated). 547 (prog (degs mins) 548 (setq degs (fix x)) 549 (setq x (times 60 (difference x degs))) 550 (setq mins (fix x)) 551 (return (list degs mins (round (times 60 (difference x mins))))))) 552 553