;;; This is a -*-Lisp-*- file.
;;;
;;; **********************************************************************
;;; This code was written as part of the Spice Lisp project at
;;; Carnegie-Mellon University, and has been placed in the public domain.
;;; Spice Lisp is currently incomplete and under active development.
;;; If you want to use this code or any part of Spice Lisp, please contact
;;; Scott Fahlman (FAHLMAN@CMUC).
;;; **********************************************************************
;;;
;;; Functions to implement irrational functions for Spice Lisp
;;; Maintained by Steven Handerson.
;;;
;;; The irrational functions are part of the standard Spicelisp environment.
;;;
;;; **********************************************************************
;;;
;;; Note: long-float constants are in the file lfloat-consts.
;;;
(eval-when (compile)
(defmacro **short-float** (f e)
`(scale-float (decode-float (coerce ,f 'short-float)) ,e))
)
;;; From perqinit.
(defconstant most-positive-single-float
(%primitive make-immediate-type most-positive-fixnum 18))
(defconstant most-negative-single-float
(%primitive make-immediate-type 1 19))
;; This should be fixed whenever the printer is.
(defconstant least-positive-single-float
(%primitive make-immediate-type 1 18))
(defconstant least-negative-single-float
(%primitive make-immediate-type most-positive-fixnum 19))
(defconstant %single-float-exponent-length 8)
(defconstant %single-float-mantissa-length 20)
(defconstant single-float-epsilon (**short-float** 1 -20))
(defconstant single-float-negative-epsilon (**short-float** 1 -21))
(defconstant most-positive-short-float most-positive-single-float)
(defconstant least-positive-short-float least-positive-single-float)
(defconstant least-negative-short-float least-negative-single-float)
(defconstant most-negative-short-float most-negative-single-float)
(defconstant short-float-epsilon single-float-epsilon)
(defconstant short-float-negative-epsilon single-float-negative-epsilon)
(defconstant %double-float-exponent-length 11)
(defconstant %double-float-mantissa-length 53)
(defconstant %short-float-one 1.0)
(defconstant %short-float-half .5)
(defconstant %short-sqrt-half (**short-float** #o5520236 0))
(defconstant %short-pi 3.1415926535897932384)
(defconstant %short-pi/2 1.5707963267948966192)
(defconstant %short-pi/4 .7853981633974)
(defconstant %short-e 2.71828187
"Euler's number, short float version.")
(defconstant sixth-pi 0.5236)
(defconstant sqrt3 1.7305248)
(defconstant 2-sqrt3 0.26794816)
(defconstant inv-2-sqrt3 3.7320704)
(defconstant %short-tan-max-x 53526.413)
(defconstant %short-2/pi .6366197723)
(defconstant %short-tan-eps (**short-float** 1 -10))
(defconstant %short-tan-c1 (**short-float** #o1444 1))
(Defconstant %short-tan-c2 4.838067948s-4)
(defconstant %short-tan-p1 -.957017723s-1)
(defconstant %short-tan-q0 1.0s0)
(defconstant %short-tan-q1 -.429135777s0)
(defconstant %short-tan-q2 .971685835s-2)
(defconstant %short-asin-eps (**short-float** 1 -10))
(defconstant %short-asin-p1 .0833331098s0)
(defconstant %short-asin-p2 -.0450065898s0)
(defconstant %short-asin-q0 .5)
(defconstant %short-asin-q1 -.4950779396s0)
(defconstant %short-asin-q2 .0892278749s0)
(defconstant %short-log-base2-e 1.44269504088)
(defconstant %short-exp-c1 (**short-float** #o543 0))
(defconstant %short-exp-c2 -2.1219444005e-4)
(defconstant %short-exp-p0 .2499999995)
(defconstant %short-exp-p1 .4160288626e-2)
(defconstant %short-exp-q1 .49987178778e-1)
(defconstant %short-exp-array
'#(1.0
1.09050773
1.1892071
1.29683954
1.41421354
1.54221082
1.68179283
1.83400811)
"Used in the exp function. (aref %short-exp-array i) <-> (expt 2.0 (/ i 8))")
(defconstant %short-log-C0 0.7071067811)
(defconstant %short-log-C1 (**short-float** #o543 0))
(defconstant %short-log-C2 -0.0002121944400546905827679)
(defconstant %short-float-half .5)
(defconstant %short-log-a0 -0.5527074855)
(defconstant %short-log-b0 -6.632718214)
(defconstant %short-log-b1 1.0)
(defconstant %short-sin-maximum 1.0s10)
(defconstant %short-sin-epsilon (**short-float** 1 -10))
(defconstant %short-c2-sin 0.000967653589793)
(defconstant %short-sin-r1 -0.1666665668)
(defconstant %short-sin-r2 0.008333025139)
(defconstant %short-sin-r3 -0.0001980741872)
(defconstant %short-sin-r4 0.000002601903036)
(defconstant %short-zero 0.0)
(defconstant %short-c1-sin 3.140625)
(defconstant %short-expt-a1
(make-array 17. :initial-contents
(list (**short-float** 1 1)
(**short-float** #o7522256 0)
(**short-float** #o7254030 0)
(**short-float** #o7014632 0)
(**short-float** #o6564236 0)
(**short-float** #o6342220 0)
(**short-float** #o6126344 0)
(**short-float** #o5720424 0)
(**short-float** #o5520236 0)
(**short-float** #o5325406 0)
(**short-float** #o5137732 0)
(**short-float** #o4757246 0)
(**short-float** #o4603376 0)
(**short-float** #o4434172 0)
(**short-float** #o4271270 0)
(**short-float** #o4132530 0)
(**short-float** #o4 0))))
(defconstant %short-expt-a2
(make-array 8. :initial-contents
(list (**short-float** #o1505222 #o-24)
(**short-float** #o1673024 #o-24)
(**short-float** #o1405216 #o-24)
(**short-float** #o347654 #o-26)
(**short-float** #o1672440 #o-24)
(**short-float** #o230110 #o-26)
(**short-float** #o334724 #o-26)
(**short-float** #o331746 #o-26))))
(defconstant %short-expt-p1 .83357541s-1)
(defconstant %short-expt-q1 .69314675s0)
(defconstant %short-expt-q2 .24018510s0)
(defconstant %short-expt-q3 .54360383s-1)
(defconstant %short-expt-k .44269504s0)
;;; Warning: multiple evaluation in these macros is the rule,
;;; not the exception.
;;; Some floating macros.
(eval-when (compile)
(defmacro float-op (a b short-float long-float)
`(cond ((typep ,a 'long-float) (,long-float ,a (%primitive float-long ,b)))
((typep ,b 'long-float) (,long-float (%primitive float-long ,a) ,b))
(t (,short-float (%primitive float-short ,a) (%primitive float-short ,b)))))
(defmacro floatize (f)
`(if (floatp ,f) ,f (setq ,f (float ,f))))
(defmacro scale-float-by (e f) `(scale-float ,f ,e)))
(eval-when (compile load)
(defmacro poly-eval (x &rest list)
(do ((list (cdr (reverse list)) (cdr list))
(form (car (last list)) `(+ ,(car list) (* ,x ,form))))
((null list) form))))
(defun %signum (x)
(let ((res (cond ((plusp x) 1)
((zerop x) 0)
(t -1))))
(if (floatp x) (float res x) res)))
;;; Integer square root - (<= (expt result 2) input).
;;; Performs in (log input) time.
;;; Successively approximates the result using two bounds and their average,
;;; repeated until the bounds differ by at most 1 (for input=1, both start
;;; equal).
(defun isqrt (x)
(if (and (integerp x) (not (minusp x)))
(if (zerop x) 0
(do* ((p () (<= (* m m) x))
(b 1 (if p m b))
(h x (if p h m))
(m (ash x -1) (ash (+ b h) -1)))
((<= h (1+ b)) b)))
(error "Isqrt: ~S argument must be a nonnegative integer" x)))
;;; Assumes that the argument is a nonnegative long-float.
;;; Uses Newton's method with a special initial estimate, described
;;; in Cody and Waite.
(eval-when (compile)
(defmacro short-sqrt (x)
`(multiple-value-bind (f N) (decode-float ,x)
(let ((y (+ (* f .59016) .41731)))
(dotimes (i 2) (setq y (scale-float (+ y (/ f y)) -1)))
(cond ((oddp N)
(setq y (* y %short-sqrt-half))
(setq N (1+ N))))
(scale-float y (/ N 2))))))
(eval-when (compile)
(defmacro long-sqrt (x)
`(multiple-value-bind (f N) (decode-float ,x)
(let ((y (+ %long-sqrt-c1 (* %long-sqrt-c2 f))))
(dotimes (i 3) (setq y (scale-float (+ y (/ f y)) -1)))
(cond ((oddp N)
(setq y (* y %long-sqrt-half))
(setq N (1+ N))))
(scale-float y (/ N 2))))))
(defun sqrt (x)
(cond ((minusp x) (error "~S Argument out of range." x)))
(typecase x
(short-float (short-sqrt x))
(long-float (long-sqrt x))
(number (sqrt (float x)))
(t (error "Sqrt: Argument not number - ~A" x))
))
;;; Function calculates the value of x raised to the nth power.
;;; This function calculates the successive squares of base,
;;; storing them in newbase, halving n at the same time. If
;;; n is odd total is multiplied by base, at any given time (fix later)
;(declare (inline intexp))
(defun intexp (base power)
(cond ((minusp power)
(/ (intexp base (- power))))
((and (rationalp base)
(not (integerp base))) ; Could be a make-rational.
(/ (intexp (numerator base) power) (intexp (denominator base) power)))
((and (integerp base) (= base 2))
(ash 1 power))
(t (do ((nextn (ash power -1) (ash power -1))
(total (if (oddp power) base 1)
(if (oddp power) (* base total) total)))
((zerop nextn) total)
(setq base (* base base))
(setq power nextn)))))
;;; 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.
(defun expt (x n)
(cond ((and (rationalp x) (integerp n)) (intexp x n))
((zerop x) (if (plusp n) x
(error "~A to a non-positive power ~A." x n)))
((and (not (integerp n)) (minusp x))
(error "Negative number ~A to non-integral power ~A." x n))
(t (float-op x n short-expt long-expt))))
(eval-when (compile)
(defmacro short-exp (x)
`(let* ((n (round (* ,x %short-log-base2-e)))
(xn (float n 1.0))
(g (- (- ,x (* xn %short-exp-c1)) (* xn %short-exp-c2)))
(z (* g g))
(g*Pz (* g (poly-eval z %short-exp-p0 %short-exp-p1))))
(scale-float-by
(1+ n) (+ .5 (/ g*Pz (- (poly-eval z .5 %short-exp-q1) g*Pz)))))))
(eval-when (compile)
(Defmacro long-exp (x)
`(let* ((n (round (* ,x %long-log-base2-e)))
(xn (float n %long-log-base2-e))
(g (- (- ,x (* xn %long-exp-c1)) (* xn %long-exp-c2)))
(z (* g g))
(g*Pz (* g (poly-eval z %long-exp-p0 %long-exp-p1 %long-exp-p2))))
(scale-float-by
(1+ n) (+ %long-half
(/ g*Pz (- (poly-eval z %long-half %long-exp-q1 %long-exp-q2)
g*Pz)))))))
(defun exp (num)
"Calculates e to the given power, where e is the base of natural logarithms."
(if (integerp num) (setq num (float num)))
(typecase num
(short-float (short-exp num))
(long-float (long-exp num))
(number (exp (float num)))
(t (error "Exp: object not a number - ~A" num))))
;;; Hyperbolic trig functions.
;;; Each of the hyperbolic trig functions is calculated directly from
;;; their definition. Exp(x) is calculated only once for efficiency.
(defun sinh (x)
(let ((z (exp x)))
(/ (- z (/ z)) 2)))
(defun cosh (x)
(let ((z (exp x)))
(/ (+ z (/ z)) 2)))
;;; Different form than in the manual. Does much better.
(defun tanh (x)
(let* ((z (exp (* 2 x)))
(y (/ z)))
(- (/ (1+ y)) (/ (1+ z)))))
(eval-when (compile)
(Defmacro short-log (x)
`(cond ((plusp ,x)
(multiple-value-bind (f n) (decode-float ,x)
(declare (short-float f) (fixnum n))
(let (znum zden)
(cond ((> f %short-log-C0)
(setq znum (- (- f .5) .5))
(setq zden (+ .5 (* .5 f))))
(t (decf N)
(setq zden
(+ .5
(* .5
(setq znum (- f .5)))))))
(let* ((z (/ znum zden))
(w (* z z))
(rz2 (/ (* w (poly-eval w %short-log-a0))
(poly-eval w %short-log-b0 %short-log-b1)))
(big-Rz (+ z (* z rz2)))
(XN (%primitive float-short N)))
(+ (+ (* XN %short-log-c2) big-Rz)
(* XN %short-log-c1))))))
((zerop ,x) (error "Log of ~A undefined." ,x))
(t (error "Complex numbers not yet implemented - (log ~A)" ,x))))
(defmacro long-log (l)
`(cond ((zerop ,l) (error "Log: log of zero not defined"))
((plusp ,l)
(multiple-value-bind (f N) (decode-float ,l)
(let (znum zden)
(cond ((> f %long-sqrt-half)
(setq zden (+ %long-half (* %long-half f)))
(setq znum (- (- f %long-half) %long-half)))
(t (decf N)
(setq zden
(+ %long-half (* %long-half
(setq znum (- f %long-half)))))))
(let* ((z (/ znum zden))
(w (* z z))
(rz↑2 (/ (* w (poly-eval w %long-log-a0 %long-log-a1
%long-log-a2))
(poly-eval w %long-log-b0
%long-log-b1 %long-log-b2
%long-float-one)))
(bigrz (+ z (* z rz↑2)))
(xn (%primitive float-long N)))
(+ (+ (* xn %long-log-c2) bigrz)
(* xn %long-log-c1))))))
(t (error "Complex numbers not yet implemented - (log ~A)" ,l))))
)
(eval-when (compile)
(defmacro log-base (x base)
`(/ (log ,x) (log ,base))))
(defun log (number &optional (base nil base-supplied))
(cond (base-supplied (float-op number base log-base log-base))
(t (typecase number
(short-float (short-log number))
(long-float (long-log number))
(number (log (float number)))
(t (error "Log: argument not number - ~A" number))))))
(defun asinh (x)
(log (+ x (sqrt (+ (* x x) 1)))))
(defun acosh (x)
(if (plusp x)
(log (+ x (sqrt (- (* x x) 1))))
(error "~S argument out of range." x)))
(defun atanh (x)
(if (< -1 x 1)
(* 0.5 (log (/ (1+ x) (- 1 x))))
(error "~S argument out of range." x)))
(Eval-when (compile)
(Defmacro short-sin (x cos-flag)
`(let ((sgn ,(if cos-flag 1 `(if (minusp ,x) -1 1)))
(y ,(if cos-flag `(+ (abs ,x) %short-pi/2) `(abs ,x))))
(if (> y %short-sin-maximum)
,(if cos-flag
`(cerror "Cos: absolute value of argument, ~s, too large."
"Will return an imprecise result." ,x)
`(cerror "Sin: absolute value of argument, ~s, too large."
"Will return an imprecise result." ,x)))
(let* ((n (round Y %short-pi))
(xn (float n))
(f (- (- Y (* xn %short-c1-sin))
(* xn %short-c2-sin))))
(if (oddp n) (setq sgn (- sgn)))
(* (if (< (abs f) %short-sin-epsilon) f
(let ((g (* f f)))
(+ f (* f (* g (poly-eval g %short-sin-r1 %short-sin-r2
%short-sin-r3 %short-sin-r4))))))
sgn))))
(defmacro long-sin (x cos-flag)
`(let ((sgn ,(if cos-flag 1 `(if (minusp ,x) -1 1)))
(y ,(if cos-flag `(+ (abs ,x) %long-pi/2) `(abs ,x))))
(if (> y %long-sin-max)
,(if cos-flag
`(cerror "Cos: absolute value of argument, ~S too large."
"Will return an imprecise result." ,x)
`(cerror "Sin: absolute value of argument, ~S too large."
"Will return an imprecise result." ,x)))
(let* ((N (round Y pi))
(XN (float N %long-half))
(f (- (- Y (* XN %long-sin-c1))
(* XN %long-sin-c2))))
(if (oddp N) (setq sgn (- sgn)))
(* (cond ((< (abs f) %long-sin-epsilon) f)
(t (let* ((g (* f f)))
(+ f (* f (* g (poly-eval g %long-sin-r1 %long-sin-r2
%long-sin-r3 %long-sin-r4
%long-sin-r5 %long-sin-r6
%long-sin-r7 %long-sin-r8)))))))
sgn))))
)
(defun sin (x)
(if (not (floatp x)) (setq x (float x)))
(typecase x
(short-float (short-sin x nil))
(long-float (long-sin x nil))))
(defun cos (x)
(if (not (floatp x)) (setq x (float x)))
(typecase x
(short-float (short-sin x T))
(long-float (long-sin x T))))
;;; The arc tangent is determined using a power series described in
;;; Computer Approximations by Hart & Cheney p122 (6.5.15)
;;; As the radius of convergence is 1 there is a very poor performance
;;; When the value of x is near 1. This should be fixed at a later date.
;(declare (inline short-atan))
(defun atan1 (x)
(cond ((minusp x) (- (atan1 (- x))))
((< x 2-sqrt3) (- %short-pi/2 (atan2 (/ x))))
((< x 1)
(- (+ %short-pi/2 sixth-pi) (atan2 (/ (+ sqrt3 x) (1- (* x sqrt3))))))
((< x inv-2-sqrt3)
(let* ((inv (/ x)))
(- (atan2 (/ (+ sqrt3 inv) (1- (* sqrt3 inv)))) sixth-pi)))
(t (atan2 x))))
(defun atan2 (x)
(do* ((sqr (- (/ (* x x))))
(int 1 (+ 2 int))
(old 0 val)
(pow (/ x) (* pow sqr))
(val pow (+ val (/ pow int))))
((= old val)
(- %short-pi/2 val))))
(defun short-atan (y &optional x)
(cond ((not x)
(if (zerop y) 0.0 (atan1 y)))
((= y x 0)
(error "Error in double entry atan both 0." ))
((= x 0) (* (%signum y) %short-pi/2))
((= y 0) (if (plusp x) 0 %short-pi))
((and (plusp x)(plusp y)) (atan1 (/ y x)))
((plusp y) (- %short-pi (atan1 (/ (- y) x))))
((plusp x) (- (atan1 (/ y (- x)))))
(t (- (atan1 (/ y x)) %short-pi))))
;;; Assume both args are floats of the same type.
(defun long-atan (x &optional (y %long-float-one))
(cond ((zerop x)
(if (zerop y) (error "Attempt to divide ~S by ~S" x y) %long-pi/2))
;; Check for underflow and overflow in the division.
((multiple-value-bind (xf xe) (decode-float x)
(multiple-value-bind (yf ye) (decode-float y)
(multiple-value-bind (res rese)
(decode-float (/ xf yf))
(declare (ignore res))
(let ((e (+ (- xe ye) rese)))
(cond ((> e %long-float-exp-max) %long-pi/2)
((< e %long-float-exp-min) 0)
(t (setq x (/ x y))
nil)))))))
(t (let* ((f (abs x))
(N (cond ((> f %long-float-one) (setq f (/ f)) 2)
(t 0))))
(cond ((> f %long-atan-2-sqrt3)
(setq f (/ (1- (* f %long-atan-sqrt3))
(+ f %long-atan-sqrt3)))
(setq N (1+ N))))
(let ((res (cond ((< (abs f) %long-atan-epsilon) f)
(t (let* ((g (* f f))
(R (/ (* g (poly-eval g %long-atan-p0
%long-atan-p1
%long-atan-p2
%long-atan-p3))
(poly-eval g %long-atan-q0
%long-atan-q1
%long-atan-q2
%long-atan-q3
%long-atan-q4))))
(+ f (* f R)))))))
(if (> N 1) (setq res (- res)))
(setq res (+ res (svref %long-atan-vector N)))
(if (minusp x) (- res) res))))))
(defun atan (x &optional y)
(cond (y (if (or (complexp x) (complexp y))
(error "~S and ~S both complex." x y))
(float-op x y short-atan long-atan))
(t (typecase (floatize x)
(short-float (short-atan x))
(long-float (long-atan x))))))
;;;
;;; Short-float expt. Still computes (exp (* y (log x))), but
;;; does so as to minimize loss of precision to exp,
;;; since normally the error in the result depends upon the
;;; magnitude of the base.
;;;
;;; Used in the expt function. "The main requirement is that
;;; (< (abs (- x (reduce x))) 1/16)".
(eval-when (compile load)
(defmacro short-expt-reduce (x)
`(multiple-value-bind (f e s) (decode-float ,x)
(if (zerop f) f ;To prevent making thing.
(let ((off (- 16 e)))
(* s (scale-float
(if (plusp off)
(let ((thing (scale-float .5s0 off)))
(- (+ f thing) thing))
f)
e)))))))
;;; Not for negative or zero base.
;;;
;;; Algorithm from Cody and Waite.
;;; Altered for (gee wow) zero-based array accessing.
;;;
(defun short-expt (x y)
(multiple-value-bind (g m) (decode-float x)
(let ((p 1))
(if (<= g (svref %short-expt-A1 8)) (setq p 9))
(if (<= g (svref %short-expt-A1 (+ p 3))) (setq p (+ p 4)))
(if (<= g (svref %short-expt-A1 (1+ p))) (setq p (+ p 2)))
(let* ((z (scale-float-by
1 (/ (- (- g (svref %short-expt-a1 p))
;; 1- for 0-based array access.
(svref %short-expt-a2 (1- (/ (1+ p) 2))))
(+ g (svref %short-expt-a1 p)))))
(v (* z z))
(R (* z v (poly-eval v %short-expt-p1)))
(U2 (+ (+ (+ R (* R %short-expt-K))
(* z %short-expt-K))
z))
(U1 (scale-float
(%primitive float-short (- (* 16 m) p)) -4))
(Y1 (short-expt-reduce y))
(Y2 (- Y y1))
(W (+ (* u2 y) (* u1 y2)))
(W1 (short-expt-reduce W))
(W2 (- W W1)))
(setq W (+ W1 (* U1 Y1)))
(setq W1 (short-expt-reduce W))
(setq W2 (+ W2 (- W W1)))
(setq W (short-expt-reduce W2))
(let ((IW1 (truncate (scale-float-by 4 (+ W1 W)))))
(setq W2 (- W2 W))
; (cond ((> IW1 %short-expt-max)
; (error "~S to the ~S has overflowed." x y))
; ((< IW1 %short-expt-min)
; (error "~S to the ~S has underflowed." x y)))
(cond ((plusp W2) (incf IW1) (decf W2 .0625s0)))
(let* ((mprime (if (minusp IW1)
(truncate IW1 16)
(1+ (truncate IW1 16))))
(pprime (- (* 16 mprime) IW1)))
(scale-float-by
mprime (+ (svref %short-expt-a1 pprime)
(* (svref %short-expt-a1 pprime)
W2 (poly-eval W2 %short-expt-q1 %short-expt-q2
%short-expt-q3))))))))))
;;;
;;; Same thing, for long-floats.
;;;
;;; Truncating by bugout is slower, although plusp...
(eval-when (compile load)
(defmacro long-expt-reduce (x)
`(multiple-value-bind (f e s) (decode-float ,x)
(let ((off (- 49 e)))
(* s (scale-float
(if (plusp off)
(let ((thing (scale-float %long-half off)))
(- (+ f thing) thing))
f)
e))))))
;;; Not for negative or zero base.
;;;
;;; Algorithm from Cody and Waite.
;;; Altered for (gee wow) zero-based array accessing.
;;;
(defun long-expt (x y)
(multiple-value-bind (g m) (decode-float x)
(let ((p 1))
(if (<= g (svref %long-expt-A1 8)) (setq p 9))
(if (<= g (svref %long-expt-A1 (+ p 3))) (setq p (+ p 4)))
(if (<= g (svref %long-expt-A1 (1+ p))) (setq p (+ p 2)))
(let* ((z (scale-float-by
1 (/ (- (- g (svref %long-expt-a1 p))
;; 1- for 0-based array access.
(svref %long-expt-a2 (1- (/ (1+ p) 2))))
(+ g (svref %long-expt-a1 p)))))
(v (* z z))
(R (* z v (poly-eval v %long-expt-p1 %long-expt-p2
%long-expt-p3 %long-expt-p4)))
(U2 (+ (+ (+ R (* R %long-expt-K))
(* z %long-expt-K))
z))
(U1 (scale-float
(%primitive float-long (- (* 16 m) p)) -4))
(Y1 (long-expt-reduce y))
(Y2 (- Y y1))
(W (+ (* u2 y) (* u1 y2)))
(W1 (long-expt-reduce W))
(W2 (- W W1)))
(setq W (+ W1 (* U1 Y1)))
(setq W1 (long-expt-reduce W))
(setq W2 (+ W2 (- W W1)))
(setq W (long-expt-reduce W2))
(let ((IW1 (truncate (scale-float-by 4 (+ W1 W)))))
(setq W2 (- W2 W))
; (cond ((> IW1 %long-expt-max)
; (error "~S to the ~S has overflowed." x y))
; ((< IW1 %long-expt-min)
; (error "~S to the ~S has underflowed." x y)))
(cond ((plusp W2) (incf IW1) (decf W2 %long-float-sixteenth)))
(let* ((mprime (if (minusp IW1)
(truncate IW1 16)
(1+ (truncate IW1 16))))
(pprime (- (* 16 mprime) IW1)))
(scale-float-by
mprime (+ (svref %long-expt-a1 pprime)
(* (svref %long-expt-a1 pprime)
W2 (poly-eval W2 %long-expt-q1 %long-expt-q2
%long-expt-q3 %long-expt-q4
%long-expt-q5 %long-expt-q6
%long-expt-q7))))))))))
(defun short-asin (x flag)
(declare (short-float x) (fixnum flag))
(let* ((y (abs x))
(test (> y %short-float-half))
(i (cond (test (if (> y %short-float-one)
(error "~A - argument out of range" x))
(1- flag))
(t flag)))
(result (if (and (not test) (< y %short-asin-eps))
y
(let* ((g (if test
(let ((temp (scale-float
(- %short-float-one y) -1)))
(setq y (- (scale-float (sqrt temp) 1)))
temp)
(* y y))))
(+ y (* y (/ (* g (poly-eval g %short-asin-p1
%short-asin-p2))
(poly-eval g %short-asin-q0 %short-asin-q1
%short-asin-q2))))))))
(cond ((zerop flag)
(let ((res (if (zerop i)
result
(+ %short-pi/4 result %short-pi/4))))
(if (minusp x) (- res) res)))
((minusp x)
(if (zerop i) (+ %short-pi/2 result %short-pi/2)
(+ %short-pi/4 result %short-pi/4)))
((zerop i) (- result))
(t (+ (- %short-pi/4 result) %short-pi/4)))))
(defun long-asin (x flag)
(declare (long-float x) (fixnum flag))
(let* ((y (abs x))
(test (> y %long-float-half))
(i (cond (test (if (> y %long-float-one)
(error t nil "argument out of range"))
(1- flag))
(t flag)))
(result (if (and (not test) (< y %long-asin-eps))
y
(let* ((g (if test
(let ((temp (scale-float
(- %long-float-one y) -1)))
(setq y (- (scale-float (sqrt temp) 1)))
temp)
(* y y))))
(+ y (* y (/ (* g (poly-eval g %long-asin-p1
%long-asin-p2 %long-asin-p3
%long-asin-p4
%long-asin-p5))
(poly-eval g %long-asin-q0 %long-asin-q1
%long-asin-q2 %long-asin-q3
%long-asin-q4
%long-asin-q5))))))))
(cond ((zerop flag)
(let ((res (if (zerop i)
result
(+ %long-pi/4 result %long-pi/4))))
(if (minusp x) (- res) res)))
((minusp x)
(if (zerop i) (+ %long-pi/2 result %long-pi/2)
(+ %long-pi/4 result %long-pi/4)))
((zerop i) (- result))
(t (+ (- %long-pi/4 result) %long-pi/4)))))
(Defun asin (x)
(typecase x
(long-float (long-asin x 0))
(short-float (short-asin x 0))
(number (asin (float x)))))
(defun acos (x)
(typecase x
(long-float (long-asin x 1))
(short-float (short-asin x 1))
(number (acos (float x)))))
(defun short-tan (x)
(declare (short-float x))
(let ((y (abs x)))
(declare (short-float y))
(if (< %short-tan-max-x y)
(error "~s out of range." x)
(let* ((N (round (* x %short-2/pi)))
(xn (%primitive float-short N))
(f (- (- x (* xn %short-tan-c1))
(* xn %short-tan-c2))))
(if (< (abs f) %short-tan-eps)
(if (evenp n) f (- (/ f)))
(let* ((g (* f f))
(p (+ f (* f g (poly-eval g %short-tan-p1))))
(q (poly-eval g %short-tan-q0 %short-tan-q1
%short-tan-q2)))
(declare (short-float g p q))
(if (evenp N)
(/ p q)
(/ q (- p)))))))))
(defun long-tan (x)
(declare (long-float x))
(let ((y (abs x)))
(declare (long-float y))
(if (< %long-tan-max-x y)
(error "~s out of range." x)
(let* ((N (round (* x %long-2/pi)))
(xn (%primitive float-long N))
(f (- (- x (* xn %long-tan-c1))
(* xn %long-tan-c2))))
(if (< (abs f) %long-tan-eps)
(if (evenp n) f (- (/ f)))
(let* ((g (* f f))
(p (+ f (* f g (poly-eval g %long-tan-p1 %long-tan-p2
%long-tan-p3))))
(q (poly-eval g %long-tan-q0 %long-tan-q1
%long-tan-q2 %long-tan-q3 %long-tan-q4)))
(declare (long-float g p q))
(if (evenp N)
(/ p q)
(/ q (- p)))))))))
(Defun tan (x)
(if (not (floatp x)) (setq x (float x)))
(typecase x
(long-float (long-tan x))
(short-float (short-tan x))))