;; Numerical Methods in Emacs
;; Author: Tiago Charters de Azevedo (tca) <tca@diale.org>
;; 2008-03-03
;; License - GPL v. 3


;; Basic functions

(defun pow (x n)
  (exp (* (log x) n)))

(defun mapcarpow (x n)
  (let* ((xaux (reverse x)) (xp nil))
    (while xaux
      (setq xp (cons (exp (* (log (car xaux)) n)) xp))
      (setq xaux (cdr xaux)))
    xp))

(defun mapcarf (fun x)
  "Apply FUN to each element of X, and make a list of the results.
Syntax example: (mapcarf quote sin quote (1 2))"
  (mapcar fun x))


(defun mapcarexp (x)
  (mapcar 'exp x)) 

(defun mapcarlog (x)
  (mapcar 'log x))

(defun mapcarsin (x)
  (mapcar 'sin x))

(defun mapcarcos (x)
  (mapcar 'cos x))

(defun mapcartan (x)
  (mapcar 'tan x))

(defun mapcarabs (x)
  (mapcar 'abs x))


(defun polyval2 (x coef)
  (+ (car coef) (* x (cadr coef))))

(defun polyval (x coef)
  (let ((coefx coef) (px 0))
    (while coefx
      (setq px (+ (* x px)  (car coefx)))
      (setq coefx (cdr coefx)))
    px))

(defun polyderiv (coef)
  "Return the coefficients of the derivative of the polynomial whose
     coefficients are given by vector COEF"
  (cdr (prod  (linspace 0 (- (length coef) 1) (- (length coef) 1)) coef)))

;; Linear algebra functions

(defun  linspace (xmin xmax n)
  (let ((i 0) (lv nil))
    (while (<= i n)
      (setq lv (cons (+ xmin (* (* i (exp (-(log n)))) (- xmax xmin))) lv))
      (setq i (+ 1 i)))
    (reverse lv)))

(defun prod (x y)
  (let ((i 0) (z nil))
    (while (< i (length x))
      (setq z (cons (* (nth i x) (nth i y)) z))
      (setq i (+ i 1)))
    (reverse z)))


(defun inner (x y)
  "Inner product of two vectors X and Y."
  (let ((i 0) (aux 0))
    (while (< i (length x))
      (setq aux  (+ aux (* (nth i x) (nth i y))))
      (setq i (1+  i)))
    aux))

(defun sum (x)
  "Sum of element of X."
  (let ((i 0) (aux 0))
    (while (< i (length x))
      (setq aux (+ aux (nth i x)))
      (setq i (1+ i)))
    aux))

(defun sumsq (x)
  "Sum of squares of elements X."
  (let ((i 0) (aux 0))
    (while (< i (length x))
      (setq aux (+ aux (* (nth i x) (nth i x))))
      (setq i (1+ i)))
    aux))


;; Statistics 

(defun mean (x)
  "If X is a vector, compute the arithmetic mean of the elements of X

          mean (x) = SUM_i x(i) / N"
(* (sum x)(exp (-(log (length x))))))

(defun mapcarmean (x)
  "Compute the arithmetic mean for each row and return them
     in a list.")

;; Non-linear functions

(defun bisection (f a b &optional tol)
 "Solve the function f(x)=0 with bounds a and b and tolerance tol. a
and b need to bracket the solution.
Example:
 - define the function
(defun f (x) (- 2 (* x x)))
 - pply the method
(bisection 'f 1 3 .000001);1.5707969665527344
"
 (when (null tol)
   (setq tol 0.00001))
 (let* 
     ((m (/ (+ a b) 2.0))
      (fm (funcall f m))
      (fa (funcall f a))
      )
   (if (< (abs fm) tol)
        m
     (if (> (* fa fm) 0.0)
          (bisection 'f m b tol)
        (bisection 'f a m tol)))))

(defun fixedpoint (f x &optional tol)
  "Find the fixed point of f, i. e., f(p)=p.
Example:
- define the function
(defun f (x) (cos x))
- apply the method
(fixedpoint 'f 1 .001);0.7387603198742113
"
 (when (null tol)
   (setq tol 0.00001))
 (let* 
     ((oldx x)
      (x (funcall f x))
      )
   (if (< (abs (- x oldx)) tol) x (fixedpoint 'f x tol))
   ))