;;; -*- Mode:LISP; Readtable:ZL; Base:10 -*-

;;;This is the upated version of Arcsine to replace the original broken version -- DDB

(defun asin-fast-2 (x)
  "ARCSINE of X - optimized for speed of execution for real number results."
  ;;Test and coerce argument:
  (if (zerop x) x
    (setq x (cond
	      ((complexp x) x)
	      ((<= -1.0 x 1.0) (float x))
	      (t (complex x)))))
  ;;Dispatch on argument type:
  (if (complexp x)
      ;;;If X is complex use formal definition of arcsin:
      (* %%-i (log (+ (* %%i x) (sqrt (- 1 (* x x))))))
    ;;Otherwise use fast approximation for real numbers:
    (let* ((z (abs x))
	   (zz 1.0f0)
	   (asinz 1.5707963050))		;first coefficient for polynomial
      ;;Summation
      (do ((n 0 (+ n 1)))
	  ((> n 6))
	(incf asinz (* (aref '#(-0.2145988016	;rest of coefficient array for polynomial part
			         0.0889789874
			        -0.0501743046
			         0.0308918810
			        -0.0170881256
			         0.0066700901
			        -0.0012624911) n)
		       (setq zz (* zz z)))))
      (setq asinz (- %%pi%2 (* (sqrt (- 1.0f0 z)) asinz)))
      ;;Sign adjust
      (if (minusp x) (- asinz) asinz))))


;;;This is an attempt at beeting Sqrt's time, works best near 1.0 -- DDB

(defun rsqrt (x)
  (let* (
	 (div 1.0)
	 (chk x)
	 )
    (loop
      (setq chk (// x (setq div (// (+ chk div) 2.0))))
      (if (<= div chk)
	  (return-from rsqrt div)))))


;;;Some handy constants...

(defconstant pi 3.141592653589793238462643383279503
             "The mathematical constant PI.")

(defconstant %%pi%2 1.570796326794896619231321691639751)

(defconstant %%i #c(0.0 1.0)
	     "The square-root of -1.")

(defconstant %%-i #c(0.0 -1.0)
	     "Zero minus the square-root of -1.")

(defmacro ignore-some-arithmetic-errors (&rest forms)
  `(condition-bind (((sys:divide-by-zero) 'divide-by-zero-handler)
		    ((sys:zero-log) 'zero-log-handler))
     (catch 'my-random-arithmetic-catch-tag
       (catch 'zero-log
	 (catch 'divide-by-zero
	   (progn ,@forms))))))

(defun divide-by-zero-handler (&rest args)
  args
  (throw 'divide-by-zero t))

(defun zero-log-handler (&rest args)
  args
  (throw 'zero-log t))

(defun plot ()
  (do-forever 
    (let* ((depth-function-list (make-nth-derivative))
	  (xrange (prompt-and-read :number "what is xrange?"))
	  (xo (prompt-and-read :number "what is x coordinate?"))
	  (yo (prompt-and-read :number "what is y coordinate?"))
	  (function ())
	  (depth ())
	  (z 0.)
	  (y 0.)
	  (xx (- 475 (* xo (// 475 xrange))))
	  (yy (+ 325 (* yo (// 475 xrange))))
	  (test ()))
      (setq function (cadr depth-function-list))
      (setq depth (car depth-function-list))
      (if (y-or-n-p "clear screen before plotting?") (send *terminal-io* :refresh))
      (send *terminal-io* :draw-line
	    xx
	    0
	    xx
	    650)
      (send *terminal-io* :draw-line
	    0
	    yy
	    950
	    yy)
      (do ((x 0 (+ x 1)))
	  ((= x 950) nil)
        (ignore-some-arithmetic-errors
	  (setq z y)
	  (if (and (not (complexp z))
		   (not (complexp y))
		   (< 0 (setq y (findy xrange xo yo function x depth))  650)
		   (or (> z (setq test (findy xrange xo yo function (- x 0.5) depth)) y)
		       (< z test y)
		       (< (abs (- z y)) 10))	
		   )
	      (send my-win :draw-line
		    x
		    z
		    (+ x 1)
		    y)))))))

(defun findy (xrange xo yo function x depth)
  (let
    ((val
       (+ (* yo (// 475.0 xrange))
	   (- 325.0
	      (* (// 475.0 xrange)
		 (slope (+ xo (* xrange (// (- x 475.0) 475.0))) depth function))))))
    (if (not (complexp val))
	(floor val)
      (throw 'my-random-arithmetic-catch-tag 't))))


(defun make-nth-derivative (&optional (depth 0))
  (let ((function (if (= depth 0) (prompt-and-read :expression "Function to plot?")
		    (prompt-and-read :expression "slope of?"))))
    (if (not (eq function 'slope))
	`(,depth ,(compile () `(lambda (x y z) (,function x))))
      (make-nth-derivative (1+ depth)))))

(defun slope (x depth function)
  (cond ((= depth 0)
	 (funcall function x 0 0))
	(t (let ((f 'slope))
	     (// (- (funcall f (+ x (* 0.03 depth)) (1- depth) function)
		    (funcall f (- x (* 0.03 depth)) (1- depth) function))
		 (* 0.06 depth))))))
  
(defun cot (x)
  (// 1 (tan x)))
(defun sec (x)
  (// 1 (cos x)))
(defun csc (x)
  (// 1 (sin x)))
(defun coth (x)
   (// 1 (tanh x)))
(defun sech (x)
   (// 1 (cosh x)))
(defun csch (x)
   (// 1 (sinh x)))

(defun asin-acc (x &aux (i 0.0+1.0i))
  "ARCSINE of X - optimized for accuracy."
  (let* ((z)
	 (asinz))
    (if (zerop x) (return-from asin-acc x))
    (setq z (if (complexp x) x (complex x)))
    (setq asinz (- (* i (log (+ (* i z) (sqrt (- 1.0 (* z z))))))))
    (if (and (not (complexp x)) (<= (abs x) 1.0))
	(setq asinz (realpart asinz)))
    (values asinz)))

   
(defun my-asin (z)
  (cond
    ((= z 0) 0)
    ((= z 1) %%pi%2)
    ((= z -1) (- %%pi%2))
    ((> (realpart z) 0) (asin-evaluator z))
    (t (- (asin-evaluator (- z))))))
(defun asin-evaluator (z)
  (complex
    (atan (realpart z)
	  (realpart (* (sqrt (- 1 z))
		       (sqrt (+ z 1)))))
    (asinh (imagpart (* (conjugate (sqrt (- 1 z)))
			(sqrt (+ 1 z)))))))

(defun test-abs-my-asin (x)
  (- (abs (my-asin x)) (abs (asin-acc x))))
(defun test-real-my-asin (x)
  (- (realpart (my-asin x)) (realpart (asin-acc x))))
(defun test-imag-my-asin (x)
  (- (imagpart (my-asin x)) (imagpart (asin-acc x))))
(defun realasin-acc (x)
  (realpart (asin-acc x)))
(defun imagasin-acc (x)  
  (imagpart (asin-acc x)))
(defun realmy-asin (x)
  (realpart (my-asin x)))
(defun imagmy-asin (x)
  (imagpart (my-asin x)))
(defun absmy-asin (x)
  (abs (my-asin x))) 
(defun absasin-acc (x)
  (abs (asin-acc x)))
(defun fuckedimagmyasin (x)
  (imagpart (my-asin (+ x 0.0+10i))))
(defun fuckedimagasinacc (x)
  (imagpart (asin-acc (+ x 0.0+10i))))
(defun fuckedrealmyasin (x)
  (realpart (my-asin (+ x 0.0+10i))))
(defun fuckedrealasinacc (x)
  (realpart (asin-acc (+ x 0.0+10i))))
(defun imagscrodasinacc (x)
  (imagpart (asin-acc (complex 0.0 x))))
(defun imagscrodmyasin (x)
  (imagpart (my-asin (complex 0.0 x))))