(FILECREATED "28-Jul-86 14:12:13" {ERIS}<ROACH>CML>CMLFLOAT.;21 80481
changes to: (VARS CMLFLOATCOMS)
(FNS CL:SQRT %%SQRT-COMPLEX)
previous date: "23-Jul-86 18:31:49" {ERIS}<ROACH>CML>CMLFLOAT.;20)
(* 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. Section 12.10, implementation
parameters. -- By Kelly Roach. *)
(COMS (* Section 12.10, implementation parameters. The constants in this COMS are exported to
the user. *)
(* %%FLOAT allows us to recreate FLOATPs in a way that is independent of the ordinairy
reading and printing FLOATPs to files which involves loss of the last couple bits of
accuracy due to rounding effects. *)
(FNS %%FLOAT)
(CONSTANTS (MOST-POSITIVE-FIXNUM 65535)
(MOST-NEGATIVE-FIXNUM -65536)
(MOST-POSITIVE-SINGLE-FLOAT (%%FLOAT 32639 65535))
(LEAST-POSITIVE-SINGLE-FLOAT (%%FLOAT 128 0))
(LEAST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 32896 0))
(MOST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 65407 65535))
(MOST-POSITIVE-SHORT-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-SHORT-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-SHORT-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-SHORT-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(MOST-POSITIVE-DOUBLE-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-DOUBLE-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-DOUBLE-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-DOUBLE-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(MOST-POSITIVE-LONG-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-LONG-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-LONG-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-LONG-FLOAT MOST-NEGATIVE-SINGLE-FLOAT))
(* EPSILON is the smallest positive floating point number such that
(NOT (= (FLOAT 1 EPSILON)
(+ (FLOAT 1 EPSILON)
EPSILON)))
*)
(CONSTANTS (SINGLE-FLOAT-EPSILON (%%FLOAT 13312 0))
(SHORT-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(DOUBLE-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(LONG-FLOAT-EPSILON SINGLE-FLOAT-EPSILON))
(* NEGATIVE-EPSILON is the smallest negative floating point number such that
(NOT (= (FLOAT 1 NEGATIVE-EPSILON)
(- (FLOAT 1 NEGATIVE-EPSILON)
NEGATIVE-EPSILON)))
*)
(CONSTANTS (SINGLE-FLOAT-NEGATIVE-EPSILON (%%FLOAT 13184 0))
(SHORT-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(DOUBLE-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(LONG-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)))
(COMS (* Miscellaneous. *)
(DECLARE: EVAL@COMPILE DONTCOPY (FILES (LOADCOMP)
LLFLOAT))
(CONSTANTS (PI (%%FLOAT 16457 4059)))
(DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%E (%%FLOAT 16429 63572))
(%%2PI (%%FLOAT 16585 4059))
(%%PI (%%FLOAT 16457 4059))
(%%2PI/3 (%%FLOAT 16390 2706))
(%%PI/2 (%%FLOAT 16329 4059))
(%%-PI/2 (%%FLOAT 49097 4059))
(%%PI/3 (%%FLOAT 16262 2706))
(%%PI/4 (%%FLOAT 16201 4059))
(%%-PI/4 (%%FLOAT 48969 4059))
(%%PI/6 (%%FLOAT 16134 2706))
(%%2/PI (%%FLOAT 16162 63875))))
(FNS %%MAKE-ARRAY))
(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)) . *)
(VARS (%%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)) . *)
(VARS (%%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))
(COMS (* SQRT *)
(FNS CL:SQRT %%SQRT-FLOAT %%SQRT-COMPLEX))
(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 (CONSTANTS (%%SIN-EPSILON (%%FLOAT 14720 0))))
(* * %%SIN-PPOLY and %%SIN-QPOLY contain adapted P and Q coefficients of Harris
et al SIN 3374 rational approximation to (SIN X)
in interval (0 (/ PI 2))
%. The coefficients for %%SIN-PPOLY and %%SIN-QPOLY have been computed from
Harris using extended precision routines and the relations %%SIN-PPOLY =
(REVERSE (for I from 0 as ENTRY in PS collect
(/ (CL:* (EXPT (/ 2 PI)
(1+ (CL:* 2 I)))
ENTRY)
Q0)))
and %%SIN-QPOLY = (REVERSE (for I from 0 as ENTRY in QS collect
(/ (CL:* (EXPT (/ 2 PI)
(CL:* 2 I))
ENTRY)
Q0)))
*)
(VARS (%%SIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 45236 25611)
(%%FLOAT 13589 26148)
(%%FLOAT 47286 34797)
(%%FLOAT 15295 3306)
(%%FLOAT 48666 34805)
(%%FLOAT 16256 0))))
(%%SIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 11384 52865)
(%%FLOAT 12553 9550)
(%%FLOAT 13604 38385)
(%%FLOAT 14593 18841)
(%%FLOAT 15489 5549)
(%%FLOAT 16256 0))))))
(FNS CL:SIN %%SIN-FLOAT %%SIN-COMPLEX))
(COMS (* COS *)
(FNS CL:COS %%COS-COMPLEX))
(COMS (* TAN *)
(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 (CONSTANTS (%%TAN-EPSILON (%%FLOAT 14720 0))))
(* * %%TAN-PPOLY and %%TAN-QPOLY contain adapted P and Q coefficients of Harris
et al TAN 4288 rational approximation to (TAN X)
in interval (-PI/4 PI/4)
%. The coefficients for %%TAN-PPOLY and %%TAN-QPOLY have been computed from
Harris using extended precision routines and the relations %%TAN-PPOLY =
(REVERSE (for I from 0 as ENTRY in PS collect
(/ (CL:* (EXPT (/ 4 PI)
(1+ (CL:* 2 I)))
ENTRY)
Q0)))
and %%TAN-QPOLY = (REVERSE (for I from 0 as ENTRY in QS collect
(/ (CL:* (EXPT (/ 4 PI)
(CL:* 2 I))
ENTRY)
Q0)))
*)
(VARS (%%TAN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 13237 21090)
(%%FLOAT 47141 15825)
(%%FLOAT 15246 8785)
(%%FLOAT 48655 48761)
(%%FLOAT 16256 0))))
(%%TAN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 45267 36947)
(%%FLOAT 13848 46875)
(%%FLOAT 47612 53738)
(%%FLOAT 15596 52854)
(%%FLOAT 48882 35303)
(%%FLOAT 16256 0))))))
(FNS CL:TAN %%TAN-FLOAT %%TAN-COMPLEX))
(COMS (* ASIN *)
(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 (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)) . *)
(VARS (%%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 *)
(DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%SQRT3 (%%FLOAT 16349 46039))
(%%2-SQRT3 (%%FLOAT 16009 12451))
(%%INV-2-SQRT3 (%%FLOAT 16494 55788)))
(FILES (LOADCOMP)
LLFLOAT))
(FNS CL:ATAN %%ATAN-FLOAT1 %%ATAN-FLOAT2 %%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))
(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 %%ATAN-FLOAT2 %%ATAN-FLOAT1 CL:ATAN
%%ACOS-COMPLEX ACOS %%ASIN-COMPLEX ASIN %%TAN-COMPLEX CL:TAN
%%COS-COMPLEX CL:COS %%SIN-COMPLEX CL:SIN %%SQRT-PERFECT 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. Section 12.10, implementation parameters.
-- By Kelly Roach. *)
(* Section 12.10, implementation parameters. The constants in this COMS are exported to the
user. *)
(* %%FLOAT allows us to recreate FLOATPs in a way that is independent of the ordinairy reading
and printing FLOATPs to files which involves loss of the last couple bits of accuracy due to
rounding effects. *)
(DEFINEQ
(%%FLOAT
(LAMBDA (HIWORD LOWORD) (* kbr: "14-May-86 16:43")
(\FLOATBOX (\VAG2 HIWORD LOWORD))))
)
(DECLARE: EVAL@COMPILE
(RPAQQ MOST-POSITIVE-FIXNUM 65535)
(RPAQQ MOST-NEGATIVE-FIXNUM -65536)
(RPAQ MOST-POSITIVE-SINGLE-FLOAT (%%FLOAT 32639 65535))
(RPAQ LEAST-POSITIVE-SINGLE-FLOAT (%%FLOAT 128 0))
(RPAQ LEAST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 32896 0))
(RPAQ MOST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 65407 65535))
(RPAQ MOST-POSITIVE-SHORT-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(RPAQ LEAST-POSITIVE-SHORT-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(RPAQ LEAST-NEGATIVE-SHORT-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(RPAQ MOST-NEGATIVE-SHORT-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(RPAQ MOST-POSITIVE-DOUBLE-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(RPAQ LEAST-POSITIVE-DOUBLE-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(RPAQ LEAST-NEGATIVE-DOUBLE-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(RPAQ MOST-NEGATIVE-DOUBLE-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(RPAQ MOST-POSITIVE-LONG-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(RPAQ LEAST-POSITIVE-LONG-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(RPAQ LEAST-NEGATIVE-LONG-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(RPAQ MOST-NEGATIVE-LONG-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(CONSTANTS (MOST-POSITIVE-FIXNUM 65535)
(MOST-NEGATIVE-FIXNUM -65536)
(MOST-POSITIVE-SINGLE-FLOAT (%%FLOAT 32639 65535))
(LEAST-POSITIVE-SINGLE-FLOAT (%%FLOAT 128 0))
(LEAST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 32896 0))
(MOST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 65407 65535))
(MOST-POSITIVE-SHORT-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-SHORT-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-SHORT-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-SHORT-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(MOST-POSITIVE-DOUBLE-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-DOUBLE-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-DOUBLE-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-DOUBLE-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(MOST-POSITIVE-LONG-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-LONG-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-LONG-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-LONG-FLOAT MOST-NEGATIVE-SINGLE-FLOAT))
)
(* EPSILON is the smallest positive floating point number such that (NOT (= (FLOAT 1 EPSILON) (
+ (FLOAT 1 EPSILON) EPSILON))) *)
(DECLARE: EVAL@COMPILE
(RPAQ SINGLE-FLOAT-EPSILON (%%FLOAT 13312 0))
(RPAQ SHORT-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(RPAQ DOUBLE-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(RPAQ LONG-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(CONSTANTS (SINGLE-FLOAT-EPSILON (%%FLOAT 13312 0))
(SHORT-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(DOUBLE-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(LONG-FLOAT-EPSILON SINGLE-FLOAT-EPSILON))
)
(* NEGATIVE-EPSILON is the smallest negative floating point number such that (NOT (= (FLOAT 1
NEGATIVE-EPSILON) (- (FLOAT 1 NEGATIVE-EPSILON) NEGATIVE-EPSILON))) *)
(DECLARE: EVAL@COMPILE
(RPAQ SINGLE-FLOAT-NEGATIVE-EPSILON (%%FLOAT 13184 0))
(RPAQ SHORT-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(RPAQ DOUBLE-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(RPAQ LONG-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(CONSTANTS (SINGLE-FLOAT-NEGATIVE-EPSILON (%%FLOAT 13184 0))
(SHORT-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(DOUBLE-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(LONG-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON))
)
(* Miscellaneous. *)
(DECLARE: EVAL@COMPILE DONTCOPY
(FILESLOAD (LOADCOMP)
LLFLOAT)
)
(DECLARE: EVAL@COMPILE
(RPAQ PI (%%FLOAT 16457 4059))
(CONSTANTS (PI (%%FLOAT 16457 4059)))
)
(DECLARE: EVAL@COMPILE DONTCOPY
(DECLARE: EVAL@COMPILE
(RPAQ %%E (%%FLOAT 16429 63572))
(RPAQ %%2PI (%%FLOAT 16585 4059))
(RPAQ %%PI (%%FLOAT 16457 4059))
(RPAQ %%2PI/3 (%%FLOAT 16390 2706))
(RPAQ %%PI/2 (%%FLOAT 16329 4059))
(RPAQ %%-PI/2 (%%FLOAT 49097 4059))
(RPAQ %%PI/3 (%%FLOAT 16262 2706))
(RPAQ %%PI/4 (%%FLOAT 16201 4059))
(RPAQ %%-PI/4 (%%FLOAT 48969 4059))
(RPAQ %%PI/6 (%%FLOAT 16134 2706))
(RPAQ %%2/PI (%%FLOAT 16162 63875))
(CONSTANTS (%%E (%%FLOAT 16429 63572))
(%%2PI (%%FLOAT 16585 4059))
(%%PI (%%FLOAT 16457 4059))
(%%2PI/3 (%%FLOAT 16390 2706))
(%%PI/2 (%%FLOAT 16329 4059))
(%%-PI/2 (%%FLOAT 49097 4059))
(%%PI/3 (%%FLOAT 16262 2706))
(%%PI/4 (%%FLOAT 16201 4059))
(%%-PI/4 (%%FLOAT 48969 4059))
(%%PI/6 (%%FLOAT 16134 2706))
(%%2/PI (%%FLOAT 16162 63875)))
)
)
(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))))
)
(* 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: "23-Jul-86 13:45")
(* * (EXP X) for float X calculated via EXPB 1103 rational approximation of
Harris et al. *)
(PROG (RECIPFLG M N R ANSWER)
(DECLARE (GLOBALVARS %%EXP-TABLE %%EXP-POLY))
(* * First, arrange X to be in interval (0 infinity) via identity
(EXP (minus X)) = (/ (EXP X)) *)
(SETQ R X)
(COND
((< R 0)
(SETQ R (- R))
(SETQ RECIPFLG T)))
(* * Next, 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)
(TRUNCATE (CL:* %%LOG-BASE2-E R)))
(MULTIPLE-VALUE-SETQ (N R)
(TRUNCATE R .125))
(SETQ ANSWER (SCALE-FLOAT (CL:* (ELT %%EXP-TABLE N)
(/ (POLYEVAL R %%EXP-POLY 5)
(POLYEVAL (- R)
%%EXP-POLY 5)))
M))
(COND
(RECIPFLG (SETQ ANSWER (/ ANSWER))))
(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 (BASE-NUMBER POWER-NUMBER) (* lmm "27-Jun-86 22:52")
(* This function calculates BASE-NUMBER raised to the nth power.
It separates the cases by the type of POWER-NUMBER for efficiency reasons, as
powers can be calculated more efficiently if POWER-NUMBER is a positive
integer, Therefore, All integer values of POWER-NUMBER are calculated as
positive integers, and inverted if negative.
*)
(COND
((AND (RATIONALP BASE-NUMBER)
(INTEGERP POWER-NUMBER))
(%%EXPT-INTEGER BASE-NUMBER POWER-NUMBER))
((= BASE-NUMBER 0)
(CL:IF (PLUSP POWER-NUMBER)
BASE-NUMBER
(CL:ERROR "~A to a non-positive power ~A." BASE-NUMBER POWER-NUMBER)))
((AND (COMPLEXP BASE-NUMBER)
(INTEGERP POWER-NUMBER))
(%%EXPT-INTEGER BASE-NUMBER POWER-NUMBER))
((COMPLEXP POWER-NUMBER)
(CL:IF (OR (COMPLEXP BASE-NUMBER)
(PLUSP BASE-NUMBER))
(%%EXPT-COMPLEX-POWER BASE-NUMBER POWER-NUMBER)
(CL:ERROR "~A negative number to a complex power ~A." BASE-NUMBER POWER-NUMBER)))
((COMPLEXP BASE-NUMBER)
(%%EXPT-COMPLEX BASE-NUMBER POWER-NUMBER))
((AND (NOT (INTEGERP POWER-NUMBER))
(MINUSP BASE-NUMBER))
(CL:ERROR "Negative number ~A to non-integral power ~A." BASE-NUMBER POWER-NUMBER))
(T (%%EXPT-FLOAT (FLOAT BASE-NUMBER)
(FLOAT POWER-NUMBER))))))
(%%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))))
)
(* SQRT *)
(DEFINEQ
(CL:SQRT
(CL:LAMBDA (X) (* kbr: "28-Jul-86 14:06")
(COND
((COMPLEXP X)
(%%SQRT-COMPLEX X))
((MINUSP X) (* Negative real axis maps into
positive imaginary axis.
*)
(COMPLEX 0 (CL:SQRT (- X))))
(T (%%SQRT-FLOAT (FLOAT X))))))
(%%SQRT-FLOAT
(LAMBDA (X) (* kbr: "21-Jul-86 16:38")
(* * (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-COMPLEX
(LAMBDA (X) (* kbr: "28-Jul-86 14:04")
(* * (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 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))))
)
(* 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
(DECLARE: EVAL@COMPILE
(RPAQ %%SIN-EPSILON (%%FLOAT 14720 0))
(CONSTANTS (%%SIN-EPSILON (%%FLOAT 14720 0)))
)
)
(* * %%SIN-PPOLY and %%SIN-QPOLY contain adapted P and Q coefficients of Harris et al SIN 3374
rational approximation to (SIN X) in interval (0 (/ PI 2)) %. The coefficients for %%SIN-PPOLY
and %%SIN-QPOLY have been computed from Harris using extended precision routines and the
relations %%SIN-PPOLY = (REVERSE (for I from 0 as ENTRY in PS collect (/ (CL:* (EXPT (/ 2 PI) (
1+ (CL:* 2 I))) ENTRY) Q0))) and %%SIN-QPOLY = (REVERSE (for I from 0 as ENTRY in QS collect (/
(CL:* (EXPT (/ 2 PI) (CL:* 2 I)) ENTRY) Q0))) *)
(RPAQ %%SIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 45236 25611)
(%%FLOAT 13589 26148)
(%%FLOAT 47286 34797)
(%%FLOAT 15295 3306)
(%%FLOAT 48666 34805)
(%%FLOAT 16256 0))))
(RPAQ %%SIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 11384 52865)
(%%FLOAT 12553 9550)
(%%FLOAT 13604 38385)
(%%FLOAT 14593 18841)
(%%FLOAT 15489 5549)
(%%FLOAT 16256 0))))
(DEFINEQ
(CL:SIN
(CL:LAMBDA (RADIANS) (* kbr: " 7-Apr-86 22:57")
(CTYPECASE RADIANS (COMPLEX (%%SIN-COMPLEX RADIANS))
(FLOAT (%%SIN-FLOAT RADIANS NIL))
(NUMBER (%%SIN-FLOAT (FLOAT RADIANS)
NIL)))))
(%%SIN-FLOAT
(LAMBDA (X COS-FLAG) (* kbr: "23-Jul-86 15:53")
(* * 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 (minus %%PI/2 X)) = (SIN (+ %%PI/2 X)) Case out on sign of X for
improved numerical stability. Avoids unnecessary rounding and promotes
symmetric properties. (COS X) = (COS (minus X)) is guaranteed by this strategy.
*)
(SETQ R (COND
((NOT COS-FLAG)
X)
((> X 0)
(- %%PI/2 X))
(T (+ %%PI/2 X))))
(* * First range reduce to (0 infinity) by
(SIN (minus X)) = (minus (SIN X)) This strategy guarantees
(SIN (minus X)) = (minus (SIN X)) *)
(COND
((< R 0)
(SETQ SIGN -1.0)
(SETQ R (- R)))
(T (SETQ SIGN 1.0)))
(* * Next range reduce to interval (0 %%2PI) by
(SIN X) = (SIN (MOD X %%2PI)) . *)
(SETQ R (REM R %%2PI))
(* * Next range reduce to interval (0 PI) by
(SIN (+ X PI)) = (minus (SIN X)) *)
(COND
((> R %%PI)
(SETQ SIGN (- SIGN))
(SETQ R (- R %%PI))))
(* * 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 (CL:* SIGN R))))
(* * 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 (RADIANS) (* lmm " 1-Jul-86 15:55")
(CTYPECASE RADIANS (FLOAT (%%SIN-FLOAT RADIANS T))
(COMPLEX (%%COS-COMPLEX RADIANS))
(NUMBER (%%SIN-FLOAT (FLOAT RADIANS)
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 *)
(* * %%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 adapted P and Q coefficients of Harris et al TAN 4288
rational approximation to (TAN X) in interval (-PI/4 PI/4) %. The coefficients for %%TAN-PPOLY
and %%TAN-QPOLY have been computed from Harris using extended precision routines and the
relations %%TAN-PPOLY = (REVERSE (for I from 0 as ENTRY in PS collect (/ (CL:* (EXPT (/ 4 PI) (
1+ (CL:* 2 I))) ENTRY) Q0))) and %%TAN-QPOLY = (REVERSE (for I from 0 as ENTRY in QS collect (/
(CL:* (EXPT (/ 4 PI) (CL:* 2 I)) ENTRY) Q0))) *)
(RPAQ %%TAN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 13237 21090)
(%%FLOAT 47141 15825)
(%%FLOAT 15246 8785)
(%%FLOAT 48655 48761)
(%%FLOAT 16256 0))))
(RPAQ %%TAN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 45267 36947)
(%%FLOAT 13848 46875)
(%%FLOAT 47612 53738)
(%%FLOAT 15596 52854)
(%%FLOAT 48882 35303)
(%%FLOAT 16256 0))))
(DEFINEQ
(CL:TAN
(CL:LAMBDA (RADIANS) (* lmm "27-Jun-86 22:41")
(CTYPECASE RADIANS (COMPLEX (%%TAN-COMPLEX RADIANS))
(NUMBER (%%TAN-FLOAT (FLOAT RADIANS))))))
(%%TAN-FLOAT
(LAMBDA (X) (* kbr: "23-Jul-86 17:29")
(* * TAN of a FLOAT X calculated via TAN 4288 rational approximation of Harris
et al. *)
(PROG (R SIGN RECIPFLG R2 ANSWER)
(DECLARE (GLOBALVARS %%TAN-PPOLY %%TAN-QPOLY)
(TYPE FLOATP X R R2 ANSWER))
(SETQ R X)
(* * First range reduce to (0 infinity) by
(TAN (minus X)) = (minus (TAN X)) *)
(COND
((< R 0)
(SETQ SIGN -1.0)
(SETQ R (- R)))
(T (SETQ SIGN 1.0)))
(* * Next range reduce to (0 PI) *)
(SETQ R (REM R %%PI))
(* * Next, range reduce to (-PI/4 PI/4) using
(TAN X) = (TAN (minus X PI)) 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 (CL:* SIGN R))
(COND
(RECIPFLG (SETQ ANSWER (/ ANSWER))))
(RETURN ANSWER)))
(* * 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:* SIGN 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 *)
(* * %%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 (NUMBER) (* kbr: "12-Jul-86 18:05")
(CTYPECASE NUMBER (FLOAT (%%ASIN-FLOAT NUMBER))
(COMPLEX (%%ASIN-COMPLEX NUMBER))
(NUMBER (%%ASIN-FLOAT (FLOAT NUMBER))))))
(%%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 (NUMBER) (* kbr: "12-Jul-86 18:05")
(CTYPECASE NUMBER (FLOAT (%%ASIN-FLOAT NUMBER T))
(COMPLEX (%%ACOS-COMPLEX NUMBER))
(NUMBER (%%ASIN-FLOAT (FLOAT NUMBER)
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 *)
(DECLARE: EVAL@COMPILE DONTCOPY
(DECLARE: EVAL@COMPILE
(RPAQ %%SQRT3 (%%FLOAT 16349 46039))
(RPAQ %%2-SQRT3 (%%FLOAT 16009 12451))
(RPAQ %%INV-2-SQRT3 (%%FLOAT 16494 55788))
(CONSTANTS (%%SQRT3 (%%FLOAT 16349 46039))
(%%2-SQRT3 (%%FLOAT 16009 12451))
(%%INV-2-SQRT3 (%%FLOAT 16494 55788)))
)
(FILESLOAD (LOADCOMP)
LLFLOAT)
)
(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))))))
(%%ATAN-FLOAT1
(CL:LAMBDA (X) (* kbr: "23-Jul-86 18:26")
(COND
((MINUSP X) (* (ATAN (minus X)) =
(minus (ATAN X)) *)
(- (%%ATAN-FLOAT1 (- X))))
((< X %%2-SQRT3) (* (ATAN X) = (minus %%PI/2
(ATAN (/ X))) *)
(%%ATAN-FLOAT2 X))
((< X 1)
(+ %%PI/6 (%%ATAN-FLOAT2 (/ (1- (CL:* X %%SQRT3))
(+ %%SQRT3 X)))))
((< X %%INV-2-SQRT3)
(- %%PI/3 (%%ATAN-FLOAT2 (/ (- %%SQRT3 X)
(1+ (CL:* X %%SQRT3))))))
(T (- %%PI/2 (%%ATAN-FLOAT2 (/ X)))))))
(%%ATAN-FLOAT2
(CL:LAMBDA (X) (* kbr: "23-Jul-86 18:26")
(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)
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
(%%ATAN-FLOAT1 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))
(%%ATAN-FLOAT1 (/ Y X)))
((PLUSP Y)
(- %%PI (%%ATAN-FLOAT1 (/ (- Y)
X))))
((PLUSP X)
(- (%%ATAN-FLOAT1 (/ Y (- X)))))
(T (- (%%ATAN-FLOAT1 (/ 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 (RADIANS) (* kbr: "12-Jul-86 18:05")
(* Return cos (Theta) + i sin
(Theta)%. *)
(CL:IF (COMPLEXP RADIANS)
(CL:ERROR "Argument to CIS is complex: ~S" RADIANS)
(COMPLEX (CL:COS RADIANS)
(CL:SIN RADIANS)))))
)
(* SINH COSH TANH *)
(DEFINEQ
(SINH
(CL:LAMBDA (NUMBER) (* kbr: "12-Jul-86 18:05")
(* 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 NUMBER)))
(/ (- Z (/ Z))
2))))
(COSH
(CL:LAMBDA (NUMBER) (* kbr: "12-Jul-86 18:05")
(LET ((Z (EXP NUMBER)))
(/ (+ Z (/ Z))
2))))
(TANH
(CL:LAMBDA (NUMBER) (* kbr: "12-Jul-86 18:05")
(* Different form than in the manual.
Does much better. *)
(LET* ((Z (EXP (CL:* 2 NUMBER)))
(Y (/ Z)))
(- (/ (1+ Y))
(/ (1+ Z))))))
)
(* ASINH ACOSH ATANH *)
(DEFINEQ
(ASINH
(CL:LAMBDA (NUMBER) (* kbr: "12-Jul-86 18:05")
(CL:LOG (+ NUMBER (CL:SQRT (+ (CL:* NUMBER NUMBER)
1))))))
(ACOSH
(CL:LAMBDA (NUMBER) (* kbr: "12-Jul-86 18:05")
(CL:LOG (+ NUMBER (CL:SQRT (- (CL:* NUMBER NUMBER)
1))))))
(ATANH
(CL:LAMBDA (NUMBER) (* kbr: "12-Jul-86 18:05")
(CL:IF (%%ATANH-DOMAIN-CHECK NUMBER)
(CL:* .5 (CL:LOG (/ (1+ NUMBER)
(- 1 NUMBER))))
(CL:ERROR "~S argument out of range." NUMBER))))
(%%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 %%ATAN-FLOAT2 %%ATAN-FLOAT1 CL:ATAN %%ACOS-COMPLEX ACOS
%%ASIN-COMPLEX ASIN %%TAN-COMPLEX CL:TAN %%COS-COMPLEX CL:COS %%SIN-COMPLEX CL:SIN
%%SQRT-PERFECT CL:SQRT %%LOG-COMPLEX CL:LOG %%EXPT-COMPLEX-POWER %%EXPT-COMPLEX
%%EXPT-INTEGER CL:EXPT %%EXP-COMPLEX %%EXP-FLOAT EXP)
)
(PRETTYCOMPRINT CMLFLOATCOMS)
(RPAQQ CMLFLOATCOMS
((* * CMLFLOAT -- Covering sections 12.5-12.5.3 irrational, transcendental, exponential,
logarithmic, trigonometric, and hyperbolic functions. Section 12.10, implementation
parameters. -- By Kelly Roach. *)
(COMS (* Section 12.10, implementation parameters. The constants in this COMS are exported to
the user. *)
(* %%FLOAT allows us to recreate FLOATPs in a way that is independent of the ordinairy
reading and printing FLOATPs to files which involves loss of the last couple bits of
accuracy due to rounding effects. *)
(FNS %%FLOAT)
(CONSTANTS (MOST-POSITIVE-FIXNUM 65535)
(MOST-NEGATIVE-FIXNUM -65536)
(MOST-POSITIVE-SINGLE-FLOAT (%%FLOAT 32639 65535))
(LEAST-POSITIVE-SINGLE-FLOAT (%%FLOAT 128 0))
(LEAST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 32896 0))
(MOST-NEGATIVE-SINGLE-FLOAT (%%FLOAT 65407 65535))
(MOST-POSITIVE-SHORT-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-SHORT-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-SHORT-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-SHORT-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(MOST-POSITIVE-DOUBLE-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-DOUBLE-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-DOUBLE-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-DOUBLE-FLOAT MOST-NEGATIVE-SINGLE-FLOAT)
(MOST-POSITIVE-LONG-FLOAT MOST-POSITIVE-SINGLE-FLOAT)
(LEAST-POSITIVE-LONG-FLOAT LEAST-POSITIVE-SINGLE-FLOAT)
(LEAST-NEGATIVE-LONG-FLOAT LEAST-NEGATIVE-SINGLE-FLOAT)
(MOST-NEGATIVE-LONG-FLOAT MOST-NEGATIVE-SINGLE-FLOAT))
(* EPSILON is the smallest positive floating point number such that
(NOT (= (FLOAT 1 EPSILON)
(+ (FLOAT 1 EPSILON)
EPSILON)))
*)
(CONSTANTS (SINGLE-FLOAT-EPSILON (%%FLOAT 13312 0))
(SHORT-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(DOUBLE-FLOAT-EPSILON SINGLE-FLOAT-EPSILON)
(LONG-FLOAT-EPSILON SINGLE-FLOAT-EPSILON))
(* NEGATIVE-EPSILON is the smallest negative floating point number such that
(NOT (= (FLOAT 1 NEGATIVE-EPSILON)
(- (FLOAT 1 NEGATIVE-EPSILON)
NEGATIVE-EPSILON)))
*)
(CONSTANTS (SINGLE-FLOAT-NEGATIVE-EPSILON (%%FLOAT 13184 0))
(SHORT-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(DOUBLE-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)
(LONG-FLOAT-NEGATIVE-EPSILON SINGLE-FLOAT-NEGATIVE-EPSILON)))
(COMS (* Miscellaneous. *)
(DECLARE: EVAL@COMPILE DONTCOPY (FILES (LOADCOMP)
LLFLOAT))
(CONSTANTS (PI (%%FLOAT 16457 4059)))
(DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%E (%%FLOAT 16429 63572))
(%%2PI (%%FLOAT 16585 4059))
(%%PI (%%FLOAT 16457 4059))
(%%2PI/3 (%%FLOAT 16390 2706))
(%%PI/2 (%%FLOAT 16329 4059))
(%%-PI/2 (%%FLOAT 49097 4059))
(%%PI/3 (%%FLOAT 16262 2706))
(%%PI/4 (%%FLOAT 16201 4059))
(%%-PI/4 (%%FLOAT 48969 4059))
(%%PI/6 (%%FLOAT 16134 2706))
(%%2/PI (%%FLOAT 16162 63875))))
(FNS %%MAKE-ARRAY))
(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)) . *)
(VARS (%%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)) . *)
(VARS (%%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))
(COMS (* SQRT *)
(FNS CL:SQRT %%SQRT-FLOAT %%SQRT-COMPLEX))
(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 (CONSTANTS (%%SIN-EPSILON (%%FLOAT 14720 0))))
(* * %%SIN-PPOLY and %%SIN-QPOLY contain adapted P and Q coefficients of Harris
et al SIN 3374 rational approximation to (SIN X)
in interval (0 (/ PI 2))
%. The coefficients for %%SIN-PPOLY and %%SIN-QPOLY have been computed from
Harris using extended precision routines and the relations %%SIN-PPOLY =
(REVERSE (for I from 0 as ENTRY in PS collect
(/ (CL:* (EXPT (/ 2 PI)
(1+ (CL:* 2 I)))
ENTRY)
Q0)))
and %%SIN-QPOLY = (REVERSE (for I from 0 as ENTRY in QS collect
(/ (CL:* (EXPT (/ 2 PI)
(CL:* 2 I))
ENTRY)
Q0)))
*)
(VARS (%%SIN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 45236 25611)
(%%FLOAT 13589 26148)
(%%FLOAT 47286 34797)
(%%FLOAT 15295 3306)
(%%FLOAT 48666 34805)
(%%FLOAT 16256 0))))
(%%SIN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 11384 52865)
(%%FLOAT 12553 9550)
(%%FLOAT 13604 38385)
(%%FLOAT 14593 18841)
(%%FLOAT 15489 5549)
(%%FLOAT 16256 0))))))
(FNS CL:SIN %%SIN-FLOAT %%SIN-COMPLEX))
(COMS (* COS *)
(FNS CL:COS %%COS-COMPLEX))
(COMS (* TAN *)
(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 (CONSTANTS (%%TAN-EPSILON (%%FLOAT 14720 0))))
(* * %%TAN-PPOLY and %%TAN-QPOLY contain adapted P and Q coefficients of Harris
et al TAN 4288 rational approximation to (TAN X)
in interval (-PI/4 PI/4)
%. The coefficients for %%TAN-PPOLY and %%TAN-QPOLY have been computed from
Harris using extended precision routines and the relations %%TAN-PPOLY =
(REVERSE (for I from 0 as ENTRY in PS collect
(/ (CL:* (EXPT (/ 4 PI)
(1+ (CL:* 2 I)))
ENTRY)
Q0)))
and %%TAN-QPOLY = (REVERSE (for I from 0 as ENTRY in QS collect
(/ (CL:* (EXPT (/ 4 PI)
(CL:* 2 I))
ENTRY)
Q0)))
*)
(VARS (%%TAN-PPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 13237 21090)
(%%FLOAT 47141 15825)
(%%FLOAT 15246 8785)
(%%FLOAT 48655 48761)
(%%FLOAT 16256 0))))
(%%TAN-QPOLY (%%MAKE-ARRAY (LIST (%%FLOAT 45267 36947)
(%%FLOAT 13848 46875)
(%%FLOAT 47612 53738)
(%%FLOAT 15596 52854)
(%%FLOAT 48882 35303)
(%%FLOAT 16256 0))))))
(FNS CL:TAN %%TAN-FLOAT %%TAN-COMPLEX))
(COMS (* ASIN *)
(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 (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)) . *)
(VARS (%%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 *)
(DECLARE: EVAL@COMPILE DONTCOPY (CONSTANTS (%%SQRT3 (%%FLOAT 16349 46039))
(%%2-SQRT3 (%%FLOAT 16009 12451))
(%%INV-2-SQRT3 (%%FLOAT 16494 55788)))
(FILES (LOADCOMP)
LLFLOAT))
(FNS CL:ATAN %%ATAN-FLOAT1 %%ATAN-FLOAT2 %%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))
(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 %%ATAN-FLOAT2 %%ATAN-FLOAT1 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)))))
(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 %%ATAN-FLOAT2 %%ATAN-FLOAT1 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))
(DECLARE: DONTCOPY
(FILEMAP (NIL (18287 18447 (%%FLOAT 18297 . 18445)) (23075 23929 (%%MAKE-ARRAY 23085 . 23927)) (25114
27838 (EXP 25124 . 25631) (%%EXP-FLOAT 25633 . 27490) (%%EXP-COMPLEX 27492 . 27836)) (27858 31958 (
CL:EXPT 27868 . 29594) (%%EXPT-INTEGER 29596 . 30770) (%%EXPT-FLOAT 30772 . 31082) (%%EXPT-COMPLEX
31084 . 31509) (%%EXPT-COMPLEX-POWER 31511 . 31956)) (32978 35512 (CL:LOG 32988 . 33542) (%%LOG-FLOAT
33544 . 35309) (%%LOG-COMPLEX 35311 . 35510)) (35532 38884 (CL:SQRT 35542 . 36078) (%%SQRT-FLOAT 36080
. 37612) (%%SQRT-COMPLEX 37614 . 38882)) (40773 44100 (CL:SIN 40783 . 41107) (%%SIN-FLOAT 41109 .
43738) (%%SIN-COMPLEX 43740 . 44098)) (44119 44820 (CL:COS 44129 . 44448) (%%COS-COMPLEX 44450 . 44818
)) (46651 49743 (CL:TAN 46661 . 46895) (%%TAN-FLOAT 46897 . 49215) (%%TAN-COMPLEX 49217 . 49741)) (
51408 54434 (ASIN 51418 . 51700) (%%ASIN-FLOAT 51702 . 54143) (%%ASIN-COMPLEX 54145 . 54432)) (54454
55093 (ACOS 54464 . 54781) (%%ACOS-COMPLEX 54783 . 55091)) (55478 58766 (CL:ATAN 55488 . 55798) (
%%ATAN-FLOAT1 55800 . 56718) (%%ATAN-FLOAT2 56720 . 57087) (%%ATAN-DOMAIN-CHECK 57089 . 57591) (
%%ATAN-FLOAT 57593 . 58382) (%%ATAN-COMPLEX 58384 . 58764)) (58785 59271 (CIS 58795 . 59269)) (59301
60449 (SINH 59311 . 59816) (COSH 59818 . 60016) (TANH 60018 . 60447)) (60482 61774 (ASINH 60492 .
60713) (ACOSH 60715 . 60936) (ATANH 60938 . 61263) (%%ATANH-DOMAIN-CHECK 61265 . 61772)))))
STOP