(FILECREATED "14-Aug-85 22:27:27" {ERIS}<LISPCORE>SOURCES>AARITH.;12 20745 changes to: (VARS AARITHCOMS) previous date: "29-Jul-85 15:14:29" {ERIS}<JANDERSSON>CMLFFT>AARITH.;3) (* Copyright (c) 1981, 1983, 1984, 1985 by Xerox Corporation. All rights reserved.) (PRETTYCOMPRINT AARITHCOMS) (RPAQQ AARITHCOMS ((FNS LOG ANTILOG SIN ARCSIN COS ARCCOS TAN ARCTAN ARCTAN2 ATAN FEXPT) (VARS \ANTILOGARRAY \ANTILOGCARRAY \ARCTANARRAY \LOGARRAY \SINARRAY1 \SINARRAY2 \TANARRAY \ATANARRAY) (DECLARE: EVAL@COMPILE DONTCOPY (MACROS (* now obsolete - use POLYEVAL instead) HORNERIFY FLEQ FGEQ) (FILES (LOADCOMP) LLFLOAT) (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) (* hdj "11-Feb-85 17:14") (DECLARE (GLOBALVARS \LOGARRAY)) (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 (POLYEVAL SX \LOGARRAY 8)) (* * 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) (* JAS "19-Jul-85 11:55") (DECLARE (GLOBALVARS \ANTILOGARRAY \ANTILOGCARRAY)) (PROG ((XX (FLOAT X)) FRAC IP SSUM YY) (DECLARE (TYPE FLOATING XX FRAC SSUM YY)) (SETQ YY (FABS XX)) (COND ((GREATERP YY 88.7) (COND ((LESSP XX 0) (RETURN 0.0)) (T (ERROR "FLOATING OVERFLOW" X))))) (SETQ FRAC (FDIFFERENCE YY (FTIMES (CONSTANT (LOG 2)) (SETQ IP (FIX (FTIMES YY (CONSTANT (FQUOTIENT 1.0 (LOG 2))))))))) (SETQ SSUM (POLYEVAL FRAC \ANTILOGARRAY 7)) (* * Polynomial from Handbook of Mathematical Functions (edited by Aramowitz) page 71 accuracy 32 bits of the 24 available!) (* SSUM is in the range .5 to 2 (and series produced .5 to 1)) (SETQ SSUM (QUOTIENT SSUM (ELT \ANTILOGCARRAY (IPLUS IP 127)))) (RETURN (if (FGREATERP XX 0.0) then (FQUOTIENT 1.0 SSUM) else SSUM))))) (SIN (LAMBDA (X RADIANSFLG) (* lmm " 8-Mar-85 07:29") (DECLARE (GLOBALVARS \SINARRAY1 \SINARRAY2)) (PROG ((XX X) X2) (DECLARE (TYPE FLOATP XX 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!)) (SETQ X2 (FTIMES XX (POLYEVAL X2 \SINARRAY1 5))) 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) (SETQ X2 (FTIMES XX (POLYEVAL X2 \SINARRAY2 3)))))))) (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) (* lmm " 8-Mar-85 07:32") (PROGN (DECLARE (GLOBALVARS \SINARRAY1 \SINARRAY2)) (PROG ((XX X) X2) (DECLARE (TYPE FLOATP XX 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))) (SETQ XX (FDIFFERENCE PI/2 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!)) (SETQ X2 (FTIMES XX (POLYEVAL X2 \SINARRAY1 5))) 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) (SETQ X2 (FTIMES XX (POLYEVAL X2 \SINARRAY2 3))))))))) (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) (* lmm " 8-Mar-85 08:29") (DECLARE (GLOBALVARS \TANARRAY)) (PROG ((XX X) FLIPPED Y X2) (DECLARE (TYPE FLOATP XX Y X2)) (SETQ XX (if RADIANSFLG then (FREMAINDER XX PI) else (FTIMES (FREMAINDER XX 180.0) PI/180))) (DECLARE (TYPE FLOATP XX Y X2)) (* 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 (POLYEVAL X2 \TANARRAY 6))) (* POLYNOMIAL APPROXIMATION FROM HANDBOOK OF MATHEMATICAL FUNCTIONS (EDITED BY ABRAMOWITZ) PAGE 76 GOOD TO ALMOST 26 BITS) (RETURN (if FLIPPED then (SETQ Y (FQUOTIENT 1.0 Y)) else Y))))) (ARCTAN (LAMBDA (X RADIANSFLG) (* hdj "11-Feb-85 17:24") (DECLARE (GLOBALVARS \ARCTANARRAY)) (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 (POLYEVAL X2 \ARCTANARRAY 8))) (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) (* hdj "11-Feb-85 17:26") (* 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))) (DECLARE (GLOBALVARS \ATANARRAY)) (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. (POLYEVAL R \ATANARRAY 7))) (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) (* JAS "29-Jul-85 15:13") (* In addition to coercing the args to floating-point, this handles the case of negative values for N) (COND ((EQP A 0.0) 0.0) (T (ANTILOG (FTIMES (OR (FLOATP N) (FLOAT N)) (LOG (OR (FLOATP A) (FLOAT A))))))))) ) (RPAQ \ANTILOGARRAY (READARRAY 8 (QUOTE FLOATP) 0)) (-.0001413161 .001329882 -.00830136 .04165735 -.1666653 .4999999 -1.0 1.0 NIL ) (RPAQ \ANTILOGCARRAY (READARRAY 255 (QUOTE FLOATP) 0)) (5.877474E-39 1.175494E-38 2.350992E-38 4.70198E-38 9.40396E-38 1.880791E-37 3.761582E-37 7.523175E-37 1.504633E-36 3.009266E-36 6.018532E-36 1.203706E-35 2.407413E-35 4.814832E-35 9.629656E-35 1.92593E-34 3.85186E-34 7.70372E-34 1.540746E-33 3.081488E-33 6.162979E-33 1.232595E-32 2.465191E-32 4.930382E-32 9.860764E-32 1.972152E-31 3.944305E-31 7.888613E-31 1.577722E-30 3.155448E-30 6.310891E-30 1.262178E-29 2.524355E-29 5.04871E-29 1.009743E-28 2.019484E-28 4.038968E-28 8.077936E-28 1.615587E-27 3.231174E-27 6.462349E-27 1.29247E-26 2.58494E-26 5.169879E-26 1.033976E-25 2.067952E-25 4.135903E-25 8.271806E-25 1.654361E-24 3.308724E-24 6.617445E-24 1.323489E-23 2.646978E-23 5.293956E-23 1.058791E-22 2.117583E-22 4.235165E-22 8.47033E-22 1.694066E-21 3.388132E-21 6.776264E-21 1.355253E-20 2.710506E-20 5.421011E-20 1.084202E-19 2.168405E-19 4.336809E-19 8.67363E-19 1.734724E-18 3.469447E-18 6.938904E-18 1.387779E-17 2.775558E-17 5.551115E-17 1.110223E-16 2.220446E-16 4.440892E-16 8.881784E-16 1.776359E-15 3.552714E-15 7.105429E-15 1.421086E-14 2.842171E-14 5.684342E-14 1.136868E-13 2.273737E-13 4.547474E-13 9.094948E-13 1.818989E-12 3.637979E-12 7.275959E-12 1.455192E-11 2.910383E-11 5.820766E-11 1.164153E-10 2.328307E-10 4.656613E-10 9.31324E-10 1.862645E-9 3.72529E-9 7.450583E-9 1.490116E-8 2.980232E-8 5.960465E-8 1.192093E-7 2.384186E-7 4.768372E-7 9.536744E-7 1.907349E-6 3.814697E-6 7.629397E-6 .00001525879 .00003051759 .00006103516 .0001220703 .0002441406 .0004882813 .0009765626 .001953125 .00390625 .0078125 .015625 .03125 .0625 .125 .25 .5 1.0 2.0 4.0 8.0 16.0 32.0 64.0 128.0 256.0 512.0 1024.0 2048.0 4096.0 8192.0 16384.0 32768.0 65536.0 131072.0 262144.0 524288.0 1048576.0 2097152.0 4194304.0 8388608.0 16777220.0 33554430.0 67108870.0 1.342177E8 2.684355E8 5.368709E8 1.073742E9 2.147484E9 4.294968E9 8.589936E9 1.717987E10 3.435974E10 6.871948E10 1.37439E11 2.748779E11 5.497558E11 1.099512E12 2.199023E12 4.398047E12 8.796094E12 1.759219E13 3.518437E13 7.036876E13 1.407375E14 2.81475E14 5.6295E14 1.1259E15 2.2518E15 4.5036E15 9.0072E15 1.80144E16 3.60288E16 7.205761E16 1.441152E17 2.882304E17 5.764609E17 1.152922E18 2.305843E18 4.611686E18 9.223372E18 1.844675E19 3.689349E19 7.378698E19 1.47574E20 2.951479E20 5.902958E20 1.180592E21 2.361183E21 4.722367E21 9.444734E21 1.888947E22 3.777893E22 7.555786E22 1.511157E23 3.022315E23 6.044629E23 1.208926E24 2.417852E24 4.835703E24 9.671406E24 1.934281E25 3.868563E25 7.737125E25 1.547425E26 3.09485E26 6.1897E26 1.23794E27 2.47588E27 4.95176E27 9.90352E27 1.980704E28 3.961408E28 7.922816E28 1.584563E29 3.169126E29 6.338253E29 1.267651E30 2.535301E30 5.070602E30 1.01412E31 2.028241E31 4.056482E31 8.112964E31 1.622593E32 3.245186E32 6.490371E32 1.298074E33 2.596148E33 5.192297E33 1.038459E34 2.076919E34 4.153837E34 8.307675E34 1.661535E35 3.32307E35 6.64614E35 1.329228E36 2.658456E36 5.316912E36 1.063382E37 2.126765E37 4.25353E37 8.50706E37 1.701412E38 NIL ) (RPAQ \ARCTANARRAY (READARRAY 9 (QUOTE FLOATP) 0)) (.002866226 -.01616574 .04290961 -.07528964 .1065626 -.142089 .1999355 -.3333315 1.0 NIL ) (RPAQ \LOGARRAY (READARRAY 9 (QUOTE FLOATP) 0)) (-.006453544 .03608849 -.0953294 .1676541 -.2407338 .331799 -.4998741 .9999964 0.0 NIL ) (RPAQ \SINARRAY1 (READARRAY 6 (QUOTE FLOATP) 0)) (-2.39E-8 2.7526E-6 -.000198409 .008333332 -.1666667 1.0 NIL ) (RPAQ \SINARRAY2 (READARRAY 4 (QUOTE FLOATP) 0)) (-.0000359544 .002490007 -.08074545 .7853982 NIL ) (RPAQ \TANARRAY (READARRAY 7 (QUOTE FLOATP) 0)) (.00951681 .002900525 .02456509 .05337406 .1333924 .3333314 1.0 NIL ) (RPAQ \ATANARRAY (READARRAY 8 (QUOTE FLOATP) 0)) (-.004054058 .02186123 -.05590989 .09642004 -.1390853 .1994654 -.3332986 .9999994 NIL ) (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)))) ) (FILESLOAD (LOADCOMP) LLFLOAT) (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 1985)) (DECLARE: DONTCOPY (FILEMAP (NIL (1032 14964 (LOG 1042 . 2616) (ANTILOG 2618 . 3859) (SIN 3861 . 5413) (ARCSIN 5415 . 6899) (COS 6901 . 8554) (ARCCOS 8556 . 8868) (TAN 8870 . 10255) (ARCTAN 10257 . 11187) (ARCTAN2 11189 . 12073) (ATAN 12075 . 14484) (FEXPT 14486 . 14962))))) STOP