;;; 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))))