1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2% 3% File: $pu/fast-math.sl 4% Description: Some useful mathematical functions for PSL 5% Author: Contributions from Galway, Griss, Irish, Morrison, and others 6% Created: 7% Modified: Wed Jan 7 08:16:45 1987 (Russ Fish) 8% Mode: Lisp 9% Package: Utilities 10% Status: Open Source: BSD License 11% Notes: 12% Compiletime: 13% Runtime: 14% 15% (c) Copyright 1982, University of Utah 16% 17% Redistribution and use in source and binary forms, with or without 18% modification, are permitted provided that the following conditions are met: 19% 20% * Redistributions of source code must retain the relevant copyright 21% notice, this list of conditions and the following disclaimer. 22% * Redistributions in binary form must reproduce the above copyright 23% notice, this list of conditions and the following disclaimer in the 24% documentation and/or other materials provided with the distribution. 25% 26% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 27% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 28% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 29% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 30% CONTRIBUTORS 31% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 32% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 33% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 34% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 35% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 36% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 37% POSSIBILITY OF SUCH DAMAGE. 38% 39%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 40% 41% Revisions: 42% 43% Fri Mar 20 1992 (Winfried Neun) 44% Added some extra fltinf for IBM RS/600 or HP snake , ... 45% Thu Jun 27 1991 (Winfried Neun) 46% Fixed problems in sin, cos ,... with non-float argument and gc 47% see comment at sin. 48% Wed Feb 13 1991 (Winfried Neun) 49% Added support for IBM RS/6000 floats 50% Tue Oct 20 18:23:35 1987 (Jim Cobb) 51% Added mathlib-style argument range checking. 52% Wed Jan 7 08:17:05 1987 (Russ Fish) 53% Converted mathlib.sl to fast-math.sl, using C math external functions. 54% Also added if_system for VAX to copyd external- routines to the names 55% actually used in the functions. 56% Mar 11 1986, Russ Fish 57% More precise definitions of floating point constants now work. 58% 9 Jan 1985 1716-PST (Cris Perdue) 59% Installed new definition of MOD with compiletime load of 60% useful to support it. Old definition was buggy. 61% 02-Dec-83 18:07:29 (Nancy Kendzierski) 62% Translated from Rlisp to Lisp. 63% 02-Dec-83 15:54:03 Nancy Kendzierski 64% Changed 'newnam stuff to setq's of fluids, since 'newnam is an rlisp-only 65% hack. 66% 01-Dec-83 19:21:50 Nancy Kendzierski 67% Changed numbers to correct round-off problem w/un-rlisp. Since this was 68% done by hand, the process was prone to error. I don't know if I got them 69% all correct. 70% 1 July 1983, Martin Griss 71% LPRIM replaced by ERRORPRINTF so works from BARE-PSL 72% MATHLIB.RED, 16-Dec-82 21:56:52, Edit by GALWAY 73% Various fixes and enhancements too numerous for me to remember. 74% Includes fixes in SQRT function, modifications of RANDOM and other 75% functions to bring them more in line with Common Lisp, addition of MOD 76% and FLOOR. 77% <PSL.UTIL>MATHLIB.RED.13, 13-Sep-82 08:49:52, Edit by BENSON 78% Bug in EXP, changed 2**N to 2.0**N 79% <PSL.UTIL>MATHLIB.RED.12, 2-Sep-82 09:22:19, Edit by BENSON 80% Changed all calls in REDERR to calls on STDERROR 81% <PSL.UTIL>MATHLIB.RED.2, 17-Jan-82 15:48:21, Edit by GRISS 82% changed for PSL 83% 84%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 85 86%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 87% Most of these routines not very heavily tested. 88% 89% Should these names be changed so that they all begin with an F or some 90% other distinguishing mark? Are they in conflict with anything? Or should 91% we wait until we have packages? 92% Consider using Sasaki's BigFloat package -- it has all this and more, to 93% arbitrary precision. The only drawback is speed. 94% 95%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 96 97(compiletime (load useful if-system)) 98 99%***************** Constants previously declared as NewNam's ***************** 100(fluid '(number2pi numberpi numberpi!/2 numberpi!/4 number3pi!/4 101 number!-2pi number!-pi number!-pi!/2 number!-pi!/4 number!-3pi!/4 102 numbere numberinversee naturallog2 naturallog10)) 103 104(setq number2pi 6.2831853071795864770) 105(setq numberpi 3.1415926535897932385) 106(setq numberpi!/2 1.57079632679489661925) 107(setq numberpi!/4 0.785398163397448309625) 108(setq number3pi!/4 2.356194490192344928875) 109 110% It's probably a mistake to have two definitions, with and without dashes... 111(setq number!-2pi 6.2831853071795864770) 112(setq number!-pi 3.1415926535897932385) 113(setq number!-pi!/2 1.57079632679489661925) 114(setq number!-pi!/4 0.785398163397448309625) 115(setq number!-3pi!/4 2.356194490192344928875) 116 117(setq numbere 2.71828182845904523536) 118(setq numberinversee 0.367879441171442321596) % 1/e 119(setq naturallog2 0.6931471805599453094) 120(setq naturallog10 2.302585092994045684) 121 122%% Do this because the Vax external functions have "external-" tacked 123%% on the front. 124(if_system VAX 125 (progn 126 (copyd 'uxsin 'external-uxsin) 127 (copyd 'uxcos 'external-uxcos) 128 (copyd 'uxtan 'external-uxtan) 129 (copyd 'uxasin 'external-uxasin) 130 (copyd 'uxacos 'external-uxacos) 131 (copyd 'uxatan 'external-uxatan) 132 (copyd 'uxsqrt 'external-uxsqrt) 133 (copyd 'uxexp 'external-uxexp) 134 (copyd 'uxlog 'external-uxlog) 135 (copyd 'uxatan2 'external-uxatan2))) 136 137(if_system RS6000 138 (compiletime (ds fltinf (x) (mkstr x)))) % crazy , isnt it? 139 140(if_system HP-RISC 141 (compiletime (ds fltinf (x) (mkvec x)))) % even more 142 143%********************* Basic functions *************************************** 144 145% Mathematically, mod(x,m) should be x-m*floor(x/m) 146(de mod (x m) 147 (let* ((pair (divide x m)) 148 (q (car pair)) 149 (r (cdr pair))) 150 % Quotient in PSL division is always truncated toward 0. We 151 % adjust Q to be the "floor" of the true quotient. 152 (if (and (lessp q 0) (not (equal r 0))) 153 (setq q (sub1 q))) 154 (difference x (times m q)))) 155 156(de floor (x) 157 % Returns the largest integer less than or equal to X. (I.e. the "greatest 158 % integer" function.) 159 (if (fixp x) 160 x 161 (prog (n) 162 (setq n (fix x)) 163 % Note the trickiness to compensate for fact that (unlike APL's 164 % "FLOOR" function) FIX truncates towards zero. 165 (return (cond ((equal x (float n)) n) 166 ((geq x 0) n) 167 (t (difference n 1))))))) 168 169(de ceiling (x) 170 % Returns the smallest integer greater than or equal to X. 171 (if (fixp x) 172 x 173 (prog (n) 174 (setq n (fix x)) 175 % Note the trickiness to compensate for fact that (unlike APL's 176 % "FLOOR" function) FIX truncates towards zero. 177 (return (cond ((equal x (float n)) n) 178 ((greaterp x 0) (plus n 1)) 179 (t n)))))) 180 181(de round (x) 182 % Rounds to the closest integer. 183 % Kind of sloppy -- it's biased when the digit causing rounding is a five, 184 185 % it's a bit weird with negative arguments, round(-2.5)= -2. 186 (if (fixp x) 187 x 188 (floor (plus x 0.5)))) 189 190%***************** Trigonometric Functions *********************************** 191 192% Trig functions are all in radians. The following few functions may be used 193 194% to convert to/from degrees, or degrees/minutes/seconds. 195(de degreestoradians (x) 196 (times x 0.0174532925199432957694)) % 2*pi/360 197 198(de radianstodegrees (x) 199 (times x 57.2957795130823208761)) % 360/(2*pi) 200 201% 360/(2*pi) 202(de radianstodms (x) 203 % Converts radians to a list of degrees, minutes, and seconds (rounded, not 204 205 % truncated, to the nearest integer). 206 (prog (degs mins) 207 (setq x (radianstodegrees x)) 208 (setq degs (fix x)) 209 (setq x (times 60 (difference x degs))) 210 (setq mins (fix x)) 211 (return (list degs mins (round (times 60 (difference x mins))))))) 212 213(de dmstoradians (degs mins sex) 214 % Converts degrees, minutes, seconds to radians. 215 (degreestoradians 216 (plus degs (quotient mins 60.0) (quotient sex 3600.0)))) 217 218(de sin (x) 219 (prog (result) 220 (setq x (float x)) % we have to insure that no gc can happen 221 % between gtfltn and mkfltn 222 (setq result (gtfltn)) 223 (uxsin (floatbase (fltinf result)) (floatbase (fltinf x))) 224% used to be (uxsin (floatbase (fltinf result)) (floatbase (fltinf (float x)))) 225 (return (mkfltn result)))) 226 227(de cos (x) 228 (prog (result) 229 (setq x (float x)) 230 (setq result (gtfltn)) 231 (uxcos (floatbase (fltinf result)) (floatbase (fltinf x))) 232 (return (mkfltn result)))) 233 234(de tan (x) 235 (prog (result) 236 (setq x (float x)) 237 (setq result (gtfltn)) 238 (uxtan (floatbase (fltinf result)) (floatbase (fltinf x))) 239 (return (mkfltn result)))) 240 241(de cot (x) 242 (quotient 1.0 (tan x))) 243 244(de sec (x) 245 (quotient 1.0 (cos x))) 246 247(de csc (x) 248 (quotient 1.0 (sin x))) 249 250(de sind (x) 251 (sin (degreestoradians x))) 252 253(de cosd (x) 254 (cos (degreestoradians x))) 255 256(de tand (x) 257 (tan (degreestoradians x))) 258 259(de cotd (x) 260 (cot (degreestoradians x))) 261 262(de secd (x) 263 (sec (degreestoradians x))) 264 265(de cscd (x) 266 (csc (degreestoradians x))) 267 268(de asin (x) 269 (prog (result) 270 (when (greaterp (abs x) 1.0) 271 (stderror (list "Argument to ASIN has magnitude too large:" x))) 272 (setq x (float x)) 273 (setq result (gtfltn)) 274 (uxasin (floatbase (fltinf result)) (floatbase (fltinf x))) 275 (return (mkfltn result)))) 276 277(de acos (x) 278 (prog (result) 279 (when (greaterp (abs x) 1.0) 280 (stderror (list "Argument to ACOS has magnitude too large:" x))) 281 (setq x (float x)) 282 (setq result (gtfltn)) 283 (uxacos (floatbase (fltinf result)) (floatbase (fltinf x))) 284 (return (mkfltn result)))) 285 286(de atan (x) 287 (prog (result) 288 (setq x (float x)) 289 (setq result (gtfltn)) 290 (uxatan (floatbase (fltinf result)) (floatbase (fltinf x))) 291 (return (mkfltn result)))) 292 293(de acot (x) 294 (cond ((minusp x) (if (lessp x -1.0) 295 (minus (atan (quotient -1.0 x))) 296 (plus number!-pi!/2 (atan (minus x))))) 297 ((greaterp x 1.0) (atan (quotient 1.0 x))) 298 (t (difference numberpi!/2 (atan x))))) 299 300(de asec (x) 301 (acos (quotient 1.0 x))) 302 303(de acsc (x) 304 (asin (quotient 1.0 x))) 305 306(de asind (x) 307 (radianstodegrees (asin x))) 308 309(de acosd (x) 310 (radianstodegrees (acos x))) 311 312(de atand (x) 313 (radianstodegrees (atan x))) 314 315(de acotd (x) 316 (radianstodegrees (acot x))) 317 318(de asecd (x) 319 (radianstodegrees (asec x))) 320 321(de acscd (x) 322 (radianstodegrees (acsc x))) 323 324%****************** Roots and such ******************************************* 325 326(de sqrt (x) 327 (prog (result) 328 (when (lessp x 0.00000E+000) 329 (stderror (list "SQRT given negative argument:" x))) 330 (setq x (float x)) 331 (setq result (gtfltn)) 332 (uxsqrt (floatbase (fltinf result)) (floatbase (fltinf x))) 333 (return (mkfltn result)))) 334 335%******************** Logs and Exponentials ********************************** 336 337(de exp (x) 338 (prog (result) 339 (setq x (float x)) 340 (setq result (gtfltn)) 341 (uxexp (floatbase (fltinf result)) (floatbase (fltinf x))) 342 (return (mkfltn result)))) 343 344(de log (x) 345 (prog (result) 346 (when (leq x 0.00000E+000) 347 (stderror (list "LOG given non-positive argument:" x))) 348 (setq x (float x)) 349 (setq result (gtfltn)) 350 (uxlog (floatbase (fltinf result)) (floatbase (fltinf x))) 351 (return (mkfltn result)))) 352 353(de log2 (x) 354 (quotient (log x) naturallog2)) 355 356(de log10 (x) 357 (quotient (log x) naturallog10)) 358 359%********************* Random Number Generator ******************************* 360 361% The declarations below constitute a linear, congruential 362% random number generator (see Knuth, "The Art of Computer 363% Programming: Volume 2: Seminumerical Algorithms", pp9-24). 364% With the given constants it has a period of 392931 and 365% potency 6. To have deterministic behaviour, set 366% RANDOMSEED. 367% 368% Constants are: 6 2 369% modulus: 392931 = 3 * 7 * 11 370% multiplier: 232 = 3 * 7 * 11 + 1 371% increment: 65537 is prime ` 372% 373% Would benefit from being recoded in SysLisp, when full word integers should 374% be used with "automatic" modular arithmetic (see Knuth). Perhaps we should 375% have a longer period version? 376% By E. Benson, W. Galway and M. Griss 377(fluid '(randomseed randommodulus)) 378 379(setq randommodulus 392931) 380 381(setq randomseed (remainder (time) randommodulus)) 382 383(de next!-random!-number () 384 % Returns a pseudo-random number between 0 and RandomModulus-1 (inclusive). 385 386 (setq randomseed 387 (remainder (plus (times 232 randomseed) 65537) randommodulus))) 388 389(de random (n) 390 % Return a pseudo-random number uniformly selected from the range 0..N-1. 391 % NOTE that this used to be called RandomMod(N). Needs to be made more 392 % compatible with Common LISP's random? 393 (fix (quotient (times (float n) (next!-random!-number)) randommodulus))) 394 395(de factorial (n) 396 % Simple factorial 397 (prog (m) 398 (setq m 1) 399 (for (from i 1 n 1) (do (setq m (times m i)))) 400 (return m))) 401 402% Some functions from ALPHA_1 users 403(de atan2d (y x) 404 (radianstodegrees (atan2 y x))) 405 406(de atan2 (y x) 407 (prog (result) 408 (setq x (float x)) 409 (setq y (float y)) 410 (setq result (gtfltn)) 411 (uxatan2 (floatbase (fltinf result)) 412 (floatbase (fltinf y)) 413 (floatbase (fltinf x))) 414 (return (mkfltn result)))) 415 416(de transfersign (s val) 417 % Transfers the sign of S to Val by returning abs(Val) if S >= 0, 418 % otherwise -abs(Val). 419 (if (geq s 0) 420 (abs val) 421 (minus (abs val)))) 422 423(de dmstodegrees (degs mins sex) 424 % Converts degrees, minutes, seconds to degrees 425 (plus degs (quotient mins 60.0) (quotient sex 3600.0))) 426 427(de degreestodms (x) 428 % Converts degrees to a list of degrees, minutes, and seconds (all integers, 429 % rounded, not truncated). 430 (prog (degs mins) 431 (setq degs (fix x)) 432 (setq x (times 60 (difference x degs))) 433 (setq mins (fix x)) 434 (return (list degs mins (round (times 60 (difference x mins))))))) 435 436