head     1.1;
access   ;
symbols  ;
locks    ; strict;
comment  @@;


1.1
date     90.02.22.16.43.50;  author ram;  state Exp;
branches ;
next     ;


desc
@Working code version as of 2/22/90.
@



1.1
log
@Initial revision
@
text
@;;; -*- Log: code.log; Package: Lisp -*-
;;;
;;; **********************************************************************
;;; This code was written as part of the Spice Lisp project at
;;; Carnegie-Mellon University, and has been placed in the public domain.
;;; 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.
;;;
(in-package 'lisp)
(export '(most-positive-short-float
	  least-positive-short-float least-negative-short-float
	  most-negative-short-float most-positive-single-float
	  least-positive-single-float least-negative-single-float
	  most-negative-single-float short-float-epsilon single-float-epsilon
	  short-float-negative-epsilon single-float-negative-epsilon
	  exp expt log sqrt isqrt
	  sin cos tan
	  asin acos atan
	  sinh cosh tanh asinh acosh atanh))


(in-package 'system)
(export '(short-float-radix single-float-radix
			    double-float-radix long-float-radix))
(in-package "EXTENSIONS")
(export '*intexp-maximum-exponenent*)

(in-package 'lisp)


(eval-when (compile)
(defmacro **short-float** (f e)
  `(scale-float
    (multiple-value-bind (m e s)
			 (decode-float (coerce ,f 'short-float))
      (declare (ignore e))
      (* m s))
    ,e))

(defmacro **single-float** (f e)
  `(scale-float
    (multiple-value-bind (m e s)
			 (decode-float (coerce ,f 'single-float))
      (declare (ignore e))
      (* m s))
    ,e))
)

;;; From perqinit.

(defconstant most-positive-single-float
  (%primitive make-immediate-type #x07f7ffff %short-+-float-type)
  "Positive single floating-point number closest to positive infinity.")
(defconstant most-negative-single-float
  (%primitive make-immediate-type #x07F7ffff %short---float-type)
  "Negative single floating-point number closest to negative infinity.")
;; This should be fixed whenever the printer is.
(defconstant least-positive-single-float
  (%primitive make-immediate-type #x00080000 %short-+-float-type)
  "Positive single floating-point number closest (but not equal) to zero.")
(defconstant least-negative-single-float
  (%primitive make-immediate-type #x00080000 %short---float-type)
  "Negative single floating-point number closest (but not equal) to zero.")

(defconstant %single-float-exponent-length 8)
(defconstant %single-float-mantissa-length 20)
(defconstant single-float-epsilon (**short-float** 1 -20)
  "Smallest positive single floating-point number (e) such that
   (not (= (float 1 e) (+ (float 1 e) e))) returns true.")
(defconstant single-float-negative-epsilon (**short-float** -1 -21)
  "Smallest positive single floating-point number (e) such that 
   (not (= (float 1 e) (- (float 1 e) e))) returns true.")

(defconstant most-positive-short-float most-positive-single-float
  "Positive short floating-point number closest to positive infinity.")
(defconstant least-positive-short-float least-positive-single-float
  "Positive short floating-point number closest (but not equal) to zero.")
(defconstant least-negative-short-float least-negative-single-float
  "Negative short floating-point number closest (but not equal) to zero.")
(defconstant most-negative-short-float most-negative-single-float
  "Negative single floating-point number closest to negative infinity.")
(defconstant short-float-epsilon single-float-epsilon
  "Smallest positive short floating-point number (e) such that
   (not (= (float 1 e) (+ (float 1 e) e))) returns true.")
(defconstant short-float-negative-epsilon single-float-negative-epsilon
  "Smallest positive short floating-point number (e) such that 
   (not (= (float 1 e) (- (float 1 e) e))) returns true.")

(defconstant %double-float-exponent-length 11)
(defconstant %double-float-mantissa-length 53)
(defconstant %long-float-exp-min (multiple-value-bind (f e)
				   (decode-float least-positive-long-float)
				   (declare (ignore f))
				   (- e 3)))
(defconstant %long-float-exp-max (multiple-value-bind (f e)
				   (decode-float most-positive-long-float)
				   (declare (ignore f))
				   (+ e 3)))

(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 single-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))
#|
	   ((typep ,a 'single-float)
	    (,single-float ,a (%primitive float-single ,b)))
	   ((typep ,b 'single-float)
	    (,single-float (%primitive float-single ,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))))

(proclaim '(inline %signum %multbyi %negate-complex))
(defun %signum (x)
  (let ((res (cond ((plusp x) 1)
		   ((zerop x) 0)	
		   (t -1))))
    (if (floatp x) (float res x) res)))


;;;Useful functions for complex numbers

;magnitude of a complex number
(defun %complex-mag (z)
  (let* ((x (realpart z))
	 (y (imagpart z)))
    (sqrt (+ (* x x) (* y y)))))


;multiplying i (the square root of -1) times a number
(defun %multbyi (z)
  (%primitive make-complex (- (imagpart z)) (realpart z)))

(defun %negate-complex (z)
  (%primitive make-complex (- (realpart z)) (- (imagpart z))))


;;; 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)))
	 (cond ((oddp N)
		(setq y (* y %short-sqrt-half))
		(setq N (1+ N))))
	 (setq y (scale-float y (/ N 2)))
 	 (dotimes (i 2) (setq y (scale-float (+ y (/ ,x y)) -1)))
	 y))))

#|
(eval-when (compile)
  (defmacro single-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)
  "Returns the principle square root of x."
  (system:%primitive sqrt x))

(defun %sp-sqrt-short (x)
  (if (minusp x)
      (complex 0.0 (short-sqrt (abs x)))
      (short-sqrt x)))

(defun %sp-sqrt-long (x)
  (if (minusp x)
      (complex 0.0l0 (long-sqrt (abs x)))
      (long-sqrt x)))

;the square root of a complex number returns two roots
;this function returns the primary root
(defun %sp-sqrt-complex (z)
  (* (exp (/ (%multbyi (phase z)) 2)) (sqrt (%complex-mag z))))

    
;;; 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))

(defparameter *intexp-maximum-exponent* 10000)

(defun intexp (base power)
  (when (> (abs power) *intexp-maximum-exponent*)
    (cerror "Continue with calculation."
	    "The absolute value of ~S exceeds ~S."
	    power '*intexp-maximum-exponent* 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)))))



;compute (complex)^n where n is an integer
;some round-off error if n is not a fixnum

(defun complex-expt (z n)
  (* (expt (%complex-mag z) n) (cis (* n (phase z)))))

;this function computes z ^ w where w is a complex number
;it can also be used for any positive real number.

(defun complex-expt-compow (z w)
  (exp (* w (log z))))


;;; 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)
  "Returns x raised to the nth power."
  (cond ((= n 0)
	 (setq x (realpart x))
	 (cond ((rationalp x) 1)
	       ((floatp x) (float 1.0 x))
	       (T 1)))
	((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 (complexp x) (integerp n)) (intexp x n))
	((complexp n) (if (or (complexp x) (plusp x)) (complex-expt-compow x n)
		  (error "~A negative number to a complex power ~A." x n)))
	((complexp x) (complex-expt x n))
	((integerp n) (intexp x n))
	((and (not (integerp n)) (minusp x))
	 (exp (* n (log x))))
	(t (float-op x n short-expt single-expt long-expt))))




(eval-when (compile)
  (defmacro short-exp (x)
    `(let* ((n (round (* ,x %short-log-base2-e)))
	    (xn (float n 1.0s0))
	    (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) (+ .5s0 (/ g*Pz (- (poly-eval z .5s0 %short-exp-q1) g*Pz)))))))


;;; This is bogus, but can stay until we are over to the new machine.
#|
(eval-when (compile)
  (defmacro single-exp (x)
    `(let* ((n (round (* ,x %short-log-base2-e)))
	    (xn (float n 1.0f0))
	    (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) (+ .5f0 (/ g*Pz (- (poly-eval z .5f0 %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."
  (%primitive exp num))

(defun %sp-exp-short (x)
  (short-exp x))

(defun %sp-exp-long (x)
  (long-exp x))


;compute exp(z) where z is complex
;is called by exp
(defun %sp-exp-complex (z)
  (let* ((x (realpart z))
	 (y (imagpart z)))
	 (* (exp x) (cis y))))) 

;;; 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.

(defun sinh (x)
  "Returns the hyperbolic sine of x."
  (let ((z (exp x)))
    (/ (- z (/ z)) 2)))

(defun cosh (x)
  "Returns the hyperbolic cosine of x."
  (let ((z (exp x)))
    (/ (+ z (/ z)) 2)))

;;; Different form than in the manual.  Does much better.
(defun tanh (x)
  "Returns the hyperbolic tangent of x."
  (let* ((z (exp (* 2 x)))
	 (y (/ z)))
    (- (/ (1+ y)) (/ (1+ z)))))


;;; Exp of a number.

(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 (complex (log (abs ,x)) %short-pi))))

#|
  (Defmacro single-log (x) 
    `(cond ((plusp ,x)
	    (multiple-value-bind (f n) (decode-float ,x)
	      (declare (single-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-single N)))
		  (+ (+ (* XN %short-log-c2) big-Rz)
		     (* XN %short-log-c1))))))
	   ((zerop ,x) (error "Log of ~A undefined." ,x))
	   (t (complex (log (abs ,x)) %short-pi))))
|#

  (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 (complex (log (abs ,l)) pi))))
)

;ATTENTION:
;log(exp(z)) <> z most of the time  when z is complex 
;exp(log(z)) = z all of the time.   when z is complex

(defun log (number &optional (base nil base-supplied))
  "Returns the natural logarithm of a number.  If an optional base
   is specified then it returns the log of the number to that base."
  (cond (base-supplied (/ (log number) (log base)))
	(T (%primitive log number))))

(defun %sp-log-short (x)
  (short-log x))

(defun %sp-log-long (x)
  (long-log x))

;natural log of a complex number
(defun %sp-log-complex (z)
  (complex (log (%complex-mag z)) (phase z)))

;;; Hyperbolic functions

;;; define the domain for complex hyperbolic functions

(defun in-asinh-domain (z)
  (cond ((< (- %short-pi/2) (imagpart z) %short-pi/2) t)
	((= (imagpart z) (- %short-pi/2))
	 (if (<= (realpart z) 0) t
	     nil))
	((= (imagpart z) %short-pi/2)
	 (if (>= (realpart z) 0) t
	     nil))
	(t nil)))


(defun in-acosh-domain (x)
  (cond ((plusp (realpart x))
		(if (and (<= (- %short-pi) x) (< x %short-pi))
		    t
		    nil))
	       ((= (realpart x) 0)
		(if (<= 0 (imagpart x) %short-pi)
		    t
		    nil))
	       (t  nil)))



(defun in-atanh-domain (z)
  (cond ((< (- %short-pi/2) (imagpart z) %short-pi/2) t)
	((= (imagpart z) (- %short-pi/2))
	 (if (minusp (realpart z)) t
	     nil))
	((= (imagpart z) %short-pi/2)
	 (if (plusp (realpart z)) t
	     nil))
	(t nil)))



(defun asinh (x)
  "Returns the hyperbolic arc sine of x."
  (cond ((complexp x)
	 (if (in-asinh-domain x)
	     (log (+ x (sqrt (+ (* x x) 1))))
	     (error "~S argument out of range." x)))
	(t (log (+ x (sqrt (+ (* x x) 1)))))))


(defun acosh (x) 
  "Returns the hyperbolic arc cosine of x."
  (cond ((complexp x) (if (in-acosh-domain x)
			  (log (+ x (sqrt (- (* x x) 1))))
			  (error "~S argument out of range." x)))
	((plusp x) (log (+ x (sqrt (- (* x x) 1)))))
	(t (error "~S argument out of range." x))))


(defun atanh (x)
  "Returns the hyperbolic arc tangent of x."
  (cond ((complexp x) (if (in-atanh-domain x)
	  (* 0.5 (log (/ (1+ x) (- 1 x))))
	  (error "~S argument out of range." x)))
	((< -1 x 1) (* 0.5 (log (/ (1+ x) (- 1 x)))))
	(t (error "~S argument out of range." x))))



;;; Support for sin and cos.

(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 (truncate (+ (/ Y %short-pi) 0.5)))
	      (xn (float n))
	      (f (- (- Y (* xn %short-c1-sin))
		    (* xn %short-c2-sin))))
	 (if (oddp n) (setq sgn (- sgn)))
	 (setq f (* (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))
	 (when (> f 1.0) (setq f 1.0))
	 (when (< f -1.0) (setq f -1.0))
	 f)))

#|
  (Defmacro single-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 (truncate (+ (/ Y %short-pi) 0.5)))
	      (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 (truncate (+ (/ Y pi) 0.5L0)))
	      (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))))
)


;;; Sine of a number.
(defun sin (x)
  "Returns the sine of x."
  (system:%primitive sin x))

;;; The following functions are invoked by the cos miscop when necessary.

(defun %sp-sin-short (x)
  (short-sin x nil))

(defun %sp-sin-long (x)
  (long-sin x nil))

(defun %sp-sin-complex (z)
  (let* ((x (realpart z))
	 (y (imagpart z)))
    (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))


;;; Cosine of a number.
(defun cos (x)
  "Returns the cosine of x."
  (system:%primitive cos x))


;;; The following functions are invoked by the cos miscop when necessary.

(defun %sp-cos-short (x)
  (short-sin x t))

(defun %sp-cos-long (x)
  (long-sin x t))

(defun %sp-cos-complex (z)
  (let* ((x (realpart z))
	 (y (imagpart z)))
    (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))

;;; 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.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))))


#|
(defun single-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))
  (let ((original-x x))
    (cond ((zerop x)
	   (if (zerop y)
	       (error "Attempt to divide ~S by ~S" x y)
	       (if (plusp y) 0.0l0 pi)))
	  ;; Check for underflow and overflow in the division.
	  ((zerop y)
	   (if (plusp x) %long-pi/2 (- %long-pi/2)))
	  ((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 ((exponent (+ (- xe ye) rese)))
		   (cond ((> exponent %long-float-exp-max) %long-pi/2)
			 ((< exponent %long-float-exp-min) 0.0l0)
			 (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 (plusp original-x) res (- res))))))))

(defun atan (x &optional y)
  "Returns the arc tangent of x.  If an optional y
   is given, then the arc tangent of y/x is returned."
  (cond ((null y)
	 (%primitive atan x))
	((or (complexp x) (complexp y))
	 (error "~S and ~S both complex." x y))
	(T (float-op x y short-atan single-atan long-atan))))

(defun %sp-atan-short (x)
  (short-atan x))

(defun %sp-atan-long (x)
  (long-atan x))

(defun %sp-atan-complex (z)
  (if (in-atan-domain z)
      (let ((i (%multbyi z)))
	(%negate-complex (%multbyi (* .5 (log (/ (+ 1 i) (- 1 i)))))))
      (error "Argument not in domain for atan. ~S" z)))


;define the domain for atan
(defun in-atan-domain (z)
  (cond ((< (- %short-pi/2) (realpart z) %short-pi/2) t)
	((= (realpart z) (- %short-pi/2))
	 (if (plusp (imagpart z)) t
	     nil))
	((= (realpart z) %short-pi/2)
	 (if (minusp (imagpart z)) t
	     nil))
	(t nil)))

;;;
;;; 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 ex s) (decode-float ,x)
       (if (zerop f) f				;To prevent making thing.
	   (let ((off (- 16 ex)))
	     (* s (scale-float
		   (if (plusp off)
		       (let ((thing (scale-float .5s0 off)))
			 (- (+ f thing) thing))
		       f)
		   ex)))))))


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



;;; Used in the expt function.  "The main requirement is that
;;; (< (abs (- x (reduce x))) 1/16)".
#|

(eval-when (compile load)
  (defmacro single-expt-reduce (x)
    `(multiple-value-bind (f ex s) (decode-float ,x)
       (if (zerop f) f				;To prevent making thing.
	   (let ((off (- 20 ex)))
	     (* s (scale-float
		   (if (plusp off)
		       (let ((thing (scale-float .5f0 off)))
			 (- (+ f thing) thing))
		       f)
		   ex)))))))

;;; Not for negative or zero base.
;;;
;;; Algorithm from Cody and Waite.
;;; Altered for (gee wow) zero-based array accessing.
;;;

(defun single-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-single (- (* 16 m) p)) -4))
	     (Y1 (single-expt-reduce y))
	     (Y2 (- Y y1))
	     (W (+ (* u2 y) (* u1 y2)))
	     (W1 (single-expt-reduce W))
	     (W2 (- W W1)))
	(setq W (+ W1 (* U1 Y1)))
	(setq W1 (single-expt-reduce W))
	(setq W2 (+ W2 (- W W1)))
	(setq W (single-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 ex s) (decode-float ,x)
       (if (zerop f) f
	   (let ((off (- 49 ex)))
	     (* s (scale-float
		   (if (plusp off)
		       (let ((thing (scale-float %long-half off)))
			 (- (+ f thing) thing))
		       f)
		   ex)))))))


;;; 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 not in domain" 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 single-asin (x flag) 
  (declare (single-float x) (fixnum flag))
  (let* ((y (abs x))
	 (test (> y %short-float-half))
	 (i (cond (test (if (> y %short-float-one)
			    (error "~A - argument not in domain" 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 not in domain"))
			(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)))))


;define the domains for the inverse trig functions

(defun in-asin-domain (z)
  (cond ((< (- %short-pi/2) (realpart z) %short-pi/2) t)
	((= (realpart z) (- %short-pi/2))
	 (if (>= (imagpart z) 0) t
	     nil))
	((= (realpart z) %short-pi/2)
	 (if (<= (imagpart z) 0) t
	     nil))
	(t nil)))

(defun in-acos-domain (z)
  (cond ((< 0 (realpart z) %short-pi) t)
	((= 0 (realpart z))
	 (if (>= (imagpart z) 0)
	     t
	     nil))
	  ((= (realpart z) %short-pi)
	   (if (<= (imagpart z) 0) 
	       t
	       nil))
	(t nil)))

(defun complex-asin (z)
  (if (in-asin-domain z)
      (let* ((i (%multbyi z)))
	(%negate-complex (%multbyi (log (+ i (sqrt (- 1 (* z z))))))))
      (error "Argument not in domain for asin. ~S" z)))


(defun complex-acos (z)
  (if (in-acos-domain z)
      (%negate-complex (%multbyi (log (+ z (%multbyi (sqrt (- 1 (* z z))))))))
      (error "Argument not in domain for acos. ~S" z)))


(defun asin (x)
  "Returns the arc sine of x."
  (typecase x
    (long-float (long-asin x 0))
;    (single-float (single-asin x 0))
    (short-float (short-asin x 0))
    (complex (complex-asin x))
    (number (asin (float x)))))

(defun acos (x)
  "Returns the arc cosine of x."
  (typecase x
    (long-float (long-asin x 1))
    (short-float (short-asin x 1))
    (complex (complex-acos x))
    (number (acos (float x)))))


(defun tan (x)
  "Returns the tangent of x."
  (%primitive tan x))

(defun %sp-tan-short (x)
  (declare (short-float x))
  (let ((y (abs x)))
    (declare (short-float y))
    (if (< %short-tan-max-x y)
	(error "~s not in domain." 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 single-tan (x)
  (declare (single-float x))
  (let ((y (abs x)))
    (declare (single-float y))
    (if (< %short-tan-max-x y)
	(error "~s not in domain." x)
	(let* ((N (round (* x %short-2/pi)))
	       (xn (%primitive float-single 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 (single-float g p q))
		(if (evenp N)
		    (/ p q)
		    (/ q (- p)))))))))
|#

(defun %sp-tan-long (x)
  (declare (long-float x))
  (let ((y (abs x)))
    (declare (long-float y))
    (if (< %long-tan-max-x y)
	(error "~s not in domain." 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)))))))))

;tan of a complex number
;there was a nicer algorithm but it turned out not to work so well.
(defun %sp-tan-complex (z)
  (let* ((num (sin z))
	 (denom (cos z)))
    (if (zerop denom) (error "~S undefined tangent." z)
	(/ num denom))))

@