; test of math functions
(setf eps 1.e-7 x  (list 1 2 3 4 -1 -2)	y (list -4.5 6.03 -7.7 2.00 3))

(check #'= (abs 4.7) (abs 4.7))
(check #'= (abs -6.3) (abs 6.3))
(check #'= (floor x) (ceiling x))
(check #'= (- (ceiling (+ x eps)) (floor (+ x eps))) 1)
(check #'= (truncate (* x (+ 1 eps))) x)
(check #'= (cumsum (combine (select x 0) (difference x))) x)
(check #'= (min x y (uniform-rand 77)) (min (min x) (min y)))
(check #'= (max x y (uniform-rand 77)) (max (max x) (max y)))
(check #'= (pmax (iseq 1 6) (iseq 0 5) (iseq -3 2)) (iseq 1 6))
(check #'= 
       (pmin (+ (iseq 1 6) eps) (+ (iseq 0 5) eps) (+ (iseq -3 2) eps)) 
       (+ (iseq -3 2) eps))
(check #'= (median y) (quantile y .5))
(check #'< (abs (- (sum y) (first (last (cumsum y))))) eps)
(check #'= (prod x) (apply #'* x))

(check #'< (abs (- (^ 3 (log x 3)) x)) eps)
(check #'< (abs (- (^ 10 (log y 10)) y)) eps)
(check #'< (abs (- (sqrt y) (^ 10 (* (/ 1 2) (log y 10))))) eps)
(check #'< (abs (- (exp 1) (/ 1 (exp -1)))) eps)
(check #'< (abs (- (cumsum (log (iseq 1 10))) (log-gamma (iseq 2 11)))) eps)

(labels ((cis (theta) (exp (* #c(0 1) theta)))
         (norm (z) (sqrt (/ (realpart (sum (* z (conjugate z)))) (length z))))
         (texp (z) 
	   (combine (- (log (abs (exp z))) (realpart z))
		    (- (phase (exp z)) (phase (cis (imagpart z))))))
         (tlog (z) 
           (combine (- (log (abs z)) (realpart (log z)))
		    (- (phase z) (imagpart (log z)))))
         (tpow (z w) (- (^ z w) (exp (* w (log z)))))
         (reldiff (x y) (/ (- x y) (pmax 1 (abs x) (abs y))))
         (tsin (z)
           (reldiff (sin z) (/ (- (exp (* #c(0 1) z))
				  (exp (* #c(0 -1) z))) #c(0 2))))
         (tcos (z)
           (reldiff (cos z) (/ (+ (exp (* #c(0 1) z))
				  (exp (* #c(0 -1) z))) 2)))
         (ttan (z) (- (tan z) (/ (sin z) (cos z))))
         (tasin (z)
           (+ (asin z) 
	      (* #c(0 1) (log (+ (* #c(0 1) z) (sqrt (- 1 (* z z))))))))
         (tacos (z) 
           (+ (acos z)
	      (* #c(0 1) (log (+ z (* #c(0 1) (sqrt (- 1 (* z z)))))))))
         (tatan (z)
	   (+ (atan z)
	      (* #c(0 1) 
		 (log (* (+ 1 (* #c(0 1) z)) (sqrt (/ 1 (+ 1 (* z z))))))))))
  (let* ((n 500)
         (z1 (+ (/ (normal-rand n) (uniform-rand n))
		(* #c(0 1) (/ (normal-rand n) (uniform-rand n)))))
         (z2 (+ (- (* (uniform-rand n) 100) 50)
		(* #c(0 1) (- (* (uniform-rand n) 100) 50))))
         (z3 (+ (normal-rand n) (* #c(0 1) (normal-rand n))))
         (w3 (+ (normal-rand n) (* #c(0 1) (normal-rand n))))
         (r (tan (* (uniform-rand n) (/ pi 2))))
         (e2 (exp z2))
         (l1 (log z1))
         (b1 (list (+ 1 r) (- -1 r)))
         (b2 (list -1 1))
         (eps 1e-9))
    (check #'eq
	   (combine (= (phase 0) 0) (= (norm (- (phase (- r)) pi)) 0)
		    (>= (phase z1) (- pi)) (<= (phase z1) pi)
		    (< (norm (- (* (abs z1) (cis (phase z1))) z1)) eps))
	   t)
    (check #'< (norm (texp z2)) eps)
    (check #'< (norm (tpow z3 w3)) eps) 
    (check #'< (combine (norm (tsin z2)) (norm (tcos z2)) (norm (ttan z2)))
	   eps)
    (check #'< (combine (norm (tsin z3)) (norm (tcos z3)) (norm (ttan z3))) 
	   eps)
    (check #'< (combine (norm (tasin z1)) (norm (tacos z1)) (norm (tatan z1)))
	   eps)
    (check #'< (combine (norm (tasin b1)) (norm (tacos b1)) (norm (tatan b1)))
	   eps)
    (check #'<  (combine (norm (tasin b2)) (norm (tacos b2))) eps)))