(defun phi (x)
  "Adopted from CLASP 1.4.3, see copyright notice at http://eksl-www.cs.umass.edu/clasp.html"
 ;(test-variables (x number))
  (setf x (coerce x 'double-float))
  (locally (declare (type double-float x))
    (* 0.5d0 (+ 1.0d0 (error-function (/ x (sqrt 2.0d0)))))))

(defun sign (x)
  (cond ((minusp x) -1)
	((plusp x) 1)
	((zerop x) 0)
	(t nil)))



(defun gamma-incomplete (a x)
  "Adopted from CLASP 1.4.3, <a href=\"http://eksl-www.cs.umass.edu/clasp.html\">http://eksl-www.cs.umass.edu/clasp.html</a>"
  (declare (optimize (safety 3)))
  (setq a (coerce a 'double-float))
  (let ((gln (the double-float (gamma-ln a))))
    (when (= x 0.0)
      (return-from gamma-incomplete (values 0.0d0 gln)))
    (if (< x (+ a 1.0d0))
        ;; Use series representation.  The following is the code of what
        ;; Numerical Recipes in C calls ``GSER'
        (let* ((itmax 1000)
               (eps   3.0d-7)
               (ap    a)
               (sum   (/ 1d0 a))
               (del sum))
          (declare (type double-float ap sum del))
          (dotimes (i itmax)
            (incf ap 1.0d0)
            (setf del (* del (/ x ap)))
            (incf sum del)
            (if (< (abs del) (* eps (abs sum)))
                (let ((result (underflow-goes-to-zero
                               (* sum (safe-exp (- (* a (log x)) x gln))))))
                  (return-from gamma-incomplete (values result gln)))))
          (error "Series didn't converge:~%~
                  Either a=~s is too large, or ITMAX=~d is too small." a itmax))
        ;; Use the continued fraction representation.  The following is the
        ;; code of what Numerical Recipes in C calls ``GCF.'' Their code
        ;; computes the complement of the desired result, so we subtract from
        ;; 1.0 at the end.
        (let ((itmax 1000)
              (eps   3.0e-7)
              (gold 0d0) (g 0d0) (fac 1d0) (b1 1d0) (b0 0d0)
              (anf 0d0) (ana 0d0) (an 0d0) (a1 x) (a0 1d0))
          (declare (type double-float gold g fac b1 b0 anf ana an a1 a0))
          (dotimes (i itmax)
            (setf an  (float (1+ i))
                  ana (- an a)
                  a0  (* fac (+ a1 (* a0 ana)))
                  b0  (* fac (+ b1 (* b0 ana)))
                  anf (* fac an)
                  a1  (+ (* x a0) (* anf a1))
                  b1  (+ (* x b0) (* anf b1)))
            (unless (zerop a1)
              (setf fac (/ 1.0d0 a1)
                    g   (* b1 fac))
              (if (< (abs (/ (- g gold) g)) eps)
                  (let ((result (underflow-goes-to-zero
                                 (* (safe-exp (- (* a (log x)) x gln)) g))))
                    (return-from
                     gamma-incomplete (values (- 1.0 result) gln)))
                  (setf gold g))))
          (error "Continued Fraction didn't converge:~%~
                  Either a=~s is too large, or ITMAX=~d is too small." a
                  itmax)))))


(defmacro square (x)
  `(* ,x ,x))

(defun error-function (x)
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  ;; per CMUCL type-checker, the results of this function may not be
  ;; correct, if the input isn't a double-float.  For now, I think
  ;; it's easiest to coerce, but later it would be better to ensure
  ;; that the callers do the right thing."
  (let ((erf (gamma-incomplete 0.5d0 
			       (square (coerce x 'double-float)))))
    (if (>= x 0.0d0) erf (- erf))))


(defun wilcoxon-signed-rank-test (differences &optional (tails :both))
  (let* ((nonzero-differences (remove 0 differences :test #'=))
         (sorted-list (sort (mapcar #'(lambda (dif)
                                        (list (abs dif) (sign dif)))
                                    nonzero-differences)
                            #'<
                            :key #'first))
         (distinct-values (delete-duplicates (mapcar #'first sorted-list)))
         (ties nil))

    (when (< (length nonzero-differences) 16)
      (error "This Wilcoxon Signed-Rank Test (normal approximation method) requires nonzero N > 15"))

    (unless (member tails '(:positive :negative :both))
      (error "tails must be one of :positive, :negative or :both, not ~a" tails))
    
    ; add avg-rank to the sorted values
    
    (dolist (value distinct-values)
      (let ((first (position value sorted-list :key #'first))
            (last (position value sorted-list :key #'first :from-end t)))
        (if (= first last)
            (nconc (find value sorted-list :key #'first) (list (1+ first)))
            (let ((number-tied (1+ (- last first)))
                  (avg-rank (1+ (/ (+ first last) 2)))) ; +1 since 0 based
              (push number-tied ties)
              (dotimes (i number-tied)
                (nconc (nth (+ first i) sorted-list) (list avg-rank)))))))
    (setq ties (nreverse ties))
    (let* ((direction (if (eq tails :negative) -1 1))
           (r1 (reduce #'+ (mapcar #'(lambda (entry)
                                       (if (= (second entry) direction)
                                           (third entry)
                                           0))
                                   sorted-list)))
           (n (length nonzero-differences))
           (expected-r1 (/ (* n (1+ n)) 4))
           (ties-factor (if ties
                            (/ (reduce #'+ (mapcar #'(lambda (ti)
                                                       (- (* ti ti ti) ti))
                                                   ties))
                               48)
                            0))
           (var-r1 (- (/ (* n (1+ n) (1+ (* 2 n))) 24) ties-factor))
           (T-score (/ (- (abs (- r1 expected-r1)) 1/2) (sqrt var-r1))))
      (* (if (eq tails :both) 2 1) (- 1 (phi T-score))))))