(FILECREATED " 4-Jul-86 01:28:23" {ERIS}<LISPCORE>LIBRARY>CMLFLOAT.;34      changes to:  (VARS CMLFLOATCOMS)      previous date: " 2-Jul-86 10:29:42" {ERIS}<LISPCORE>LIBRARY>CMLFLOAT.;33)(* Copyright (c) 1986 by Xerox Corporation. All rights reserved.)(PRETTYCOMPRINT CMLFLOATCOMS)(RPAQQ CMLFLOATCOMS        [(* * CMLFLOAT -- Covering sections 12.5-12.5.3 irrational, transcendental, exponential,            logarithmic, trigonometric, and hyperbolic functions.)        (COMS (FNS %%MAKE-ARRAY %%FLOAT))        (DECLARE: EVAL@COMPILE DONTCOPY               (VARIABLES %%SINGLE-FLOAT-EXPONENT-LENGTH %%SINGLE-FLOAT-MANTISSA-LENGTH %%FLOAT-ONE                       %%FLOAT-HALF %%SQRT-HALF %%E SQRT3 2-SQRT3 INV-2-SQRT3 %%ASIN-EPS %%ASIN-P1                       %%ASIN-P2 %%ASIN-Q0 %%ASIN-Q1 %%ASIN-Q2 %%FLOAT-HALF %%ZERO)               (FILES (LOADCOMP)                      LLFLOAT))        (VARIABLES PI)        (COMS (* Random old spice junk. *)              (FUNCTIONS SCALE-FLOAT-BY POLY-EVAL))        (COMS (* EXP *)              [COMS (CONSTANTS (%%LOG-BASE2-E (%%FLOAT 16312 43579)))                    (* * %%EXP-POLY contains P and Q coefficients of Harris et al EXPB 1103 rational                        approximation to (EXPT 2 X)                       in interval (0 .125)                       %. %%EXP-TABLE contains values of powers (EXPT 2 (/ N 8)) . *)                    (INITVARS [%%EXP-POLY (%%MAKE-ARRAY (LIST (%%FLOAT 15549 17659)                                                              (%%FLOAT 16256 0)                                                              (%%FLOAT 16801 38273)                                                              (%%FLOAT 17257 7717)                                                              (%%FLOAT 17597 11739)                                                              (%%FLOAT 17800 30401]                           (%%EXP-TABLE (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                                            (%%FLOAT 16267 38338)                                                            (%%FLOAT 16280 14320)                                                            (%%FLOAT 16293 65239)                                                            (%%FLOAT 16309 1267)                                                            (%%FLOAT 16325 26410)                                                            (%%FLOAT 16343 17661)                                                            (%%FLOAT 16362 49351]              (FNS EXP %%EXP-FLOAT %%EXP-COMPLEX))        (COMS (* EXPT *)              (FNS CL:EXPT %%EXPT-INTEGER %%EXPT-FLOAT %%EXPT-COMPLEX %%EXPT-COMPLEX-POWER))        (COMS (* LOG *)              [COMS [DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%LOG2 (%%FLOAT 16177 29208))                                                           (%%SQRT2 (%%FLOAT 16309 1267]                    (* * %%LOG-PPOLY and %%LOG-QPOLY contain P and Q coefficients of Harris et al                        LOGE 2707 rational approximation to (LOG X)                       in interval ((SQRT .5)                                    (SQRT 2)) . *)                    (INITVARS [%%LOG-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16042 22803)                                                               (%%FLOAT 49484 23590)                                                               (%%FLOAT 17044 17982)                                                               (%%FLOAT 49926 37153)                                                               (%%FLOAT 17046 5367]                           (%%LOG-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                                            (%%FLOAT 49512 9103)                                                            (%%FLOAT 16992 42274)                                                            (%%FLOAT 49823 38048)                                                            (%%FLOAT 16918 5367]              (FNS CL:LOG %%LOG-FLOAT %%LOG-COMPLEX))        (FNS CL:SQRT)        (DECLARE: EVAL@COMPILE DONTCOPY (VARIABLES %%2PI %%PI %%PI/2 %%-PI/2 %%PI/4 %%-PI/4 SIXTH-PI                                                %%2/PI))        (COMS (* SIN *)              (COMS (* * %%SIN-EPSILON is sufficiently small that (SIN X)                       = X for X in interval (0 %%SIN-EPSILON)                       %. It suffices to take %%SIN-EPSILON a little bit smaller than                       (SQRT (CL:* 6 SINGLE-FLOAT-EPSILON))                       which we get by the Taylor series expansion (SIN X)                       =                       (+ X (/ (EXPT X 3)                               6)                          ...)                       (The relative error caused by ommitting (/ (EXPT X 3)                                                                  6)                            isn't observable.)                       Comparison against %%SIN-EPSILON is used to avoid POLYEVAL microcode underflow                        when computing SIN. *)                    (DECLARE: EVAL@COMPILE DONTCOPY (VARIABLES %%SIN-EPSILON))                    (* * %%SIN-PPOLY and %%SIN-QPOLY contain P and Q coefficients of Harris et al SIN                        3374 rational approximation to (SIN X)                       in interval (0 (/ PI 2)) . *)                    (VARIABLES %%SIN-PPOLY %%SIN-QPOLY))              (FNS CL:SIN %%SIN-FLOAT %%SIN-COMPLEX))        (COMS (* COS *)              (FNS CL:COS %%COS-COMPLEX))        (COMS [COMS (* * %%TAN-EPSILON is sufficiently small that (TAN X)                       = X for X in interval (0 %%TAN-EPSILON)                       %. It suffices to take %%TAN-EPSILON a little bit smaller than                       (SQRT (CL:* 3 SINGLE-FLOAT-EPSILON))                       which we get by the Taylor series expansion (TAN X)                       =                       (+ X (/ (EXPT X 3)                               3)                          ...)                       (The relative error caused by ommitting (/ (EXPT X 3)                                                                  3)                            isn't observable.)                       Comparison against %%TAN-EPSILON is used to avoid POLYEVAL microcode underflow                        when computing TAN. *)                    (DECLARE: EVAL@COMPILE DONTCOPY (VARIABLES %%TAN-EPSILON))                    (* * %%TAN-PPOLY and %%TAN-QPOLY contain P and Q coefficients of Harris et al TAN                        4288 rational approximation to (TAN X)                       in interval (-PI/4 PI/4) . *)                    (INITVARS [%%TAN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 49803 44451)                                                               (%%FLOAT 18254 23362)                                                               (%%FLOAT 51983 57389)                                                               (%%FLOAT 19947 57965)                                                               (%%FLOAT 53162 16965]                           (%%TAN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                                            (%%FLOAT 50453 51515)                                                            (%%FLOAT 18760 65370)                                                            (%%FLOAT 52376 39517)                                                            (%%FLOAT 20221 24840)                                                            (%%FLOAT 53208 51139]              (FNS CL:TAN %%TAN-FLOAT %%TAN-COMPLEX))        (COMS [COMS (* * %%ASIN-EPSILON is sufficiently small that (ASIN X)                       = X for X in interval (0 %%ASIN-EPSILON)                       %. It suffices to take %%ASIN-EPSILON a little bit smaller than                       (CL:* 2 SINGLE-FLOAT-EPSILON)                       which we get by the Taylor series expansion (ASIN X)                       =                       (+ X (/ (EXPT X 3)                               6)                          ...)                       (The relative error caused by ommitting (/ (EXPT X 3)                                                                  6)                            isn't observable.)                       Comparison against %%ASIN-EPSILON is used to avoid POLYEVAL microcode                        underflow when computing SIN. *)                    (DECLARE: EVAL@COMPILE DONTCOPY (VARIABLES %%ASIN-EPSILON))                    (* * %%ASIN-PPOLY and %%ASIN-QPOLY contain P and Q coefficients of Harris et al                        ARCSN 4671 rational approximation to (ASIN X)                       in interval (0 (SQRT .5)) . *)                    (INITVARS [%%ASIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16007 50045)                                                                (%%FLOAT 49549 8020)                                                                (%%FLOAT 17236 15848)                                                                (%%FLOAT 50285 63464)                                                                (%%FLOAT 17650 31235)                                                                (%%FLOAT 50403 62852)                                                                (%%FLOAT 17440 39471]                           (%%ASIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                                             (%%FLOAT 49672 25817)                                                             (%%FLOAT 17308 55260)                                                             (%%FLOAT 50326 38098)                                                             (%%FLOAT 17674 22210)                                                             (%%FLOAT 50417 22451)                                                             (%%FLOAT 17440 39471]              (FNS ASIN %%ASIN-FLOAT %%ASIN-COMPLEX))        (COMS (* ACOS *)              (FNS ACOS %%ACOS-COMPLEX))        (COMS (* ATAN *)              (FNS CL:ATAN %%ATAN1 %%ATAN2 %%ATAN-DOMAIN-CHECK %%ATAN-FLOAT %%ATAN-COMPLEX))        (FNS CIS)        (FNS SINH COSH TANH)        (FNS ASINH ACOSH ATANH %%ATANH-DOMAIN-CHECK)        (PROP FILETYPE CMLFLOAT)        (DECLARE: DONTEVAL@LOAD DOEVAL@COMPILE DONTCOPY COMPILERVARS               (ADDVARS (NLAMA)                      (NLAML)                      (LAMA %%ATANH-DOMAIN-CHECK ATANH ACOSH ASINH TANH COSH SINH CIS %%ATAN-COMPLEX                             %%ATAN-FLOAT %%ATAN-DOMAIN-CHECK %%ATAN2 %%ATAN1 CL:ATAN %%ACOS-COMPLEX                             ACOS %%ASIN-COMPLEX ASIN %%TAN-COMPLEX CL:TAN %%COS-COMPLEX CL:COS                             %%SIN-COMPLEX CL:SIN CL:SQRT %%LOG-COMPLEX CL:LOG %%EXPT-COMPLEX-POWER                             %%EXPT-COMPLEX %%EXPT-INTEGER CL:EXPT %%EXP-COMPLEX %%EXP-FLOAT EXP])(* * CMLFLOAT -- Covering sections 12.5-12.5.3 irrational, transcendental, exponential, logarithmic, trigonometric, and hyperbolic functions.)(DEFINEQ(%%MAKE-ARRAY  [LAMBDA (LIST)                                             (* kbr: "15-May-86 18:09")                                                             (* Function to build Interlisp arrays.                                                             I would prefer to build Common Lisp                                                              arrays, but I don't think POLYEVAL                                                              opcode will understand them.                                                             *)    (PROG (ARRAY)          (SETQ ARRAY (ARRAY (LENGTH LIST)                             (QUOTE FLOATP)                             0.0 0))          (for ELEMENT in LIST as I from 0 do (SETA ARRAY I ELEMENT))          (RETURN ARRAY])(%%FLOAT  [LAMBDA (HIWORD LOWORD)                                    (* kbr: "14-May-86 16:43")    (\FLOATBOX (\VAG2 HIWORD LOWORD]))(DECLARE: EVAL@COMPILE DONTCOPY (DEFCONSTANT %%SINGLE-FLOAT-EXPONENT-LENGTH 8)(DEFCONSTANT %%SINGLE-FLOAT-MANTISSA-LENGTH 23)(DEFCONSTANT %%FLOAT-ONE 1.0)(DEFCONSTANT %%FLOAT-HALF .5)(DEFCONSTANT %%SQRT-HALF (%%FLOAT 16181 1267) )(DEFCONSTANT %%E (%%FLOAT 16429 63530) )(DEFCONSTANT SQRT3 1.730525)(DEFCONSTANT 2-SQRT3 .2679482)(DEFCONSTANT INV-2-SQRT3 3.732071)(DEFCONSTANT %%ASIN-EPS (**FLOAT** 1 -10) )(DEFCONSTANT %%ASIN-P1 .08333311)(DEFCONSTANT %%ASIN-P2 -.04500659)(DEFCONSTANT %%ASIN-Q0 .5)(DEFCONSTANT %%ASIN-Q1 -.495078)(DEFCONSTANT %%ASIN-Q2 .08922788)(DEFCONSTANT %%FLOAT-HALF .5)(DEFCONSTANT %%ZERO 0.0)(FILESLOAD (LOADCOMP)       LLFLOAT))(DEFCONSTANT PI (%%FLOAT 16457 4059) )(* Random old spice junk. *)(DEFMACRO SCALE-FLOAT-BY (E F) (BQUOTE (SCALE-FLOAT (\, F)                                              (\, E))))(DEFMACRO POLY-EVAL (X &REST LIST) (CL:DO ((LIST (CDR (CL:REVERSE LIST))                                                 (CDR LIST))                                           (FORM (CAR (LAST LIST))                                                 (BQUOTE (+ (\, (CAR LIST))                                                            (CL:* (\, X)                                                                  (\, FORM))))))                                          ((NULL LIST)                                           FORM)))(* EXP *)(DECLARE: EVAL@COMPILE (RPAQ %%LOG-BASE2-E (%%FLOAT 16312 43579))(CONSTANTS (%%LOG-BASE2-E (%%FLOAT 16312 43579))))(* * %%EXP-POLY contains P and Q coefficients of Harris et al EXPB 1103 rational approximation to (EXPT 2 X) in interval (0 .125) %. %%EXP-TABLE contains values of powers (EXPT 2 (/ N 8)) . *)(RPAQ? %%EXP-POLY (%%MAKE-ARRAY (LIST (%%FLOAT 15549 17659)                                      (%%FLOAT 16256 0)                                      (%%FLOAT 16801 38273)                                      (%%FLOAT 17257 7717)                                      (%%FLOAT 17597 11739)                                      (%%FLOAT 17800 30401))))(RPAQ? %%EXP-TABLE (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                       (%%FLOAT 16267 38338)                                       (%%FLOAT 16280 14320)                                       (%%FLOAT 16293 65239)                                       (%%FLOAT 16309 1267)                                       (%%FLOAT 16325 26410)                                       (%%FLOAT 16343 17661)                                       (%%FLOAT 16362 49351))))(DEFINEQ(EXP  [CL:LAMBDA (NUMBER)                                        (* lmm "27-Jun-86 22:40")                                                             (* Calculates e to the given power,                                                              where e is the base of natural                                                              logarithms. *)         (CTYPECASE NUMBER (COMPLEX (%%EXP-COMPLEX NUMBER))                (NUMBER (%%EXP-FLOAT (FLOAT NUMBER])(%%EXP-FLOAT  (CL:LAMBDA (X)                                             (* kbr: "30-May-86 15:50")                    (* * (EXP X) for float X calculated via EXPB 1103 rational approximation of           Harris et al. *)         (PROG (M N R ANSWER)               (DECLARE (GLOBALVARS %%EXP-TABLE %%EXP-POLY))                    (* * First, the problem of (EXP X) is converted into a problem          (EXPT 2 Y) where Y = (CL:* %%LOG-BASE2-E X)%.          Then range reduction is accomplished via          (EXPT 2 Y) = (CL:* (EXPT 2 M) (EXPT 2 (/ N 8))          (EXPT 2 R)) where M and N are integers and R is a float in the interval          (0.0 .125) After M N R are determined, (EXPT 2 M) is effected by scaling,          (EXPT 2 (/ N 8)) is found by table lookup, and          (EXPT 2 R) is calculated by rational approximation EXPB 1103 of Harris et al.          *)               (MULTIPLE-VALUE-SETQ (M R)                      (CL:FLOOR (CL:* %%LOG-BASE2-E X)))               (MULTIPLE-VALUE-SETQ (N R)                      (CL:FLOOR R .125))               (SETQ ANSWER (SCALE-FLOAT (CL:* (ELT %%EXP-TABLE N)                                               (/ (POLYEVAL R %%EXP-POLY 5)                                                  (POLYEVAL (- R)                                                         %%EXP-POLY 5)))                                   M))               (RETURN ANSWER))))(%%EXP-COMPLEX  [CL:LAMBDA (Z)                                             (* compute exp (z) where z is complex                                                              is called by exp *)         (LET* ((X (REALPART Z))                (Y (IMAGPART Z)))               (CL:* (EXP X)                     (CIS Y]))(* EXPT *)(DEFINEQ(CL:EXPT  [CL:LAMBDA (X N)                                           (* lmm "27-Jun-86 22:52")                    (* This function calculates X raised to the nth power.          It separates the cases by the type of N for efficiency reasons, as powers can           be calculated more efficiently if N is a positive integer, Therefore, All           integer values of N are calculated as positive integers, and inverted if           negative. *)         (COND            ((AND (RATIONALP X)                  (INTEGERP N))             (%%EXPT-INTEGER X N))            ((= X 0)             (CL:IF (PLUSP N)                    X                    (CL:ERROR "~A to a non-positive power ~A." X N)))            ((AND (COMPLEXP X)                  (INTEGERP N))             (%%EXPT-INTEGER X N))            ((COMPLEXP N)             (CL:IF (OR (COMPLEXP X)                        (PLUSP X))                    (%%EXPT-COMPLEX-POWER X N)                    (CL:ERROR "~A negative number to a complex power ~A." X N)))            ((COMPLEXP X)             (%%EXPT-COMPLEX X N))            ((AND (NOT (INTEGERP N))                  (MINUSP X))             (CL:ERROR "Negative number ~A to non-integral power ~A." X N))            (T (%%EXPT-FLOAT (FLOAT X)                      (FLOAT N])(%%EXPT-INTEGER  [CL:LAMBDA (BASE POWER)                                    (* kbr: "28-May-86 15:58")                    (* * (EXPT BASE POWER) where POWER is an integer.          *)         (COND            [(MINUSP POWER)             (/ (%%EXPT-INTEGER BASE (- POWER]            ((AND (RATIONALP BASE)                  (NOT (INTEGERP BASE)))             (%%MAKE-RATIO (%%EXPT-INTEGER (NUMERATOR BASE)                                  POWER)                    (%%EXPT-INTEGER (DENOMINATOR BASE)                           POWER)))            ((AND (INTEGERP BASE)                  (= BASE 2))             (ASH 1 POWER))            (T (CL:DO ((NEXTN (ASH POWER -1)                              (ASH POWER -1))                       (TOTAL (CL:IF (ODDP POWER)                                     BASE 1)                              (CL:IF (ODDP POWER)                                     (CL:* BASE TOTAL)                                     TOTAL)))                      ((ZEROP NEXTN)                       TOTAL)                      (SETQ BASE (CL:* BASE BASE))                      (SETQ POWER NEXTN])(%%EXPT-FLOAT  [LAMBDA (X Y)                                              (* kbr: "29-May-86 00:03")                    (* * (EXPT X Y) where X is a nonnegative float and Y is a float.          *)    (COND       ((= X 0.0)        0.0)       (T (EXP (CL:* Y (CL:LOG X])(%%EXPT-COMPLEX  [CL:LAMBDA (Z N)                                           (* compute (complex) ^n where n is an                                                              integer some round-off error if n is                                                              not a fixnum *)         (CL:* (CL:EXPT (%%COMPLEX-ABS Z)                      N)               (CIS (CL:* N (PHASE Z])(%%EXPT-COMPLEX-POWER  [CL:LAMBDA (Z W)                                           (* this function computes z ^ w where                                                              w is a complex number it can also be                                                              used for any positive real number.                                                             *)         (EXP (CL:* W (CL:LOG Z]))(* LOG *)(DECLARE: EVAL@COMPILE DONTCOPY (DECLARE: EVAL@COMPILE (RPAQ %%LOG2 (%%FLOAT 16177 29208))(RPAQ %%SQRT2 (%%FLOAT 16309 1267))(CONSTANTS (%%LOG2 (%%FLOAT 16177 29208))       (%%SQRT2 (%%FLOAT 16309 1267)))))(* * %%LOG-PPOLY and %%LOG-QPOLY contain P and Q coefficients of Harris et al LOGE 2707 rational approximation to (LOG X) in interval ((SQRT .5) (SQRT 2)) . *)(RPAQ? %%LOG-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16042 22803)                                       (%%FLOAT 49484 23590)                                       (%%FLOAT 17044 17982)                                       (%%FLOAT 49926 37153)                                       (%%FLOAT 17046 5367))))(RPAQ? %%LOG-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                       (%%FLOAT 49512 9103)                                       (%%FLOAT 16992 42274)                                       (%%FLOAT 49823 38048)                                       (%%FLOAT 16918 5367))))(DEFINEQ(CL:LOG  [CL:LAMBDA (NUMBER &OPTIONAL (BASE NIL BASE-SUPPLIED))     (* lmm "27-Jun-86 22:40")         (COND            (BASE-SUPPLIED (/ (CL:LOG NUMBER)                              (CL:LOG BASE)))            (T (CTYPECASE NUMBER (COMPLEX (%%LOG-COMPLEX NUMBER))                      (NUMBER (CL:IF (MINUSP NUMBER)                                     (COMPLEX (%%LOG-FLOAT (FLOAT (- NUMBER)))                                            %%PI)                                     (%%LOG-FLOAT (FLOAT NUMBER])(%%LOG-FLOAT  [LAMBDA (X)                                                (* kbr: "30-May-86 00:04")                    (* * (EXP X) for float X calculated via EXPB 1103 rational approximation of           Harris et al. *)                    (* * (LOG X) for nonnegative float X. *)    (PROG (R EXP Z Z2 ANSWER)          (DECLARE (GLOBALVARS %%LOG-PPOLY %%LOG-QPOLY))                    (* * NOTE: Probably best not to declare type of variables inside this routine           for now. I've found that FLOATP smashing combined with FLOATP declarations           compiles wrong. *)          (COND             ((NOT (FGREATERP X 0.0))              (ERROR "LOG OF ZERO:" X)))                    (* * Range reduce to an R in interval ((SQRT .5)          (SQRT 2)) via identity (LOG X) = (+ (LOG R)          (CL:* %%LOG-2 EXP)) for a suitable integer EXP.          This reduction depends on and smashes the FLOATP representation of R.          *)          (SETQ R (\FLOAT.BOX X))          (SETQ EXP (- (fetch (FLOATP EXPONENT) of R)                       \EXPONENT.BIAS))          (replace (FLOATP EXPONENT) of R with \EXPONENT.BIAS)          [COND             ((> R %%SQRT2)              (SETQ EXP (1+ EXP))              (SETQ R (/ R 2.0]                    (* * (LOG R) is calculated by rational approximation LOGE 2707 of Harris et al.          *)          (SETQ Z (/ (1- R)                     (1+ R)))          (SETQ Z2 (CL:* Z Z))          (SETQ ANSWER (+ (CL:* Z (/ (POLYEVAL Z2 %%LOG-PPOLY 4)                                     (POLYEVAL Z2 %%LOG-QPOLY 4)))                          (CL:* %%LOG2 EXP)))          (RETURN ANSWER])(%%LOG-COMPLEX  (CL:LAMBDA (Z)                                             (* natural log of a complex number *)         (COMPLEX (CL:LOG (%%COMPLEX-ABS Z))                (PHASE Z)))))(DEFINEQ(CL:SQRT  [CL:LAMBDA (X)                                             (* lmm " 1-Jul-86 15:56")                                                             (* currently defined in terms of                                                              Interlisp's SQRT)         (COND            ((COMPLEXP X)             (PROG (ABS PHASE A B C D E ANSWER)                   (DECLARE (TYPE FLOATP ABS PHASE A B C D E))                   (SETQ A (FLOAT (COMPLEX-REALPART X)))                   (SETQ B (FLOAT (COMPLEX-IMAGPART X)))                    (* * Make initial guess. *)                   [SETQ ABS (SQRT (SQRT (+ (CL:* A A)                                            (CL:* B B]                   (SETQ PHASE (/ (PHASE X)                                  2.0))                   (SETQ C (CL:* ABS (CL:COS PHASE)))                   (SETQ D (CL:* ABS (CL:SIN PHASE)))                    (* * Newton's method. *)                   [for I from 1 to 4 do (SETQ E (+ (CL:* C C)                                                    (CL:* D D)))                                         [SETQ C (CL:* .5 (+ C (/ (+ (CL:* A C)                                                                     (CL:* B D))                                                                  E]                                         (SETQ D (CL:* .5 (+ D (/ (- (CL:* B C)                                                                     (CL:* A D))                                                                  E]                   (SETQ ANSWER (COMPLEX C D))                   (RETURN ANSWER)))            [(MINUSP X)             (COMPLEX 0 (SQRT (- X]            (T (SQRT (FLOAT X]))(DECLARE: EVAL@COMPILE DONTCOPY (DEFCONSTANT %%2PI (%%FLOAT 16585 4059) )(DEFCONSTANT %%PI (%%FLOAT 16457 4059) )(DEFCONSTANT %%PI/2 (%%FLOAT 16329 4059) )(DEFCONSTANT %%-PI/2 (%%FLOAT 49097 4059) )(DEFCONSTANT %%PI/4 (%%FLOAT 16201 4059) )(DEFCONSTANT %%-PI/4 (%%FLOAT 48969 4059) )(DEFCONSTANT SIXTH-PI .5236)(DEFCONSTANT %%2/PI (%%FLOAT 16162 63875) ))(* SIN *)(* * %%SIN-EPSILON is sufficiently small that (SIN X) = X for X in interval (0 %%SIN-EPSILON) %. It suffices to take %%SIN-EPSILON a little bit smaller than (SQRT (CL:* 6 SINGLE-FLOAT-EPSILON)) which we get by the Taylor series expansion (SIN X) = (+ X (/ (EXPT X 3) 6) ...) (The relative error caused by ommitting (/ (EXPT X 3) 6) isn't observable.) Comparison against %%SIN-EPSILON is used to avoid POLYEVAL microcode underflow when computing SIN. *)(DECLARE: EVAL@COMPILE DONTCOPY (DEFCONSTANT %%SIN-EPSILON (%%FLOAT 14720 0) ))(* * %%SIN-PPOLY and %%SIN-QPOLY contain P and Q coefficients of Harris et al SIN 3374 rational approximation to (SIN X) in interval (0 (/ PI 2)) . *)(DEFVAR %%SIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 50193 50772)                                        (%%FLOAT 18371 47202)                                        (%%FLOAT 51905 54216)                                        (%%FLOAT 19748 29063)                                        (%%FLOAT 52951 41140)                                        (%%FLOAT 20368 50691))) )(DEFVAR %%SIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                        (%%FLOAT 17252 49821)                                        (%%FLOAT 18142 34466)                                        (%%FLOAT 18957 45135)                                        (%%FLOAT 19685 22067)                                        (%%FLOAT 20280 21713))) )(DEFINEQ(CL:SIN  [CL:LAMBDA (X)                                             (* kbr: " 7-Apr-86 22:57")         (CTYPECASE X (COMPLEX (%%SIN-COMPLEX X))                (FLOAT (%%SIN-FLOAT X NIL))                (NUMBER (%%SIN-FLOAT (FLOAT X)                               NIL])(%%SIN-FLOAT  [LAMBDA (X COS-FLAG)                                       (* kbr: "28-May-86 21:40")                    (* * SIN of a FLOAT X calculated via SIN 3374 rational approximation of Harris           et al. *)    (PROG (R SIGN R2 ANSWER)          (DECLARE (GLOBALVARS %%SIN-PPOLY %%SIN-QPOLY)                 (TYPE FLOATP X R SIGN R2 ANSWER))                    (* * If this function called by COS then use          (COS X) = (SIN (+ X %%PI/2)) *)          (SETQ R (CL:IF COS-FLAG (+ X %%PI/2)                         X))                    (* * First range reduce to interval (0 %%2PI) by          (SIN X) = (SIN (MOD X %%2PI)) . *)          (SETQ R (CL:MOD R %%2PI))                    (* * Next range reduce to interval (0 PI) by          (SIN (+ X PI)) = (minus (SIN X)) *)          (COND             ((> R %%PI)              (SETQ SIGN -1.0)              (SETQ R (- R %%PI)))             (T (SETQ SIGN 1.0)))                    (* * Next range reduce to interval (0 %%PI/2) by          (SIN (+ X %%PI/2)) = (SIN (minus %%PI/2 X)) *)          [COND             ((> R %%PI/2)              (SETQ R (- %%PI R]          (COND             ((< R %%SIN-EPSILON)                    (* * If R is in the interval (0 %%SIN-EPSILON) then          (SIN R) = R to the precision that we can offer.          Return R because (1) it is desirable that          (SIN R) = R exactly for small R and (2) microcode POLYEVAL will underflow on           sufficiently small positive R. *)              (RETURN R)))                    (* * TBW: Should get rid of following (SETQ R          (/ R %%PI/2)) line by adjusting %%SIN-PPOLY %%SIN-QPOLY to compensate.          What we have is correct but slightly suboptimal.          *)          (SETQ R (/ R %%PI/2))                    (* * Now use SIN 3374 rational approximation of Harris et al.          which works on interval (0 %%PI/2) *)          (SETQ R2 (CL:* R R))          [SETQ ANSWER (CL:* SIGN R (/ (POLYEVAL R2 %%SIN-PPOLY 5)                                       (POLYEVAL R2 %%SIN-QPOLY 5]          (RETURN ANSWER])(%%SIN-COMPLEX  [CL:LAMBDA (Z)                                             (* sin of a complex number *)         (LET* ((X (REALPART Z))                (Y (IMAGPART Z)))               (COMPLEX (CL:* (CL:SIN X)                              (COSH Y))                      (CL:* (CL:COS X)                            (SINH Y]))(* COS *)(DEFINEQ(CL:COS  [CL:LAMBDA (X)                                             (* lmm " 1-Jul-86 15:55")         (CTYPECASE X (FLOAT (%%SIN-FLOAT X T))                (COMPLEX (%%COS-COMPLEX X))                (NUMBER (%%SIN-FLOAT (FLOAT X)                               T])(%%COS-COMPLEX  [CL:LAMBDA (Z)                                             (* cosine of a complex number *)         (LET* ((X (REALPART Z))                (Y (IMAGPART Z)))               (COMPLEX (CL:* (CL:COS X)                              (COSH Y))                      (- (CL:* (CL:SIN X)                               (SINH Y]))(* * %%TAN-EPSILON is sufficiently small that (TAN X) = X for X in interval (0 %%TAN-EPSILON) %. It suffices to take %%TAN-EPSILON a little bit smaller than (SQRT (CL:* 3 SINGLE-FLOAT-EPSILON)) which we get by the Taylor series expansion (TAN X) = (+ X (/ (EXPT X 3) 3) ...) (The relative error caused by ommitting (/ (EXPT X 3) 3) isn't observable.) Comparison against %%TAN-EPSILON is used to avoid POLYEVAL microcode underflow when computing TAN. *)(DECLARE: EVAL@COMPILE DONTCOPY (DEFCONSTANT %%TAN-EPSILON (%%FLOAT 14720 0) ))(* * %%TAN-PPOLY and %%TAN-QPOLY contain P and Q coefficients of Harris et al TAN 4288 rational approximation to (TAN X) in interval (-PI/4 PI/4) . *)(RPAQ? %%TAN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 49803 44451)                                       (%%FLOAT 18254 23362)                                       (%%FLOAT 51983 57389)                                       (%%FLOAT 19947 57965)                                       (%%FLOAT 53162 16965))))(RPAQ? %%TAN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                       (%%FLOAT 50453 51515)                                       (%%FLOAT 18760 65370)                                       (%%FLOAT 52376 39517)                                       (%%FLOAT 20221 24840)                                       (%%FLOAT 53208 51139))))(DEFINEQ(CL:TAN  [CL:LAMBDA (X)                                             (* lmm "27-Jun-86 22:41")         (CTYPECASE X (COMPLEX (%%TAN-COMPLEX X))                (NUMBER (%%TAN-FLOAT (FLOAT X])(%%TAN-FLOAT  [LAMBDA (X)                                                (* kbr: "28-May-86 22:36")                    (* * TAN of a FLOAT X calculated via TAN 4288 rational approximation of Harris           et al. *)    (PROG (R R2 RECIPFLG ANSWER)          (DECLARE (GLOBALVARS %%TAN-PPOLY %%TAN-QPOLY)                 (TYPE FLOATP X R R2 ANSWER))                    (* * First, range reduce to (0 PI) *)          (SETQ R (CL:MOD X PI))                    (* * Next, range reduce to (-PI/4 PI/4) using          (TAN (+ X PI/2)) = (TAN (minus PI/2 X)) to get into interval          (-PI/2 PI/2) and then (TAN X) = (/ (TAN (minus PI/2 X))) to get into interval          (-PI/4 PI/4) *)          [COND             [(> R %%PI/2)              (SETQ R (- R %%PI))              (COND                 ((< R %%-PI/4)                  (SETQ RECIPFLG T)                  (SETQ R (- %%-PI/2 R]             (T (COND                   ((> R %%PI/4)                    (SETQ RECIPFLG T)                    (SETQ R (- %%PI/2 R]          (COND             ((< (ABS R)               %%TAN-EPSILON)                    (* * If R is in the interval (0 %%TAN-EPSILON) then          (TAN R) = R to the precision that we can offer.          Return R because (1) it is desirable that          (TAN R) = R exactly for small R and (2) microcode POLYEVAL will underflow on           sufficiently small positive R. *)              (SETQ ANSWER (COND                              (RECIPFLG (/ R))                              (T R)))              (RETURN ANSWER)))                    (* * TBW: Should get rid of following (SETQ R          (/ R %%PI/2)) line by adjusting %%SIN-PPOLY %%SIN-QPOLY to compensate.          What we have is correct but slightly suboptimal.          *)          (SETQ R (/ R %%PI/4))                    (* * Now use TAN 4288 rational approximation of Harris et al.          which works on interval (0 %%PI/4) *)          (SETQ R2 (CL:* R R))          [SETQ ANSWER (CL:* R (/ (POLYEVAL R2 %%TAN-PPOLY 5)                                  (POLYEVAL R2 %%TAN-QPOLY 6]          [COND             (RECIPFLG (SETQ ANSWER (/ ANSWER]          (RETURN ANSWER])(%%TAN-COMPLEX  [CL:LAMBDA (Z)                                             (* tan of a complex number there was a                                                              nicer algorithm but it turned out not                                                              to work so well. *)         (LET* ((NUM (CL:SIN Z))                (DENOM (CL:COS Z)))               (CL:IF (ZEROP DENOM)                      (CL:ERROR "~S undefined tangent." Z)                      (/ NUM DENOM]))(* * %%ASIN-EPSILON is sufficiently small that (ASIN X) = X for X in interval (0 %%ASIN-EPSILON) %. It suffices to take %%ASIN-EPSILON a little bit smaller than (CL:* 2 SINGLE-FLOAT-EPSILON) which we get by the Taylor series expansion (ASIN X) = (+ X (/ (EXPT X 3) 6) ...) (The relative error caused by ommitting (/ (EXPT X 3) 6) isn't observable.) Comparison against %%ASIN-EPSILON is used to avoid POLYEVAL microcode underflow when computing SIN. *)(DECLARE: EVAL@COMPILE DONTCOPY (DEFCONSTANT %%ASIN-EPSILON (%%FLOAT 14720 0) ))(* * %%ASIN-PPOLY and %%ASIN-QPOLY contain P and Q coefficients of Harris et al ARCSN 4671 rational approximation to (ASIN X) in interval (0 (SQRT .5)) . *)(RPAQ? %%ASIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16007 50045)                                        (%%FLOAT 49549 8020)                                        (%%FLOAT 17236 15848)                                        (%%FLOAT 50285 63464)                                        (%%FLOAT 17650 31235)                                        (%%FLOAT 50403 62852)                                        (%%FLOAT 17440 39471))))(RPAQ? %%ASIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)                                        (%%FLOAT 49672 25817)                                        (%%FLOAT 17308 55260)                                        (%%FLOAT 50326 38098)                                        (%%FLOAT 17674 22210)                                        (%%FLOAT 50417 22451)                                        (%%FLOAT 17440 39471))))(DEFINEQ(ASIN  [CL:LAMBDA (X)                                             (* lmm "27-Jun-86 22:20")         (CTYPECASE X (FLOAT (%%ASIN-FLOAT X))                (COMPLEX (%%ASIN-COMPLEX X))                (NUMBER (%%ASIN-FLOAT (FLOAT X])(%%ASIN-FLOAT  [LAMBDA (X ACOS-FLAG)                                      (* kbr: "30-May-86 16:46")                    (* * (ASIN X) for float X calculated via ARCSN 4671 rational approximation of           Harris et al. *)    (PROG (NEGATIVE REDUCED R R2 ANSWER)          (DECLARE (GLOBALVARS %%ASIN-PPOLY %%ASIN-QPOLY))          (SETQ R X)          (COND             ((OR (< R -1.0)                  (> R 1.0))              (ERROR "ARCSIN: arg not in range" R)))          [COND             ((< R 0.0)                    (* * Range reduce to (0 1) via identity (ASIN          (minus X)) = (minus (ASIN X)) *)              (SETQ NEGATIVE T)              (SETQ R (- R]          [COND             ((> R .5)                    (* * Range reduce to (0 .5) via identity          (ASIN X) = (minus %%PI/2 (CL:* 2.0 (ASIN          (CL:SQRT (CL:* .5 (minus 1.0 R)))))) Avoids numerical instability calculating          (ASIN X) for X near one. SIN is horizontally flat near %%PI/2 so calculating          (ASIN X) by rational approximation wouldn't work well for X near          (SIN %%PI/2) = 1 *)              (SETQ REDUCED T)              (SETQ R (CL:SQRT (CL:* .5 (- 1.0 R]                    (* * R is now in range (0 .5) Use ARCSN 4671 rational approximation to           calculate (ASIN R) *)          [SETQ ANSWER (COND                          ((< R %%ASIN-EPSILON)                    (* * If R is in the interval (0 %%SIN-EPSILON) then          (ASIN R) = R to the precision that we can offer.          Return R because (1) it is desirable that          (ASIN R) = R exactly for small R and (2) microcode POLYEVAL will underflow on           sufficiently small positive R. *)                           R)                          (T (SETQ R2 (CL:* R R))                             (CL:* R (/ (POLYEVAL R2 %%ASIN-PPOLY 6)                                        (POLYEVAL R2 %%ASIN-QPOLY 6]          [COND             (REDUCED (SETQ ANSWER (- %%PI/2 (CL:* 2.0 ANSWER]          [COND             (NEGATIVE (SETQ ANSWER (- ANSWER]          [COND             (ACOS-FLAG                     (* * In case we want (ACOS X) then use identity          (ACOS X) = (minus %%PI/2 (ASIN X)) *)                    (SETQ ANSWER (- %%PI/2 ANSWER]          (RETURN ANSWER])(%%ASIN-COMPLEX  [CL:LAMBDA (Z)                                             (* kbr: "27-May-86 22:38")         (%%COMPLEX-MINUS (%%COMPLEX-TIMESI (CL:LOG (+ (%%COMPLEX-TIMESI Z)                                                       (CL:SQRT (- 1 (CL:* Z Z]))(* ACOS *)(DEFINEQ(ACOS  [CL:LAMBDA (X)                                             (* lmm "27-Jun-86 22:19")         (CTYPECASE X (FLOAT (%%ASIN-FLOAT X T))                (COMPLEX (%%ACOS-COMPLEX X))                (NUMBER (%%ASIN-FLOAT (FLOAT X)                               T])(%%ACOS-COMPLEX  [CL:LAMBDA (Z)                                             (* kbr: "27-May-86 22:31")         (%%COMPLEX-MINUS (%%COMPLEX-TIMESI (CL:LOG (+ Z (%%COMPLEX-TIMESI (CL:SQRT                                                                            (- 1 (CL:* Z Z]))(* ATAN *)(DEFINEQ(CL:ATAN  [CL:LAMBDA (X &OPTIONAL Y)                                 (* lmm "27-Jun-86 22:25")         (COND            (Y (%%ATAN-FLOAT (FLOAT X)                      (FLOAT Y)))            ((COMPLEXP X)             (%%ATAN-COMPLEX X))            (T (%%ATAN-FLOAT (FLOAT X])(%%ATAN1  [CL:LAMBDA (X)                                             (* (declare (inline short-atan)) *)         (COND            [(MINUSP X)             (- (%%ATAN1 (- X]            [(< X 2-SQRT3)             (- %%PI/2 (%%ATAN2 (/ X]            [(< X 1)             (- (+ %%PI/2 SIXTH-PI)                (%%ATAN2 (/ (+ SQRT3 X)                            (1- (CL:* X SQRT3]            ((< X INV-2-SQRT3)             (LET* ((INV (/ X)))                   (- [%%ATAN2 (/ (+ SQRT3 INV)                                  (1- (CL:* SQRT3 INV]                      SIXTH-PI)))            (T (%%ATAN2 X])(%%ATAN2  [CL:LAMBDA (X)         (CL:DO* [[SQR (- (/ (CL:* X X]                  (INT 1 (+ 2 INT))                  (OLD 0 VAL)                  (POW (/ X)                       (CL:* POW SQR))                  (VAL POW (+ VAL (/ POW INT]                ((= OLD VAL)                 (- %%PI/2 VAL])(%%ATAN-DOMAIN-CHECK  [CL:LAMBDA (Z)                                             (* kbr: "27-May-86 22:42")                                                             (* Return T if Z is in the domain of                                                              ATAN. ATAN is singular at i and -i.                                                             *)         (NOT (AND (ZEROP (REALPART Z))                   (= (ABS (IMAGPART Z))                      1])(%%ATAN-FLOAT  [CL:LAMBDA (Y &OPTIONAL X)         (COND            ((NOT X)             (CL:IF (ZEROP Y)                    0.0                    (%%ATAN1 Y)))            ((= Y X 0)             (CL:ERROR "Error in double entry atan both 0."))            ((= X 0)             (CL:* (%%SIGNUM Y)                   %%PI/2))            ((= Y 0)             (CL:IF (PLUSP X)                    0 %%PI))            ((AND (PLUSP X)                  (PLUSP Y))             (%%ATAN1 (/ Y X)))            [(PLUSP Y)             (- %%PI (%%ATAN1 (/ (- Y)                                 X]            [(PLUSP X)             (- (%%ATAN1 (/ Y (- X]            (T (- (%%ATAN1 (/ Y X))                  %%PI])(%%ATAN-COMPLEX  (CL:LAMBDA (Z)         (CL:IF (%%ATAN-DOMAIN-CHECK Z)                [LET ((I (%%COMPLEX-TIMESI Z)))                     (%%COMPLEX-MINUS (%%COMPLEX-TIMESI (CL:* .5 (CL:LOG (/ (+ 1 I)                                                                            (- 1 I]                (CL:ERROR "Argument not in domain for atan. ~S" Z)))))(DEFINEQ(CIS  [CL:LAMBDA (THETA)                                         (* Return cos (Theta) + i sin                                                             (Theta)%. *)         (CL:IF (COMPLEXP THETA)                (CL:ERROR "Argument to CIS is complex: ~S" THETA)                (COMPLEX (CL:COS THETA)                       (CL:SIN THETA]))(DEFINEQ(SINH  (CL:LAMBDA (X)                    (* Hyperbolic trig functions. Each of the hyperbolic trig functions is           calculated directly from their definition.          Exp (x) is calculated only once for efficiency.          They all work with complex arguments without modification.          *)         (LET ((Z (EXP X)))              (/ (- Z (/ Z))                 2))))(COSH  (CL:LAMBDA (X)         (LET ((Z (EXP X)))              (/ (+ Z (/ Z))                 2))))(TANH  [CL:LAMBDA (X)                                             (* Different form than in the manual.                                                             Does much better. *)         (LET* ((Z (EXP (CL:* 2 X)))                (Y (/ Z)))               (- (/ (1+ Y))                  (/ (1+ Z]))(DEFINEQ(ASINH  [CL:LAMBDA (X)                                             (* kbr: "27-May-86 22:43")         (CL:LOG (+ X (CL:SQRT (+ (CL:* X X)                                  1])(ACOSH  [CL:LAMBDA (X)                                             (* kbr: "27-May-86 22:45")         (CL:LOG (+ X (CL:SQRT (- (CL:* X X)                                  1])(ATANH  (CL:LAMBDA (X)                                             (* kbr: "27-May-86 22:46")         (CL:IF (%%ATANH-DOMAIN-CHECK X)                [CL:* .5 (CL:LOG (/ (1+ X)                                    (- 1 X]                (CL:ERROR "~S argument out of range." X))))(%%ATANH-DOMAIN-CHECK  [CL:LAMBDA (Z)                                             (* kbr: "27-May-86 22:49")                                                             (* Return T if Z is in the domain of                                                              ATANH. ATANH is singular at 1 and -1.0                                                              *)         (NOT (AND (ZEROP (IMAGPART Z))                   (= (ABS (REALPART Z))                      1]))(PUTPROPS CMLFLOAT FILETYPE COMPILE-FILE)(DECLARE: DONTEVAL@LOAD DOEVAL@COMPILE DONTCOPY COMPILERVARS (ADDTOVAR NLAMA )(ADDTOVAR NLAML )(ADDTOVAR LAMA           %%ATANH-DOMAIN-CHECK ATANH ACOSH ASINH TANH COSH SINH CIS %%ATAN-COMPLEX %%ATAN-FLOAT                 %%ATAN-DOMAIN-CHECK %%ATAN2 %%ATAN1 CL:ATAN %%ACOS-COMPLEX ACOS %%ASIN-COMPLEX ASIN                 %%TAN-COMPLEX CL:TAN %%COS-COMPLEX CL:COS %%SIN-COMPLEX CL:SIN CL:SQRT %%LOG-COMPLEX                 CL:LOG %%EXPT-COMPLEX-POWER %%EXPT-COMPLEX %%EXPT-INTEGER CL:EXPT %%EXP-COMPLEX                 %%EXP-FLOAT EXP))(PUTPROPS CMLFLOAT COPYRIGHT ("Xerox Corporation" 1986))STOP