(FILECREATED "24-Jun-86 20:37:56" {ERIS}<LISPCORE>LIBRARY>CMLFLOAT.;30 71481
changes to: (FNS %%ATAN1 %%ATAN-FLOAT \FLOAT TRUNCATE %%COMPLEX-ABS COMPLEX CL:SQRT
%%SQRT-PERFECT %%SQRT-RATIO-PERFECT CL:SIN CL:COS %%OLD-TAN-FLOAT ASIN
%%OLD-ASIN-FLOAT ACOS)
(VARS CMLFLOATCOMS CMLFLOATJUNKCOMS)
previous date: "30-May-86 16:51:42" {ERIS}<LISPCORE>LIBRARY>CMLFLOAT.;27)
(* 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. -- By Kelly Roach. *)
(COMS (* LLFLOAT addition. *)
(FNS %%FLOAT %%MAKE-ARRAY))
(COMS * CMLFLOATJUNKCOMS)
(COMS (* Constants. *)
(CONSTANTS (PI (%%FLOAT 16457 4059))))
(COMS (* CMLARITH improvements. Put these back in CMLARITH when Larry gets finished. *)
(* PHASE and TRUNCATE had bugs I missed before handing off to Larry. %%COMPLEX-ABS
changes because CL:SQRT is smarter now and also %%COMPLEX-ABS was calling SQRT
instead of CL:SQRT. COMPLEX has been fixed so that (COMPLEX 0.0 0.0)
= 0.0 \FLOAT has been upgraded so that (FLOAT (COMPLEX 2 3))
=
(COMPLEX 2.0 3.0)
*)
(FNS PHASE TRUNCATE %%COMPLEX-ABS COMPLEX \FLOAT)
(P (printout T
"Warning: CMLFLOAT redefined some CMLARITH functions that need to be fixed."
T)))
(COMS (* EXP *)
(COMS * CMLFLOATEXPCOMS)
(FNS EXP %%EXP-FLOAT %%EXP-COMPLEX))
(COMS (* EXPT *)
(FNS CL:EXPT %%EXPT-INTEGER %%EXPT-FLOAT %%EXPT-COMPLEX %%EXPT-COMPLEX-POWER))
(COMS (* LOG *)
(COMS * CMLFLOATLOGCOMS)
(FNS CL:LOG %%LOG-FLOAT %%LOG-COMPLEX %%OLD-LOG-FLOAT))
(COMS (* SQRT *)
(FNS CL:SQRT %%SQRT-PERFECT %%SQRT-INTEGER %%SQRT-INTEGER-PERFECT %%SQRT-FLOAT
%%SQRT-RATIO %%SQRT-RATIO-PERFECT %%SQRT-COMPLEX %%SQRT-COMPLEX-PERFECT
%%OLD-SQRT-COMPLEX)
(FNS ISQRT))
(DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%2PI (%%FLOAT 16585 4059))
(%%PI (%%FLOAT 16457 4059))
(%%PI/2 (%%FLOAT 16329 4059))
(%%-PI/2 (%%FLOAT 49097 4059))
(%%PI/4 (%%FLOAT 16201 4059))
(%%-PI/4 (%%FLOAT 48969 4059))
(SIXTH-PI .5236)
(%%2/PI (%%FLOAT 16162 63875))))
(COMS (* SIN *)
(COMS * CMLFLOATSINCOMS)
(FNS CL:SIN %%SIN-FLOAT %%SIN-COMPLEX))
(COMS (* COS *)
(FNS CL:COS %%COS-COMPLEX))
(COMS (* TAN *)
(COMS * CMLFLOATTANCOMS)
(FNS CL:TAN %%TAN-FLOAT %%TAN-COMPLEX))
(COMS (* ASIN *)
(COMS * CMLFLOATASINCOMS)
(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))
(COMS (* CIS *)
(FNS CIS))
(COMS (* SINH COSH TANH *)
(FNS SINH COSH TANH))
(COMS (* ASINH ACOSH ATANH *)
(FNS ASINH ACOSH ATANH %%ATANH-DOMAIN-CHECK))
(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 ISQRT %%OLD-SQRT-COMPLEX %%SQRT-PERFECT CL:SQRT
%%LOG-COMPLEX CL:LOG %%EXPT-COMPLEX-POWER %%EXPT-COMPLEX %%EXPT-INTEGER
CL:EXPT %%EXP-COMPLEX %%EXP-FLOAT EXP COMPLEX %%COMPLEX-ABS TRUNCATE
PHASE)))))
(* * CMLFLOAT -- Covering sections 12.5-12.5.3 irrational, transcendental, exponential,
logarithmic, trigonometric, and hyperbolic functions. -- By Kelly Roach. *)
(* LLFLOAT addition. *)
(DEFINEQ
(%%FLOAT
(LAMBDA (HIWORD LOWORD) (* kbr: "14-May-86 16:43")
(\FLOATBOX (\VAG2 HIWORD LOWORD))))
(%%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))))
)
(RPAQQ CMLFLOATJUNKCOMS ((DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%SINGLE-FLOAT-EXPONENT-LENGTH
8)
(%%SINGLE-FLOAT-MANTISSA-LENGTH
23)
(%%FLOAT-ONE 1.0)
(%%FLOAT-HALF .5)
(%%SQRT-HALF (%%FLOAT 16181 1267))
(%%E (%%FLOAT 16429 63530))
(SQRT3 1.730525)
(2-SQRT3 .2679482)
(INV-2-SQRT3 3.732071)
(%%FLOAT-HALF .5)
(%%ZERO 0.0)))))
(DECLARE: EVAL@COMPILE DONTCOPY
(DECLARE: EVAL@COMPILE
(RPAQQ %%SINGLE-FLOAT-EXPONENT-LENGTH 8)
(RPAQQ %%SINGLE-FLOAT-MANTISSA-LENGTH 23)
(RPAQQ %%FLOAT-ONE 1.0)
(RPAQQ %%FLOAT-HALF .5)
(RPAQ %%SQRT-HALF (%%FLOAT 16181 1267))
(RPAQ %%E (%%FLOAT 16429 63530))
(RPAQQ SQRT3 1.730525)
(RPAQQ 2-SQRT3 .2679482)
(RPAQQ INV-2-SQRT3 3.732071)
(RPAQQ %%FLOAT-HALF .5)
(RPAQQ %%ZERO 0.0)
(CONSTANTS (%%SINGLE-FLOAT-EXPONENT-LENGTH 8)
(%%SINGLE-FLOAT-MANTISSA-LENGTH 23)
(%%FLOAT-ONE 1.0)
(%%FLOAT-HALF .5)
(%%SQRT-HALF (%%FLOAT 16181 1267))
(%%E (%%FLOAT 16429 63530))
(SQRT3 1.730525)
(2-SQRT3 .2679482)
(INV-2-SQRT3 3.732071)
(%%FLOAT-HALF .5)
(%%ZERO 0.0))
)
)
(* Constants. *)
(DECLARE: EVAL@COMPILE
(RPAQ PI (%%FLOAT 16457 4059))
(CONSTANTS (PI (%%FLOAT 16457 4059)))
)
(* CMLARITH improvements. Put these back in CMLARITH when Larry gets finished. *)
(* PHASE and TRUNCATE had bugs I missed before handing off to Larry. %%COMPLEX-ABS changes
because CL:SQRT is smarter now and also %%COMPLEX-ABS was calling SQRT instead of CL:SQRT.
COMPLEX has been fixed so that (COMPLEX 0.0 0.0) = 0.0 \FLOAT has been upgraded so that (FLOAT
(COMPLEX 2 3)) = (COMPLEX 2.0 3.0) *)
(DEFINEQ
(PHASE
(CL:LAMBDA (NUMBER) (* kbr: "14-May-86 11:26")
(* Returns the angle part of the polar
representation of NUMBER *)
(COND
((COMPLEXP NUMBER)
(CL:ATAN (IMAGPART NUMBER)
(REALPART NUMBER)))
((MINUSP NUMBER)
%%PI)
(T 0))))
(TRUNCATE
(CL:LAMBDA (NUMBER &OPTIONAL (DIVISOR 1)) (* kbr: "24-Jun-86 19:30")
(* Returns number (or number/divisor)
as an integer, rounded toward 0.0 The
second returned value is the
remainder. *)
(PROG (TRU REM)
(SETQ TRU (COND
((EQ DIVISOR 1)
(TYPECASE NUMBER (INTEGER NUMBER)
(FLOAT (\FIXP.FROM.FLOATP NUMBER))
(RATIO (IQUOTIENT (RATIO-NUMERATOR NUMBER)
(RATIO-DENOMINATOR NUMBER)))
(T (LISPERROR "ILLEGAL ARG" NUMBER))))
(T (TYPECASE NUMBER (INTEGER (TYPECASE DIVISOR (INTEGER (IQUOTIENT NUMBER
DIVISOR))
(FLOAT (FQUOTIENT NUMBER DIVISOR))
(RATIO (RETURN (TRUNCATE (/ NUMBER
DIVISOR))
))
(T (LISPERROR "ILLEGAL ARG" DIVISOR))
))
(FLOAT (TYPECASE DIVISOR ((OR INTEGER FLOAT)
(\FIXP.FROM.FLOATP (FQUOTIENT NUMBER
DIVISOR)))
(RATIO (RETURN (TRUNCATE (/ NUMBER DIVISOR))))
(T (LISPERROR "ILLEGAL ARG" DIVISOR))))
(RATIO (RETURN (TRUNCATE (/ NUMBER DIVISOR))))
(T (LISPERROR "ILLEGAL ARG" NUMBER))))))
(SETQ REM (- NUMBER (CL:* TRU DIVISOR)))
(RETURN (VALUES TRU REM)))))
(%%COMPLEX-ABS
(CL:LAMBDA (Z) (* kbr: "24-Jun-86 19:09")
(* Magnitude of a complex number *)
(PROG (X Y R A)
(SETQ X (COMPLEX-REALPART Z))
(SETQ Y (COMPLEX-IMAGPART Z))
(SETQ R (+ (CL:* X X)
(CL:* Y Y)))
(SETQ A (CL:SQRT R))
(RETURN A))))
(COMPLEX
(CL:LAMBDA (REALPART &OPTIONAL (IMAGPART 0)) (* kbr: "29-May-86 20:17")
(* Builds a complex number from the
specified components.
*)
(* TBW: Bug in TYPECASE? Had to change
RATIO to (OR INTEGER RATIO) in two
places. *)
(CL:IF (= IMAGPART 0)
(TYPECASE IMAGPART ((OR INTEGER RATIO)
REALPART)
(FLOAT (FLOAT REALPART))
(T (LISPERROR "ILLEGAL ARG" IMAGPART)))
(TYPECASE REALPART ((OR INTEGER RATIO)
(TYPECASE IMAGPART ((OR INTEGER RATIO)
(%%MAKE-COMPLEX REALPART IMAGPART))
(FLOAT (%%MAKE-COMPLEX (FLOAT REALPART)
IMAGPART))
(T (LISPERROR "ILLEGAL ARG" IMAGPART))))
(FLOAT (%%MAKE-COMPLEX REALPART (FLOAT IMAGPART)))
(T (LISPERROR "ILLEGAL ARG" REALPART))))))
(\FLOAT
(LAMBDA (X) (* kbr: "24-Jun-86 20:13")
(OR (FLOATP X)
(TYPECASE X (INTEGER
(SELECTC (NTYPX X)
(\FIXP (LET ((HI (fetch (FIXP HINUM) of X))
(LO (fetch (FIXP LONUM) of X))
(SIGN 0))
(COND
((IGREATERP HI MAX.POS.HINUM)
(.NEGATE. HI LO)
(SETQ SIGN 1)))
(\MAKEFLOAT SIGN (IPLUS \EXPONENT.BIAS 31)
HI LO T)))
(\SMALLP (LET* ((HI 0)
(SIGN 0)
(LO (COND
((IGEQ X 0)
X)
(T (SETQ SIGN 1)
(* X is negative--negate it)
(COND
((EQ 0 (\LOLOC X))
(* Min small integer)
(SETQ HI 1)
0)
(T (ADD1 (IDIFFERENCE MAX.SMALL.INTEGER
(\LOLOC X)))))))))
(\MAKEFLOAT SIGN (IPLUS \EXPONENT.BIAS 31)
HI LO T)))
(\BIGNUM.TO.FLOAT X)))
(RATIO (FQUOTIENT (RATIO-NUMERATOR X)
(RATIO-DENOMINATOR X)))
(COMPLEX (%%MAKE-COMPLEX (\FLOAT (REALPART X))
(\FLOAT (IMAGPART X))))
(T (\FLOAT (LISPERROR "NON-NUMERIC ARG" X T)))))))
)
(printout T "Warning: CMLFLOAT redefined some CMLARITH functions that need to be fixed." T)
(* EXP *)
(RPAQQ CMLFLOATEXPCOMS ((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)))))))
(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) (* kbr: "28-May-86 15:09")
(* Calculates e to the given power,
where e is the base of natural
logarithms. *)
(TYPECASE NUMBER (COMPLEX (%%EXP-COMPLEX NUMBER))
(NUMBER (%%EXP-FLOAT (FLOAT NUMBER)))
(T (LISPERROR "NON-NUMERIC ARG" 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) (* kbr: "28-May-86 15:12")
(* 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))
((ZEROP X)
(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 *)
(RPAQQ CMLFLOATLOGCOMS ((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)))))))
(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)) (* kbr: "30-May-86 00:19")
(COND
(BASE-SUPPLIED (/ (CL:LOG NUMBER)
(CL:LOG BASE)))
(T (TYPECASE NUMBER (COMPLEX (%%LOG-COMPLEX NUMBER))
(NUMBER (CL:IF (MINUSP NUMBER)
(COMPLEX (%%LOG-FLOAT (FLOAT (- NUMBER)))
%%PI)
(%%LOG-FLOAT (FLOAT NUMBER))))
(T (LISPERROR "NON-NUMERIC ARG" 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))))
(%%OLD-LOG-FLOAT
(LAMBDA (X) (* kbr: "29-May-86 00:02")
(* * (LOG X) for nonnegative float X. *)
(PROG (SX EXP SSUM)
(DECLARE (GLOBALVARS \LOGARRAY))
(COND
((NOT (FGREATERP X 0.0))
(ERROR "LOG OF ZERO:" X)))
(SETQ SX (\FLOAT.BOX X))
(SETQ EXP (IDIFFERENCE (fetch (FLOATP EXPONENT) of SX)
\EXPONENT.BIAS))
(* * 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 %%LOG-2 EXP))))))
)
(* SQRT *)
(DEFINEQ
(CL:SQRT
(CL:LAMBDA (X) (* kbr: "29-May-86 16:35")
(CL:IF (AND (NOT (COMPLEXP X))
(MINUSP X))
(COMPLEX 0 (CL:SQRT (ABS X)))
(TYPECASE X (INTEGER (%%SQRT-INTEGER X))
(FLOAT (%%SQRT-FLOAT X))
(RATIO (%%SQRT-RATIO X))
(COMPLEX (%%SQRT-COMPLEX X))
(T (LISPERROR "NON-NUMERIC ARG" X))))))
(%%SQRT-PERFECT
(CL:LAMBDA (X) (* kbr: "29-May-86 19:26")
(* * Perfect (SQRT X)%. Returns square root of X if X is a perfect square of a
rational number (includes integers) or complex number with rational
coefficients and otherwise NIL. *)
(CL:IF (AND (NOT (COMPLEXP X))
(MINUSP X))
(LET ((S (%%SQRT-PERFECT (ABS X))))
(CL:IF S (COMPLEX 0 S)))
(TYPECASE X (INTEGER (%%SQRT-INTEGER-PERFECT X))
(FLOAT NIL)
(RATIO (%%SQRT-RATIO-PERFECT X))
(COMPLEX (%%SQRT-COMPLEX-PERFECT X))
(T (LISPERROR "NON-NUMERIC ARG" X))))))
(%%SQRT-INTEGER
(LAMBDA (R) (* kbr: "29-May-86 17:35")
(* * (SQRT R) where R is a nonnegative integer.
*)
(* * A little bit of cleverness to return a ratio if R is a perfect square.
Otherwise returns a float. *)
(OR (%%SQRT-INTEGER-PERFECT R)
(%%SQRT-FLOAT (FLOAT R)))))
(%%SQRT-INTEGER-PERFECT
(LAMBDA (R) (* kbr: "29-May-86 17:39")
(* * Perfect (SQRT R) where R is a nonnegative integer.
If R is a perfect square, return square root of R, otherwise return NIL.
*)
(PROG (A)
(SETQ A (ISQRT R))
(COND
((NOT (EQ (CL:* A A)
R))
(RETURN NIL)))
(RETURN A))))
(%%SQRT-FLOAT
(LAMBDA (X) (* kbr: "29-May-86 20:32")
(* * (SQRT X) for nonnegative float X. *)
(PROG (V)
(DECLARE (TYPE FLOATP X V))
(COND
((NOT (FGREATERP X 0.0)) (* Trichotomy ==> X = 0.0)
(RETURN 0.0)))
(SETQ V (create FLOATP
EXPONENT ← (LOGAND (IPLUS \EXPONENT.BIAS
(LRSH (LOGAND (IDIFFERENCE (fetch (FLOATP EXP)
of X)
\EXPONENT.BIAS)
(MASK.1'S 0 BITSPERWORD))
1))
\MAX.EXPONENT)
HIFRACTION ← (fetch (FLOATP HIFRAC) of X)))
(* Exponent is stored as excess \EXPONENT.BIAS and although the LRSH doesn't
really do division by 2 (e.g., when the arg is negative) at least the low-order
8 bits will be right. It doesn't even matter that it may be off-by-one, due to
the infamous "Arithmetic Shifting Considered Harmful" since it is only an
estimate.)
(FRPTQ 4 (SETQ V (CL:* .5 (+ V (/ X V)))))
(RETURN V))))
(%%SQRT-RATIO
(LAMBDA (R) (* kbr: "29-May-86 17:41")
(* * (SQRT R) where R is a nonnegative ratio.
*)
(* * A little bit of cleverness to return a ratio if R is a perfect square.
Otherwise returns a float. *)
(OR (%%SQRT-RATIO-PERFECT R)
(%%SQRT-FLOAT (FLOAT R)))))
(%%SQRT-RATIO-PERFECT
(LAMBDA (R) (* kbr: "29-May-86 17:40")
(* * Perfect (SQRT R) where R is a nonnegative ratio.
If R is a perfect square, return square root of R, otherwise NIL.
*)
(PROG (RNUM RDEN ANUM ADEN)
(* * A little bit of cleverness to return a ratio if R is a perfect square.
Otherwise returns a float. *)
(SETQ RNUM (RATIO-NUMERATOR R))
(SETQ RDEN (RATIO-DENOMINATOR R))
(SETQ ANUM (%%SQRT-INTEGER-PERFECT RNUM))
(COND
((NULL ANUM)
(RETURN NIL)))
(SETQ ADEN (%%SQRT-INTEGER-PERFECT RDEN))
(COND
((NULL ADEN)
(RETURN NIL)))
(RETURN (%%MAKE-RATIO ANUM ADEN)))))
(%%SQRT-COMPLEX
(LAMBDA (X) (* kbr: "29-May-86 21:16")
(* * (SQRT X) for complex X. *)
(PROG (ABS PHASE A B C D E ANSWER)
(DECLARE (TYPE FLOATP ABS PHASE A B C D E))
(SETQ ANSWER (%%SQRT-COMPLEX-PERFECT X))
(COND
(ANSWER (RETURN ANSWER)))
(SETQ A (FLOAT (REALPART X)))
(SETQ B (FLOAT (IMAGPART X)))
(* * Make initial guess. *)
(SETQ ABS (CL:SQRT (ABS X)))
(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))))
(%%SQRT-COMPLEX-PERFECT
(LAMBDA (R) (* kbr: "29-May-86 20:06")
(* * Perfect (SQRT R) where R is a complex number.
If R is a perfect square, return square root of R, otherwise NIL.
*)
(PROG (A B M N DELTA)
(* * First of all, let's say that R = M+iN.
*)
(SETQ M (REALPART R))
(COND
((FLOATP M)
(RETURN NIL)))
(SETQ N (IMAGPART R))
(* * If M is 0, then R is a pure imaginary and life is simple.
*)
(COND
((= M 0)
(COND
((PLUSP N) (* Looking for R = Ni =
(CL:* 2 (EXPT A 2) i) =
(CL:* A+Ai A+Ai) *)
(SETQ A (%%SQRT-PERFECT (/ N 2)))
(COND
((NULL A)
(RETURN NIL)))
(RETURN (%%MAKE-COMPLEX A A)))
(T (* N must be negative because we
cannot have M and N = 0
simultaneously. *)
(* Looking for R = Ni =
(CL:* -2 (EXPT A 2) i) =
(CL:* A-Ai A-Ai) *)
(SETQ A (%%SQRT-PERFECT (/ N -2)))
(COND
((NULL A)
(RETURN NIL)))
(RETURN (%%MAKE-COMPLEX A (- A)))))))
(* * At this point neither M or N is 0 %.
We need to solve for A+Bi such that (EXPT A+Bi 2) = M+iN implies M =
(minus (EXPT A 2) (EXPT B 2)) and N = (CL:* 2 A B) These requirements lead to a
quartic with four possible solutions. Two of the solutions are eliminated by
the requirement that B must be real and the third solution is eliminated by its
being on the wrong branch of SQRT. I leave the development of this quartic and
its solution as an exercise for the reader.
*)
(SETQ DELTA (%%SQRT-PERFECT (+ (CL:* M M)
(CL:* N N))))
(COND
((NULL DELTA)
(RETURN NIL)))
(SETQ A (%%SQRT-PERFECT (/ (+ M DELTA)
2)))
(COND
((NULL A)
(RETURN NIL)))
(SETQ B (/ N (CL:* 2 A)))
(RETURN (%%MAKE-COMPLEX A B)))))
(%%OLD-SQRT-COMPLEX
(CL:LAMBDA (Z) (* kbr: "29-May-86 19:42")
(* The square root of a complex number
returns two roots. This function
returns the primary root.
*)
(OR (%%SQRT-COMPLEX-PERFECT Z)
(LET ((ABS (CL:SQRT (ABS Z)))
(PHASE (/ (PHASE Z)
2.0)))
(COMPLEX (CL:* ABS (CL:COS PHASE))
(CL:* ABS (CL:SIN PHASE)))))))
)
(DEFINEQ
(ISQRT
(CL:LAMBDA (N) (* kbr: "27-May-86 20:40")
(* ISQRT: Integer square root -
isqrt (n) **2 <= n Upper and lower bounds on the result are estimated using
integer-length. On each iteration, one of the bounds is replaced by their mean.
The lower bound is returned when the bounds meet or differ by only 1.0 Initial
bounds guarantee that lg (sqrt (n)) = lg
(n) /2 iterations suffice. *)
(CL:IF (AND (INTEGERP N)
(NOT (MINUSP N)))
(CL:DO* ((LG (INTEGER-LENGTH N))
(LO (ASH 1 (ASH (1- LG)
-1)))
(HI (+ LO (ASH LO (CL:IF (ODDP LG)
-1 0)))))
((<= (1- HI)
LO)
LO)
(LET ((MID (ASH (+ LO HI)
-1)))
(CL:IF (<= (CL:* MID MID)
N)
(SETQ LO MID)
(SETQ HI MID))))
(CL:ERROR "ISQRT: ~S argument must be a nonnegative integer." N))))
)
(DECLARE: EVAL@COMPILE DONTCOPY
(DECLARE: EVAL@COMPILE
(RPAQ %%2PI (%%FLOAT 16585 4059))
(RPAQ %%PI (%%FLOAT 16457 4059))
(RPAQ %%PI/2 (%%FLOAT 16329 4059))
(RPAQ %%-PI/2 (%%FLOAT 49097 4059))
(RPAQ %%PI/4 (%%FLOAT 16201 4059))
(RPAQ %%-PI/4 (%%FLOAT 48969 4059))
(RPAQQ SIXTH-PI .5236)
(RPAQ %%2/PI (%%FLOAT 16162 63875))
(CONSTANTS (%%2PI (%%FLOAT 16585 4059))
(%%PI (%%FLOAT 16457 4059))
(%%PI/2 (%%FLOAT 16329 4059))
(%%-PI/2 (%%FLOAT 49097 4059))
(%%PI/4 (%%FLOAT 16201 4059))
(%%-PI/4 (%%FLOAT 48969 4059))
(SIXTH-PI .5236)
(%%2/PI (%%FLOAT 16162 63875)))
)
)
(* SIN *)
(RPAQQ CMLFLOATSINCOMS ((* * %%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 (CONSTANTS (%%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)) . *)
(INITVARS (%%SIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 50193 50772)
(%%FLOAT 18371 47202)
(%%FLOAT 51905 54216)
(%%FLOAT 19748 29063)
(%%FLOAT 52951 41140)
(%%FLOAT 20368 50691))))
(%%SIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 16256 0)
(%%FLOAT 17252 49821)
(%%FLOAT 18142 34466)
(%%FLOAT 18957 45135)
(%%FLOAT 19685 22067)
(%%FLOAT 20280 21713)))))))
(* * %%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
(DECLARE: EVAL@COMPILE
(RPAQ %%SIN-EPSILON (%%FLOAT 14720 0))
(CONSTANTS (%%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)) . *)
(RPAQ? %%SIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 50193 50772)
(%%FLOAT 18371 47202)
(%%FLOAT 51905 54216)
(%%FLOAT 19748 29063)
(%%FLOAT 52951 41140)
(%%FLOAT 20368 50691))))
(RPAQ? %%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")
(CL:IF (AND (NOT (FLOATP X))
(NOT (COMPLEXP X)))
(SETQ X (FLOAT X)))
(TYPECASE X (FLOAT (%%SIN-FLOAT X NIL))
(COMPLEX (%%SIN-COMPLEX X)))))
(%%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) (* kbr: "31-Mar-86 17:34")
(CL:IF (AND (NOT (FLOATP X))
(NOT (COMPLEXP X)))
(SETQ X (FLOAT X)))
(TYPECASE X (FLOAT (%%SIN-FLOAT X T))
(COMPLEX (%%COS-COMPLEX X)))))
(%%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 *)
(RPAQQ CMLFLOATTANCOMS ((* * %%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 (CONSTANTS (%%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) . *)
(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)))))))
(* * %%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
(DECLARE: EVAL@COMPILE
(RPAQ %%TAN-EPSILON (%%FLOAT 14720 0))
(CONSTANTS (%%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) (* kbr: "28-May-86 20:56")
(TYPECASE X (COMPLEX (%%TAN-COMPLEX X))
(NUMBER (%%TAN-FLOAT (FLOAT X)))
(T (LISPERROR "NON-NUMERIC ARG" 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 *)
(RPAQQ CMLFLOATASINCOMS ((* * %%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 (CONSTANTS (%%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)) . *)
(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)))))))
(* * %%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
(DECLARE: EVAL@COMPILE
(RPAQ %%ASIN-EPSILON (%%FLOAT 14720 0))
(CONSTANTS (%%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) (* kbr: "30-May-86 16:49")
(TYPECASE X (FLOAT (%%ASIN-FLOAT X))
(COMPLEX (%%ASIN-COMPLEX X))
(NUMBER (ASIN (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) (* kbr: "30-May-86 16:49")
(TYPECASE X (FLOAT (%%ASIN-FLOAT X T))
(COMPLEX (%%ACOS-COMPLEX X))
(NUMBER (ACOS (FLOAT X))))))
(%%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) (* kbr: "14-May-86 11:41")
(COND
(Y (CL:IF (OR (COMPLEXP X)
(COMPLEXP Y))
(CL:ERROR "~S and ~S both complex." X Y))
(%%ATAN-FLOAT (FLOAT X)
(FLOAT Y)))
((COMPLEXP X)
(%%ATAN-COMPLEX X))
(T (%%ATAN-FLOAT (FLOAT X))))))
(%%ATAN1
(CL:LAMBDA (X) (* kbr: "24-Jun-86 20:35")
(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) (* kbr: "24-Jun-86 20:37")
(COND
((NULL 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))))
)
(* CIS *)
(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)))))
)
(* SINH COSH TANH *)
(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))))))
)
(* ASINH ACOSH ATANH *)
(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)))))
)
(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 ISQRT
%%OLD-SQRT-COMPLEX %%SQRT-PERFECT CL:SQRT %%LOG-COMPLEX CL:LOG %%EXPT-COMPLEX-POWER
%%EXPT-COMPLEX %%EXPT-INTEGER CL:EXPT %%EXP-COMPLEX %%EXP-FLOAT EXP COMPLEX
%%COMPLEX-ABS TRUNCATE PHASE)
)
(PUTPROPS CMLFLOAT COPYRIGHT ("Xerox Corporation" 1986))
(DECLARE: DONTCOPY
(FILEMAP (NIL (4760 5764 (%%FLOAT 4770 . 4918) (%%MAKE-ARRAY 4920 . 5762)) (8201 15355 (PHASE 8211 .
8714) (TRUNCATE 8716 . 11182) (%%COMPLEX-ABS 11184 . 11659) (COMPLEX 11661 . 13148) (\FLOAT 13150 .
15353)) (18239 20654 (EXP 18249 . 18817) (%%EXP-FLOAT 18819 . 20306) (%%EXP-COMPLEX 20308 . 20652)) (
20674 24425 (CL:EXPT 20684 . 22053) (%%EXPT-INTEGER 22055 . 23229) (%%EXPT-FLOAT 23231 . 23541) (
%%EXPT-COMPLEX 23543 . 23976) (%%EXPT-COMPLEX-POWER 23978 . 24423)) (26810 30513 (CL:LOG 26820 . 27445
) (%%LOG-FLOAT 27447 . 29212) (%%LOG-COMPLEX 29214 . 29425) (%%OLD-LOG-FLOAT 29427 . 30511)) (30533
40762 (CL:SQRT 30543 . 31066) (%%SQRT-PERFECT 31068 . 31899) (%%SQRT-INTEGER 31901 . 32331) (
%%SQRT-INTEGER-PERFECT 32333 . 32819) (%%SQRT-FLOAT 32821 . 34349) (%%SQRT-RATIO 34351 . 34775) (
%%SQRT-RATIO-PERFECT 34777 . 35640) (%%SQRT-COMPLEX 35642 . 37028) (%%SQRT-COMPLEX-PERFECT 37030 .
39977) (%%OLD-SQRT-COMPLEX 39979 . 40760)) (40763 42124 (ISQRT 40773 . 42122)) (46631 49570 (CL:SIN
46641 . 46971) (%%SIN-FLOAT 46973 . 49204) (%%SIN-COMPLEX 49206 . 49568)) (49589 50303 (CL:COS 49599
. 49927) (%%COS-COMPLEX 49929 . 50301)) (54002 57128 (CL:TAN 54012 . 54284) (%%TAN-FLOAT 54286 .
56600) (%%TAN-COMPLEX 56602 . 57126)) (61390 64391 (ASIN 61400 . 61657) (%%ASIN-FLOAT 61659 . 64100) (
%%ASIN-COMPLEX 64102 . 64389)) (64411 64992 (ACOS 64421 . 64680) (%%ACOS-COMPLEX 64682 . 64990)) (
65012 68168 (CL:ATAN 65022 . 65476) (%%ATAN1 65478 . 66131) (%%ATAN2 66133 . 66457) (
%%ATAN-DOMAIN-CHECK 66459 . 66961) (%%ATAN-FLOAT 66963 . 67784) (%%ATAN-COMPLEX 67786 . 68166)) (68187
68573 (CIS 68197 . 68571)) (68603 69482 (SINH 68613 . 69034) (COSH 69036 . 69150) (TANH 69152 . 69480
)) (69515 70747 (ASINH 69525 . 69726) (ACOSH 69728 . 69929) (ATANH 69931 . 70236) (
%%ATANH-DOMAIN-CHECK 70238 . 70745)))))
STOP