(FILECREATED "17-May-84 18:11:40" {PHYLUM}<LISPCORE>SOURCES>AARITH.;9 14550
changes to: (FNS ANTILOG LOG)
previous date: " 1-May-84 21:45:24" {PHYLUM}<LISPCORE>SOURCES>AARITH.;7)
(* Copyright (c) 1981, 1983, 1984 by Xerox Corporation)
(PRETTYCOMPRINT AARITHCOMS)
(RPAQQ AARITHCOMS ((FNS LOG ANTILOG SIN ARCSIN COS ARCCOS TAN ARCTAN ARCTAN2 ATAN FEXPT)
(DECLARE: EVAL@COMPILE DONTCOPY (MACROS HORNERIFY FLEQ FGEQ)
(CONSTANTS (\EXPONENT.BIAS 127)
(2PI 6.283185)
(PI 3.141593)
(-PI -3.141593)
(-PI/2 -1.570796)
(PI/2 1.570796)
(4/PI 1.273239)
(3PI/2 4.712389)
(PI/4 .7853982)
(-PI/4 -.7853982)
(PI/180 .01745329)
(180/PI 57.29578)
(-PI/2 -1.570796)
(LN2 .6931472)
(2↑-126 1.175494E-38)))))
(DEFINEQ
(LOG
(LAMBDA (X) (* JonL "17-May-84 17:26")
(PROG ((SX (OR (FLOATP X)
(FLOAT X)))
(EXP 0)
SSUM)
(if (NOT (FGREATERP SX 0.0))
then (ERROR "LOG OF NON-POSITIVE NUMBER:" X))
(if (EQ 0 (fetch (FLOATP EXPONENT) of SX))
then
(* * Don't really need to consider unnormalized numbers, but there is a bug in Interlisp-D's floating point
arithmetic as of 3/17/84 regarding zero exponent.)
(SETQ EXP (while (FLESSP SX 2↑-126) count (SETQ SX (FTIMES SX 2.0)))))
(if (EQ SX X)
then (* Need smashable copy)
(SETQ SX (\FLOAT.BOX X)))
(SETQ EXP (IDIFFERENCE (IDIFFERENCE (fetch (FLOATP EXPONENT) of SX)
\EXPONENT.BIAS)
EXP))
(* * Depends on Interlisp-D's use of IEEE 32 bit float format internally and smashes the number to the range 1 to
2 and saves the exponent)
(replace (FLOATP EXPONENT) of SX with \EXPONENT.BIAS)
(SETQ SX (FDIFFERENCE SX 1.0))
(SETQ SSUM (HORNERIFY -.006453544 SX .03608849 -.0953294 .1676541 -.2407338 .331799
-.4998741 .9999964 0.0))
(* * Polynomial from Handbook of Mathematical Functions (edited by Aramowitz) page 69 accuracy 28 bits
(of the 24 available!))
(RETURN (FPLUS SSUM (FTIMES LN2 EXP))))))
(ANTILOG
(LAMBDA (X) (* JonL "17-May-84 18:10")
(PROG ((XX (OR (FLOATP X)
(FLOAT X)))
FRAC SSUM EXP)
(SETQ FRAC (ABS (FREMAINDER XX .6931472)))
(SETQ SSUM (HORNERIFY -.0001413161 FRAC .001329882 -.00830136 .04165735 -.1666653 .4999999
-1.0 1.0))
(* * Polynomial from Handbook of Mathematical Functions (edited by Aramowitz) page 71 accuracy 32 bits of the 24
available!)
(if (FGREATERP XX 0.0)
then (SETQ SSUM (FQUOTIENT 1.0 SSUM))) (* SSUM is in the range .5 to 2
(and series produced .5 to 1))
(SETQ EXP (fetch (FLOATP EXPONENT) of SSUM)) (* 1.442695 is (QUOTIENT 1 (LOG 2)))
(RETURN (if (AND (ILESSP 0 (add EXP (FIX (FTIMES XX 1.442695))))
(ILEQ EXP \MAX.EXPONENT))
then (* When it doesn't cause overflow of the exponent field,
we just adjust according to the INTEGERLENGTH of XX)
(replace (FLOATP EXPONENT) of SSUM with EXP)
SSUM
else (* Foo, must call \MAKEFLOAT, in order to get the
checking for \OVERFLOW and \UNDERFLOW)
(\MAKEFLOAT (fetch (FLOATP SIGNBIT) of SSUM)
(IPLUS EXP (CONSTANT (IDIFFERENCE BITSPERWORD (INTEGERLENGTH
\HIDDENBIT))))
(LOGOR \HIDDENBIT (fetch (FLOATP HIFRACTION) of SSUM))
(fetch (FLOATP LOFRACTION) of SSUM)
T SSUM))))))
(SIN
(LAMBDA (X RADIANSFLG) (* JonL " 1-May-84 21:36")
(PROG ((XX (OR (FLOATP X)
(FLOAT X)))
X2)
(if RADIANSFLG
then (if (OR (FGEQ XX 2PI)
(FLEQ XX (CONSTANT (MINUS 2PI))))
then (SETQ XX (FREMAINDER XX 2PI)))
else (if (OR (FGEQ XX 360.0)
(FLEQ XX -360.0))
then (SETQ XX (FREMAINDER XX 360.0)))
(SETQ XX (FTIMES PI/180 XX)))
(if (FLESSP XX -PI/2)
then (SETQ XX (FPLUS XX 2PI)))
(if (FGREATERP XX 3PI/2)
then (SETQ XX (FDIFFERENCE XX 2PI))
elseif (FGREATERP XX PI/2)
then (SETQ XX (FDIFFERENCE PI XX))) (* Range-reduce to between 0 and PI/2)
(RETURN (if (FGEQ XX PI/4)
then (SETQ X2 (FTIMES XX XX)) (* Polynomial from Handbook of Mathematical Functions
(edited by Aramowitz) page 76 accuracy 26 bits
(of the 24 available!))
(FTIMES XX (HORNERIFY -2.39E-8 X2 2.7526E-6 -.000198409 .008333332
-.1666667 1.0))
else (SETQ XX (FTIMES 4/PI XX))
(SETQ X2 (FTIMES XX XX)) (* Chebyshev approximation from Computer Evaluation of
Mathematical Functions (by C.
T. Fike) Page 117)
(FTIMES XX (HORNERIFY -.0000359544 X2 .002490007 -.08074545 .7853982)))))))
(ARCSIN
(LAMBDA (X RADIANSFLG) (* JonL "30-Mar-84 23:59")
(PROG ((XX (OR (FLOATP X)
(FLOAT X)))
SSUM NEGP REDUCED Z Q1 Q2)
(if (OR (FLESSP XX -1.0)
(FGREATERP XX 1.0))
then (ERROR "ARCSIN: arg not in range" XX)
elseif (FLESSP XX 0.0)
then (SETQ NEGP T)
(SETQ XX (FDIFFERENCE 0.0 XX)))
(if (FGREATERP XX .5)
then (SETQ REDUCED T)
(SETQ XX (SQRT (FTIMES .5 (FDIFFERENCE 1.0 XX)))))
(* Special case for small magnitude arguments, from
Computer Evaluation of Mathematical Funcitons
(by C. T. Fike) page 57)
(SETQ Z (FTIMES XX XX))
(SETQ Q1 (FTIMES .5315066 Z))
(SETQ Q2 (FTIMES (SETQ Q2 (FDIFFERENCE Q1 .08982446))
Q2))
(SETQ Q2 (FPLUS .3697723 (FTIMES Q2 (FPLUS .4918762 Q2))))
(SETQ SSUM (FTIMES XX (FPLUS .7533057 (FTIMES Q2 (FPLUS .6599526 Q1)))))
(if REDUCED
then (SETQ SSUM (FDIFFERENCE PI/2 (FTIMES 2.0 SSUM))))
(if NEGP
then (SETQ SSUM (FDIFFERENCE 0.0 SSUM)))
(RETURN (if RADIANSFLG
then SSUM
else (FTIMES SSUM 180/PI))))))
(COS
(LAMBDA (X RADIANSFLG) (* JonL "17-Mar-84 19:38")
(SIN (FDIFFERENCE (if RADIANSFLG
then PI/2
else 90.0)
(OR (FLOATP X)
(FLOAT X)))
RADIANSFLG)))
(ARCCOS
(LAMBDA (X RADIANSFLG) (* JonL "30-Mar-84 20:21")
(PROG ((XX (OR (FLOATP X)
(FLOAT X))))
(RETURN (FDIFFERENCE (if RADIANSFLG
then PI/2
else 90.0)
(ARCSIN XX RADIANSFLG))))))
(TAN
(LAMBDA (X RADIANSFLG) (* JonL "17-Mar-84 19:16")
(PROG ((XX (OR (FLOATP X)
(FLOAT X)))
FLIPPED Y X2)
(SETQ XX (if RADIANSFLG
then (FREMAINDER XX PI)
else (FTIMES (FREMAINDER XX 180.0)
PI/180))) (* First, normalize to between -PI and PI)
(if (FGREATERP XX PI/2)
then (SETQ XX (FDIFFERENCE XX PI))
elseif (FLESSP XX -PI/2)
then (SETQ XX (FPLUS XX PI))) (* Then normalize to between -PI/2 and PI/2)
(SETQ Y (if (FGREATERP XX PI/4)
then (SETQ FLIPPED T)
(FDIFFERENCE PI/2 XX)
elseif (FLESSP XX -PI/4)
then (SETQ FLIPPED T)
(FDIFFERENCE -PI/2 XX)
else XX))
(SETQ X2 (FTIMES Y Y))
(SETQ Y (FTIMES Y (HORNERIFY .00951681 X2 .002900525 .02456509 .05337406 .1333924 .3333314
1.0))) (* POLYNOMIAL APPROXIMATION FROM HANDBOOK OF
MATHEMATICAL FUNCTIONS (EDITED BY ABRAMOWITZ) PAGE 76
GOOD TO ALMOST 26 BITS)
(RETURN (if FLIPPED
then (FQUOTIENT 1.0 Y)
else Y)))))
(ARCTAN
(LAMBDA (X RADIANSFLG) (* JonL "17-Mar-84 17:07")
(PROG ((XX (FPLUS X 0.0))
(SSUM .002866226)
X2 FLIPPED) (* POLYNOMIAL FROM HANDBOOK OF MATHEMATICAL FUNCTIONS
(EDITED BY ARAMOWITZ) PAGE 81 ACCURACY 28 BITS)
(if (OR (FLESSP XX -1.0)
(FGREATERP XX 1.0))
then (SETQ FLIPPED (if (FLESSP XX 0.0)
then -PI/2
else PI/2))
(SETQ XX (FQUOTIENT 1.0 XX)))
(SETQ X2 (FTIMES XX XX))
(SETQ SSUM (FTIMES XX (HORNERIFY .002866226 X2 -.01616574 .04290961 -.07528964 .1065626
-.142089 .1999355 -.3333315 1.0)))
(if FLIPPED
then (SETQ SSUM (FDIFFERENCE FLIPPED SSUM)))
(RETURN (if RADIANSFLG
then SSUM
else (FTIMES SSUM 180/PI))))))
(ARCTAN2
(LAMBDA (Y X RADIANSFLG) (* JonL "17-Mar-84 21:41")
(OR (FLOATP Y)
(SETQ Y (FLOAT Y)))
(OR (FLOATP X)
(SETQ X (FLOAT X)))
(PROG ((ANGLE (ARCTAN (ABS (FQUOTIENT Y X))
T)))
(SETQ ANGLE (if (FLESSP X 0.0)
then (if (FLESSP Y 0.0)
then (* Quadrant 3)
(FPLUS -PI ANGLE)
else (* Quadrant 2)
(FDIFFERENCE PI ANGLE))
else (if (FLESSP Y 0.0)
then (* Quadrant 4)
(FDIFFERENCE 0.0 ANGLE)
else (* Quadrant 1)
ANGLE)))
(RETURN (if RADIANSFLG
then ANGLE
else (FTIMES ANGLE 180/PI))))))
(ATAN
(LAMBDA (Y X RADIANSFLG) (* JonL "17-Mar-84 23:30")
(* version of arctan which returns value in radians
between 0 and 2 PI. Copied from the PDP-10 MacLisp
machine language code.)
(OR (FLOATP Y)
(SETQ Y (FLOAT Y)))
(OR (FLOATP X)
(SETQ X (FLOAT X)))
(PROG ((Y.NEGP (FLESSP Y 0.0))
(X.NEGP (FLESSP X 0.0))
ATAN.Y ATAN.X T. TT D R (ANS -.004054058))
(SETQ ATAN.Y (if Y.NEGP
then (FDIFFERENCE 0.0 Y)
else Y))
(SETQ ATAN.X (if X.NEGP
then (FDIFFERENCE 0.0 X)
else X))
(SETQ T. (FQUOTIENT (FDIFFERENCE ATAN.Y ATAN.X)
(FPLUS ATAN.Y ATAN.X)))
(SETQ R (FTIMES T. T.))
(SETQ D (FTIMES T. (HORNERIFY -.004054058 R .02186123 -.05590989 .09642004 -.1390853
.1994654 -.3332986 .9999994)))
(SETQ TT (if (FLESSP D 0.0)
then (FDIFFERENCE 0.0 D)
else D))
(SETQ D (if (OR (FGEQ TT .7855)
(FLESSP TT .7853))
then (FPLUS D .7853982)
elseif (FLESSP D 0.0)
then (* When the rational approximation is not very good, we
can patch it up by using Y/X and an approximation for
(ARCTAN Y/X))
(FQUOTIENT ATAN.Y ATAN.X)
else (* Corresponds to label ATAN.2)
(FPLUS PI/2 (FQUOTIENT (FDIFFERENCE 0.0 ATAN.X)
ATAN.Y))))
ATAN.4 (* We now have a quadrant-1 result;
patch it up to get other quadrant values.)
(SETQ D (if X.NEGP
then (if Y.NEGP
then (* Quadrant 3)
(FPLUS PI D)
else (* Quadrant 2)
(FDIFFERENCE PI D))
else (if Y.NEGP
then (* Quadrant 4)
(FDIFFERENCE 2PI D)
else (* Quadrant 1)
D)))
(RETURN (if RADIANSFLG
then D
else (FTIMES D 180/PI))))))
(FEXPT
(LAMBDA (A N) (* JonL "17-Mar-84 21:57")
(* In addition to coercing the args to floating-point,
this handles the case of negative values for N)
(ANTILOG (FTIMES (OR (FLOATP N)
(FLOAT N))
(LOG (OR (FLOATP A)
(FLOAT A)))))))
)
(DECLARE: EVAL@COMPILE DONTCOPY
(DECLARE: EVAL@COMPILE
(PUTPROPS HORNERIFY MACRO (X
(PROG ((INITIAL (CAR X))
(VARNAME (CADR X))
(COEFFICIENTS (CDDR X))
TERM)
(OR COEFFICIENTS (SHOULDNT))
(OR (AND (LITATOM VARNAME)
VARNAME
(NEQ T VARNAME))
(\ILLEGAL.ARG VARNAME))
(SETQ TERM (LIST (QUOTE FPLUS)
(LIST (QUOTE FTIMES)
INITIAL VARNAME)
(CAR COEFFICIENTS)))
(OR (CONSTANTEXPRESSIONP (CAR COEFFICIENTS))
(ARGS.COMMUTABLEP (CAR COEFFICIENTS)
(CADR TERM))
(LISPERROR X "Can't hack non-commutable coefficient expressions"))
(RETURN (COND
((NULL (CDR COEFFICIENTS))
TERM)
(T (CONS (QUOTE HORNERIFY)
(CONS TERM (CONS VARNAME (CDR COEFFICIENTS))))))))))
(PUTPROPS FLEQ MACRO ((X Y)
(NOT (FGREATERP X Y))))
(PUTPROPS FGEQ MACRO ((X Y)
(NOT (FLESSP X Y))))
)
(DECLARE: EVAL@COMPILE
(RPAQQ \EXPONENT.BIAS 127)
(RPAQQ 2PI 6.283185)
(RPAQQ PI 3.141593)
(RPAQQ -PI -3.141593)
(RPAQQ -PI/2 -1.570796)
(RPAQQ PI/2 1.570796)
(RPAQQ 4/PI 1.273239)
(RPAQQ 3PI/2 4.712389)
(RPAQQ PI/4 .7853982)
(RPAQQ -PI/4 -.7853982)
(RPAQQ PI/180 .01745329)
(RPAQQ 180/PI 57.29578)
(RPAQQ -PI/2 -1.570796)
(RPAQQ LN2 .6931472)
(RPAQQ 2↑-126 1.175494E-38)
(CONSTANTS (\EXPONENT.BIAS 127)
(2PI 6.283185)
(PI 3.141593)
(-PI -3.141593)
(-PI/2 -1.570796)
(PI/2 1.570796)
(4/PI 1.273239)
(3PI/2 4.712389)
(PI/4 .7853982)
(-PI/4 -.7853982)
(PI/180 .01745329)
(180/PI 57.29578)
(-PI/2 -1.570796)
(LN2 .6931472)
(2↑-126 1.175494E-38))
)
)
(PUTPROPS AARITH COPYRIGHT ("Xerox Corporation" 1981 1983 1984))
(DECLARE: DONTCOPY
(FILEMAP (NIL (799 12794 (LOG 809 . 2307) (ANTILOG 2309 . 3918) (SIN 3920 . 5337) (ARCSIN 5339 . 6651)
(COS 6653 . 6891) (ARCCOS 6893 . 7177) (TAN 7179 . 8416) (ARCTAN 8418 . 9308) (ARCTAN2 9310 . 10114)
(ATAN 10116 . 12392) (FEXPT 12394 . 12792)))))
STOP