(DEFINE-FILE-INFO READTABLE "XCL" PACKAGE "CL")
(il:filecreated "10-Nov-88 14:10:25" il:{qv}<idl>next>dorado-cmlfloat.\;2 23210
il:|changes| il:|to:| (il:functions %exp-float*)
il:|previous| il:|date:| "24-Aug-88 16:48:02" il:{qv}<idl>next>dorado-cmlfloat.\;1)
; Copyright (c) 1988 by Xerox Corporation. All rights reserved.
(il:prettycomprint il:dorado-cmlfloatcoms)
(il:rpaqq il:dorado-cmlfloatcoms ((il:declare\: il:dontcopy il:doeval@compile (il:* il:\; "To get constants from llfloat ") (il:files il:cmlarray-support) (il:files (il:loadcomp) il:llfloat)) (il:coms (il:* il:|;;| "Utilities") (il:functions copy-to-boxed-vector %poly-eval) (il:declare\: il:dontcopy il:doeval@compile (il:functions %ftruncate %float-unbox %umake-float)) (il:* il:|;;| "Internal constants") (il:declare\: il:dontcopy il:doeval@compile (il:variables %2pi %pi/2 %-pi/2 %pi/3 %pi/4 %-pi/4 %pi/6))) (il:coms (il:* il:|;;| "Exp (e to the power x)") (il:declare\: il:dontcopy il:doeval@compile (il:variables %log-base2-e)) (il:variables %boxed-exp-table %boxed-exp-poly) (il:functions %exp-float*)) (il:coms (il:* il:|;;| "Expt (x to the power y)") (il:functions %expt-float-integer*)) (il:coms (il:* il:|;;| "Log (log base e)") (il:declare\: il:dontcopy il:doeval@compile (il:variables %log2 %sqrt2)) (il:variables %boxed-log-ppoly %boxed-log-qpoly) (il:functions %log-float*)) (il:coms (il:* il:|;;| "Sqrt") (il:functions %sqrt-float*)) (il:coms (il:* il:|;;| "Sin and Cos") (il:declare\: il:dontcopy il:doeval@compile (il:variables %sin-epsilon)) (il:variables %boxed-sin-ppoly %boxed-sin-qpoly) (il:functions %sin-float*)) (il:coms (il:* il:|;;| "Tan") (il:coms (il:declare\: il:dontcopy il:doeval@compile (il:variables %tan-epsilon)) (il:variables %boxed-tan-ppoly %boxed-tan-qpoly) (il:functions %tan-float*))) (il:coms (il:* il:|;;| "Asin and Acos") (il:coms (il:declare\: il:dontcopy il:doeval@compile (il:variables %asin-epsilon)) (il:variables %boxed-asin-ppoly %boxed-asin-qpoly) (il:functions %asin-float*))) (il:coms (il:* il:|;;| "Atan ") (il:declare\: il:dontcopy il:doeval@compile (il:variables %sqrt3 %2-sqrt3 %inv-2-sqrt3)) (il:functions %atan-float*)) (il:* il:|;;| "Could do sinh and friends if anyone cares. ") (il:coms (il:* il:|;;| "Misc, cleanup from cmlarith") (il:functions xcl::struncate* mod* rem*) (il:* il:|;;| "Cleanup of random") (il:functions random*)) (il:coms (il:* il:|;;| "Installation") (il:variables *fns-to-replace*) (il:functions install-boxed-fns deinstall-boxed-fns) (il:declare\: il:docopy il:donteval@load il:donteval@compile (il:p (install-boxed-fns)))) (il:declare\: il:dontcopy il:doeval@compile (il:localvars . t)) (xcl:file-environments "DORADO-CMLFLOAT")))
(il:declare\: il:dontcopy il:doeval@compile
(il:filesload il:cmlarray-support)
(il:filesload (il:loadcomp) il:llfloat)
)
(il:* il:|;;| "Utilities")
(defun copy-to-boxed-vector (float-vector) (let* ((length (length float-vector)) (vector (make-array length))) (dotimes (i length) (setf (aref vector i) (aref float-vector i))) vector))
(defun %poly-eval (x coefficients degree) (il:* il:|;;| "X is a float, COEFFICIENTS is a (pointer) vector of COEFFICIENTS, DEGREE is an integer") (do* ((base (xcl:record-fetch oned-array base coefficients)) (result (il:\\getbaseptr base 0)) (limit (il:llsh degree 1)) (i 2 (il:iplus i 2))) ((il:igreaterp i limit) result) (setq result (il:fplus (il:ftimes result x) (il:\\getbaseptr base i)))))
(il:declare\: il:dontcopy il:doeval@compile
(defmacro %ftruncate (int rem float &optional divisor) (il:* il:|;;| "As in truncate. Assumes FLOAT and DIVISOR are float's. INT and REM are symbols") (if divisor (xcl:once-only (float divisor) (il:bquote (progn (setq (il:\\\, int) (truncate (il:\\\, float) (il:\\\, divisor))) (setq (il:\\\, rem) (il:fdifference (il:\\\, float) (il:ftimes (il:\\\, divisor) (float (il:\\\, int))))) nil))) (xcl:once-only (float) (il:bquote (progn (setq (il:\\\, int) (truncate (il:\\\, float))) (setq (il:\\\, rem) (il:fdifference (il:\\\, float) (float (il:\\\, int)))) nil)))))
(defmacro %float-unbox (float sign exp hi lo &optional dontshift) (il:* il:|;;| "If dontshift is T -- the floatp fields are simply unpacked (with the hiddenbit restored -- and exp set to 1 for denormalized numbers). If dontshift is NIL -- exp, hi and lo are fiddled so the high bit of hi is on.") (il:bquote (let ((flonum (float (il:\\\, float)))) (setq (il:\\\, sign) (il:|fetch| (il:floatp il:signbit) il:|of| flonum)) (setq (il:\\\, exp) (il:|fetch| (il:floatp il:exponent) il:|of| flonum)) (setq (il:\\\, hi) (il:|fetch| (il:floatp il:hifraction) il:|of| flonum)) (setq (il:\\\, lo) (il:|fetch| (il:floatp il:lofraction) il:|of| flonum)) (if (eq (il:\\\, exp) il:\\max.exponent) (il:* il:\; "might want to check for NaN's here if EXP = \\MAX.EXPONENT") (error "Not a number: ~s" flonum)) (if (eq 0 (il:\\\, exp)) (when (not (and (eq 0 (il:\\\, hi)) (eq 0 (il:\\\, lo)))) (il:* il:\; "Denormalized number") (setq (il:\\\, exp) 1) (il:\\\,@ (if (null dontshift) (il:bquote ((loop (if (not (eq 0 (logand (il:\\\, hi) il:\\hiddenbit))) (return nil)) (il:.llsh1. (il:\\\, hi) (il:\\\, lo)) (setq (il:\\\, exp) (1- (il:\\\, exp))))))))) (il:* il:\; " Restore the hidden bit") (setq (il:\\\, hi) (+ (il:\\\, hi) il:\\hiddenbit))) (il:\\\,@ (if (null dontshift) (il:bquote ((il:.llsh8. (il:\\\, hi) (il:\\\, lo)))))) nil)))
(defmacro %umake-float (sign exp hi low) (il:* il:|;;| "as in \\makefloat -- but produces an unboxed number") (il:bquote (il:\\floatbox ((il:openlambda (sign exp hi lo) (il:.lrsh8. hi lo) (setq hi (+ (ash exp 7) (logand 127 hi))) (if (eq sign 1) (setq hi (logior il:\\signbit hi))) (il:\\vag2 hi lo)) (il:\\\, sign) (il:\\\, exp) (il:\\\, hi) (il:\\\, low)))))
)
(il:* il:|;;| "Internal constants")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %2pi (%float 16585 4059))
(defconstant %pi/2 (%float 16329 4059))
(defconstant %-pi/2 (%float 49097 4059))
(defconstant %pi/3 (%float 16262 2706))
(defconstant %pi/4 (%float 16201 4059))
(defconstant %-pi/4 (%float 48969 4059))
(defconstant %pi/6 (%float 16134 2706))
)
(il:* il:|;;| "Exp (e to the power x)")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %log-base2-e (%float 16312 43579))
)
(defvar %boxed-exp-table (copy-to-boxed-vector %exp-table))
(defvar %boxed-exp-poly (copy-to-boxed-vector %exp-poly))
(defun %exp-float* (x) (il:* il:|;;| "(CL:EXP X) for float X calculated via EXPB 1103 rational approximation of Hart et al. ") (let ((fx (float x)) r m n answer recipflg) (il:* il:|;;| "First, arrange X to be in interval (0 infinity) via identity (CL:EXP (minus X)) = (/ 1.0 (CL:EXP X))") (when (il:flessp fx 0.0) (setq fx (il:fminus fx)) (setq recipflg t)) (il:* il:|;;| "Next, the problem of (CL:EXP X) is converted into a problem (EXPT 2 Y) where Y = (* %LOG-BASE2-E X). ") (il:* il:|;;| "Then range reduction is accomplished via (EXPT 2 Y) = (* (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). ") (il:* il:|;;| "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 Hart et al. ") (%ftruncate m r (il:ftimes %log-base2-e fx)) (%ftruncate n r r 0.125) (setq fx (il:ftimes (aref %boxed-exp-table n) (il:fquotient (%poly-eval r %boxed-exp-poly 5) (%poly-eval (il:fminus r) %boxed-exp-poly 5)))) (cond (recipflg (setq answer (setq fx (il:fquotient 1.0 fx))) (scale-float answer (- m) answer)) (t (setq answer fx) (scale-float answer m answer)))))
(il:* il:|;;| "Expt (x to the power y)")
(defun %expt-float-integer* (base power) (il:* il:|;;| "(EXPT BASE POWER) where BASE is a float and POWER is an integer. ") (if (il:igreaterp 0 power) (il:fquotient 1.0 (%expt-float-integer* base (il:iminus power))) (il:* il:|;;| "float to positive integer power") (let ((fbase (float base)) (total 1.0)) (loop (if (oddp power) (setq total (il:ftimes fbase total))) (setq power (il:lrsh power 1)) (if (eq 0 power) (return total)) (setq fbase (il:ftimes fbase fbase))))))
(il:* il:|;;| "Log (log base e)")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %log2 (%float 16177 29208))
(defconstant %sqrt2 (%float 16309 1267))
)
(defvar %boxed-log-ppoly (copy-to-boxed-vector %log-ppoly))
(defvar %boxed-log-qpoly (copy-to-boxed-vector %log-qpoly))
(defun %log-float* (x) (il:* il:|;;| "(LOG X) for positive float X. ") (if (<= (setq x (float x)) 0.0) (error "Log of zero: ~s" x)) (il:* il:|;;| "Range reduce to an R in interval ((SQRT 0.5) (SQRT 2)) via identity (LOG X) = (+ (LOG R) (* %LOG-2 EXP)) for a suitable integer EXP. exp is found from the exponent field of the iee floating point number representation of x.") (let (r exp) (let (sign hi lo) (%float-unbox x sign exp hi lo) (setq exp (il:idifference exp il:\\exponent.bias)) (setq r (%umake-float sign il:\\exponent.bias hi lo))) (when (il:fgreaterp r %sqrt2) (setq exp (il:iplus 1 exp)) (setq r (il:fquotient r 2.0))) (il:* il:|;;| "(LOG R) is calculated by rational approximation LOGE 2707 of Hart et al.") (let* ((z (il:fquotient (il:fdifference r 1.0) (il:fplus 1.0 r))) (z2 (il:ftimes z z))) (il:fplus (il:ftimes z (il:fquotient (%poly-eval z2 %boxed-log-ppoly 4) (%poly-eval z2 %boxed-log-qpoly 4))) (il:ftimes %log2 exp)))))
(il:* il:|;;| "Sqrt")
(defun %sqrt-float* (x) (il:* il:|;;| "(SQRT X) for nonnegative float X") (if (<= x 0.0) 0.0 (let* ((fx (float x)) (il:* il:|;;| "First guess") (v (let (sign exp hi lo) (%float-unbox x sign exp hi lo) (il:* il:|;;| "First guess") (%umake-float 0 (+ (il:lrsh exp 1) (1+ (il:lrsh il:\\exponent.bias 1))) hi lo)))) (il:* il:|;;| "Four step newton-raphson") (dotimes (i 4 v) (setq v (il:ftimes 0.5 (il:fplus v (il:fquotient fx v))))))))
(il:* il:|;;| "Sin and Cos")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %sin-epsilon (il:* il:|;;| "%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 (* 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.") (%float 14720 0))
)
(defvar %boxed-sin-ppoly (copy-to-boxed-vector %sin-ppoly))
(defvar %boxed-sin-qpoly (copy-to-boxed-vector %sin-qpoly))
(defun %sin-float* (x cos-flg) (il:* il:|;;| "SIN of a FLOAT X calculated via SIN 3374 rational approximation of Hart et al. ") (let ((theta (float x)) (sign 1.0)) (il:* il:|;;| "If this function is called by COS then use (COS X) = (SIN (-- %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 (-- X)) is guaranteed by this strategy.") (if cos-flg (if (il:fgreaterp theta 0.0) (setq theta (il:fdifference %pi/2 theta)) (setq theta (il:fplus %pi/2 theta)))) (il:* il:|;;| "First range reduce to (0 infinity) by (SIN (minus X)) = (minus (SIN X)) This strategy guarantees (SIN (minus X)) = (minus (SIN X))") (when (il:flessp theta 0.0) (setq sign -1.0) (setq theta (il:fminus theta))) (il:* il:|;;| "Next range reduce to interval (0 %2PI) by (SIN X) = (SIN (MOD X %2PI)) ") (if (>= theta %2pi) (setq theta (il:fdifference theta (il:ftimes %2pi (float (il:fix (il:fquotient theta %2pi))))))) (il:* il:|;;| "Next range reduce to interval (0 PI) by (SIN (+ X PI)) = (minus (SIN X)) ") (when (il:fgreaterp theta pi) (setq theta (il:fdifference theta pi)) (setq sign (il:fminus sign))) (il:* il:|;;| "Next range reduce to interval (0 %PI/2) by (SIN (+ X %PI/2)) = (SIN (minus %PI/2 X))") (if (il:fgreaterp theta %pi/2) (setq theta (il:fdifference pi theta))) (if (il:flessp theta %sin-epsilon) (il:* il:|;;| "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") (setq theta (il:ftimes sign theta)) (il:* il:|;;| "Now use SIN 3374 rational approximation of Harris et al. which works on interval (0 %PI/2) ") (let ((r2 (il:ftimes theta theta))) (il:ftimes sign theta (il:fquotient (%poly-eval r2 %boxed-sin-ppoly 5) (%poly-eval r2 %boxed-sin-qpoly 5)))))))
(il:* il:|;;| "Tan")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %tan-epsilon (il:* il:|;;| "%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 (* 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.") (%float 14720 0))
)
(defvar %boxed-tan-ppoly (copy-to-boxed-vector %tan-ppoly))
(defvar %boxed-tan-qpoly (copy-to-boxed-vector %tan-qpoly))
(defun %tan-float* (x) (il:* il:|;;| "TAN of a FLOAT X calculated via TAN 4288 rational approximation of Hart et al.") (let ((fx (float x)) (sign 1.0) recipflg) (il:* il:|;;| "First range reduce to (0 infinity) by (TAN (minus X)) = (minus (TAN X))") (when (il:flessp fx 0.0) (setq sign -1.0) (setq fx (il:fminus fx))) (il:* il:|;;| "Next range reduce to (0 PI)") (if (>= fx pi) (setq fx (il:fdifference fx (il:ftimes pi (float (il:fix (il:fquotient fx pi))))))) (il:* il:|;;| "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) ") (if (il:fgreaterp fx %pi/2) (progn (setq fx (il:fdifference fx pi)) (when (il:flessp fx %-pi/4) (setq recipflg t) (setq fx (il:fdifference %-pi/2 fx)))) (when (il:fgreaterp fx %pi/4) (setq recipflg t) (setq fx (il:fdifference %pi/2 fx)))) (if (il:flessp (il:fabs fx) %tan-epsilon) (il:* il:|;;| "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.") (progn (setq fx (il:ftimes sign fx)) (if recipflg (setq fx (il:fquotient 1.0 fx)) fx)) (il:* il:|;;| "Now use TAN 4288 rational approximation of Hart et al. which works on interval (0 %PI/4)") (let ((r2 (il:ftimes fx fx))) (setq fx (il:ftimes sign fx (il:fquotient (%poly-eval r2 %boxed-tan-ppoly 4) (%poly-eval r2 %boxed-tan-qpoly 5)))) (if recipflg (il:fquotient 1.0 fx) fx)))))
(il:* il:|;;| "Asin and Acos")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %asin-epsilon (il:* il:|;;| "%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 (* 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.") (%float 14720 0))
)
(defvar %boxed-asin-ppoly (copy-to-boxed-vector %asin-ppoly))
(defvar %boxed-asin-qpoly (copy-to-boxed-vector %asin-qpoly))
(defun %asin-float* (x acos-flg) (il:* il:|;;| "(ASIN X) for float X calculated via ARCSN 4671 rational approximation of Hart et al.") (if (or (il:flessp x -1.0) (il:fgreaterp x 1.0)) (error "Arg not in range: ~s" x)) (let ((fx (float x)) negative reduced) (il:* il:|;;| "Range reduce to (0 1) via identity (ASIN (minus X)) = (minus (ASIN X)) ") (when (il:flessp fx 0.0) (setq negative t) (setq fx (il:fminus fx))) (il:* il:|;;| "Range reduce to (0 0.5) via identity (ASIN X) = (minus %PI/2 (* 2.0 (ASIN (SQRT (* 0.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") (when (il:fgreaterp fx 0.5) (setq reduced t) (setq fx (sqrt (il:ftimes 0.5 (il:fdifference 1.0 fx))))) (il:* il:|;;| "R is now in range (0 0.5) Use ARCSN 4671 rational approximation to calculate (ASIN R) ") (if (il:fgreaterp fx %asin-epsilon) (il:* il:|;;| "If R is in the interval (0 %SIN-EPSILON) then (ASIN R) = R to the precision that we can offer. ") (let ((r2 (il:ftimes fx fx))) (setq fx (il:ftimes fx (il:quotient (%poly-eval r2 %boxed-asin-ppoly 6) (%poly-eval r2 %boxed-asin-qpoly 6)))))) (if reduced (setq fx (il:fdifference %pi/2 (il:ftimes 2.0 fx)))) (if negative (setq fx (il:fminus fx))) (il:* il:|;;| "In case we want (ACOS X) then use identity (ACOS X) = (minus %PI/2 (ASIN X))") (if acos-flg (setq fx (il:fdifference %pi/2 fx))) fx))
(il:* il:|;;| "Atan ")
(il:declare\: il:dontcopy il:doeval@compile
(defconstant %sqrt3 (%float 16349 46039))
(defconstant %2-sqrt3 (%float 16009 12451))
(defconstant %inv-2-sqrt3 (%float 16494 55788))
)
(defun %atan-float* (y &optional x) (let ((fy (float y)) fx farg) (il:* il:|;;| "Compute farg") (cond ((null x) (if (= y 0.0) (return-from %atan-float* 0.0) (setq farg fy))) (t (setq fx (float x)) (cond ((= x 0.0) (if (= y 0.0) (error "Both args to atan are 0.0") (return-from %atan-float* (if (il:fgreaterp y 0.0) %pi/2 (il:fminus %pi/2))))) ((= y 0.0) (return-from %atan-float* (if (il:fgreaterp x 0.0) 0.0 pi))) ((il:fgreaterp y 0.0) (if (il:fgreaterp x 0.0) (setq farg (il:fquotient fy fx)) (setq farg (il:fquotient (il:fminus fy) fx)))) ((il:fgreaterp x 0.0) (setq farg (il:fquotient fy (il:fminus fx)))) (t (setq farg (il:fquotient fy fx)))))) (il:* il:|;;| "Compute result") (let ((constant 0.0) (constant-flag t) negate-flag add-flag) (il:* il:|;;| "(ATAN (minus X)) = (minus (ATAN X)) ") (when (il:flessp farg 0.0) (setq negate-flag t) (setq farg (il:fminus farg))) (il:* il:|;;| "Range reduce to (0, 2-sqrt(3))") (cond ((>= farg %inv-2-sqrt3) (il:* il:|;;| "(ATAN X) = (minus %PI/2 (ATAN (/ X)))") (setq constant %pi/2) (setq farg (il:fquotient 1.0 farg))) ((>= farg 1.0) (setq constant %pi/3) (setq farg (il:fquotient (il:fdifference %sqrt3 farg) (il:fplus 1.0 (il:ftimes farg %sqrt3))))) ((>= farg %2-sqrt3) (setq add-flag t) (setq constant %pi/6) (setq farg (il:fquotient (il:fdifference (il:ftimes farg %sqrt3) 1.0) (il:fplus %sqrt3 farg)))) (t (setq constant-flag nil))) (il:* il:|;;| "Power series expansion cons'ed up on the fly") (let ((sqr (il:fminus (il:ftimes farg farg))) (int 1.0) (pow farg) (old 0.0)) (loop (if (= farg old) (return nil)) (setq int (il:fplus int 2.0)) (setq pow (il:ftimes pow sqr)) (setq old farg) (setq farg (il:fplus farg (il:fquotient pow int))))) (if constant-flag (if add-flag (setq farg (il:fplus constant farg)) (setq farg (il:fdifference constant farg)))) (if negate-flag (setq farg (il:fminus farg)))) (il:* il:|;;| "Fix up") (if x (cond ((il:fgreaterp fy 0.0) (if (il:flessp fx 0.0) (setq farg (il:fdifference pi farg)))) ((il:fgreaterp fx 0.0) (setq farg (il:fminus farg))) (t (setq farg (il:fdifference farg pi))))) (il:* il:|;;| "Box and return") farg))
(il:* il:|;;| "Could do sinh and friends if anyone cares. ")
(il:* il:|;;| "Misc, cleanup from cmlarith")
(defun xcl::struncate* (number &optional divisor) (il:* il:|;;| "Returns number (or number/divisor) as an integer, rounded toward 0.0 ") (if (null divisor) (typecase number (float (il:* il:|;;| "Could be (IL:FIX NUMBER), but this is slightly faster") (il:\\fixp.from.floatp number)) (ratio (il:iquotient (ratio-numerator number) (ratio-denominator number))) (integer number) (otherwise (%not-noncomplex-number-error number))) (typecase divisor (integer (il:iquotient number divisor)) (float (il:fix (il:fquotient number divisor))) (ratio (xcl::struncate* (/ number divisor))) (otherwise (%not-noncomplex-number-error divisor)))))
(defun mod* (number divisor) (il:* il:|;;| "Returns second result of FLOOR.") (if (or (floatp number) (floatp divisor)) (let ((fx (float number)) (fy (float divisor)) rem) (setq rem (il:fdifference fx (il:ftimes (float (il:fix (il:fquotient fx fy))) fy))) (if (= rem 0.0) 0.0 (if (if (il:fgreaterp 0.0 fy) (il:fgreaterp fx 0.0) (il:fgreaterp 0.0 fx)) (setq rem (il:fplus rem fy)))) rem) (let ((rem (rem number divisor))) (if (and (not (zerop rem)) (if (minusp divisor) (plusp number) (minusp number))) (+ rem divisor) rem))))
(defun rem* (number divisor) (il:* il:|;;| "Returns the second value of truncate") (cond ((and (integerp number) (integerp divisor)) (il:iremainder number divisor)) ((or (floatp number) (floatp divisor)) (let ((fx (float number)) (fy (float divisor))) (il:fdifference fx (il:ftimes (float (il:fix (il:fquotient fx fy))) fy)))) (t (- number (* divisor (xcl::struncate number divisor))))))
(il:* il:|;;| "Cleanup of random")
(defun random* (number &optional (state *random-state*)) (if (not (> number 0)) (error "Not a positive number: ~s" number)) (let ((rv (%random state))) (typecase number (fixnum (if (eq number most-positive-fixnum) rv (il:iremainder rv number))) (float (il:ftimes number (il:fquotient (float rv) (float (1+ most-positive-fixnum))))) (integer (do ((tot rv (+ (ash tot 16) (%random state))) (end (ash number -16) (ash end -16))) ((eq 0 end) (rem tot number)))) (t (error (quote xcl:type-mismatch) :expected-type (quote (or integer float)) :name number :value number :message "an integer or a float")))))
(il:* il:|;;| "Installation")
(defvar *fns-to-replace* (quote (%exp-float %expt-float-integer %log-float %sqrt-float %sin-float %tan-float %asin-float %atan-float xcl::struncate mod rem random)))
(defun install-boxed-fns (recover-p) (let ((targets *fns-to-replace*) (sources (mapcar (function (lambda (sym) (intern (concatenate (quote string) (string sym) "*") (symbol-package sym)))) *fns-to-replace*)) (saves (mapcar (function (lambda (sym) (intern (concatenate (quote string) (string sym) "-save") (symbol-package sym)))) *fns-to-replace*))) (if recover-p (do* ((target-tail targets (cdr target-tail)) (target (car target-tail) (car target-tail)) (save-tail saves (cdr save-tail)) (save (car save-tail) (car save-tail))) ((null target-tail)) (il:movd save target)) (do* ((target-tail targets (cdr target-tail)) (target (car target-tail) (car target-tail)) (source-tail sources (cdr source-tail)) (source (car source-tail) (car source-tail)) (save-tail saves (cdr save-tail)) (save (car save-tail) (car save-tail))) ((null target-tail)) (il:movd? target save) (il:movd source target)))))
(defun deinstall-boxed-fns nil (install-boxed-fns t))
(il:declare\: il:docopy il:donteval@load il:donteval@compile
(install-boxed-fns)
)
(il:declare\: il:dontcopy il:doeval@compile
(il:declare\: il:doeval@compile il:dontcopy
(il:localvars . t)
)
)
(xcl:define-file-environment "DORADO-CMLFLOAT" :package "CL" :readtable "XCL" :compiler compile-file)
(il:putprops il:dorado-cmlfloat il:copyright ("Xerox Corporation" 1988))
(il:declare\: il:dontcopy
(il:filemap (nil)))
il:stop