(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