;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; Includes information derived from the GSL (Gnu Scientific Library)
;;;   header files which are copyrighted and distributed under the GPL
;;;   by their respective authors.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(libload "gsl/gsl-config")
(libload "gsl/gsl-idx")
(libload "gsl/complex")

;; dummy function that adds GSL C header files in C file
(eval `(de dummy_libcomplex () ,gsl-cpheader-cmd ()))

#? *** complex: Complex Numbers Library
;; This is a library of functions on complex numbers.
;; a complex number, as manipulated by these functions,
;; is simply represented as an idx1 of double with 2 elements.
;; Creating a complex literal can be done with {<c> (setq z [3 4])},
;; just like with vectors literal.
;; Allocating an empty complex is done with {<c> (setq z (new-complex))},
;; which is equivalent to {<c> (setq z (double-matrix 2))}.

#? (new-complex)
;; returns a new complex (i.e. an idx1 of double of size 2).
(dmd new-complex () `(double-matrix 2))

#? (complex-matrix [<s0> [<s1> ...<sn>]])
;; returns a new complex matrix of size <s0>,<s1>,...<sn>.
;; A complex matrix is nothing more than a double matrix
;; with an extra dimension of size two. In other words,
;; {<c> (complex-matrix)} returns a double matrix of size 2
;; (equivalent to {<c> (double-matrix 2)} and to {<c> (new-complex)}.
;; Similarly, {<c> (complex-matrix 3 3 4)} returns a double matrix of size 3x3x4x2,
;; and so on.
(dmd complex-matrix p (cons 'double-matrix (append p '(2))))

#? ** Complex functions that return a real

#? (carg <z>)
;; argument of <z> (returns angle between -pi and pi).
#? (cabs <z>)
;; magnitude/modulus of <z>
#? (cabs2 <z>)
;; sqaure magnitude/square modulus of <z>
#? (clogabs<z>)
;; logarithm of magnitude of <z>.

#? ** Dyadic Complex Operators
;; these functions allocate and return a complex result.
;; example:
;; {<code>
;; ? (c* [1 2] [3 5])
;; = [-7.00 11.00 ]
;; </code>}

#? (c+ <a> <b>)
;; addition
#? (c- <a> <b>)
;; subtraction
#? (c* <a> <b>)
;; product
#? (c/ <a> <b>)
;; division
#? (c** <a> <b>)
;; power
#? (clogb <a> <b>)
;; compute log(<a>)/log(<b>), logarithms of <a> in complex base <b>.

#? ** Unary Complex Operators
;; complex functions that return a complex.
;; example:
;; {<code>
;;  ? (csqrt [0 2])
;;  = [ 1.00  1.00 ]
;; </code>}

#? (cconjugate <z>)
;; complex conjugate
#? (cinverse <z>)
;; complex inverse
#? (cnegative <z>)
;; negative
#? (csqrt <z>)
;; square root
#? (cexp <z>)
;; exponential
#? (clog <z>)
;; log
#? (clog10 <z>)
;; log base 10
#? (csin <z>)
;; sine
#? (ccos <z>)
;; cosine
#? (ctan <z>)
;; tangent
#? (csec <z>)
;; secant
#? (ccsc <z>)
;; cosecant
#? (ccot <z>)
;; cotangent
#? (carcsin <z>)
;; arcsine
#? (carccos <z>)
;; arccosine
#? (carctan <z>)
;; arctangent
#? (csinh <z>)
;; hyperbolic sine
#? (ccosh <z>)
;; hyperbolic scosine
#? (ctanh <z>)
;; hyperbolic tangent
#? (csech <z>)
;; hyperbolic secant
#? (ccsch <z>)
;; hyperbolic cosecant
#? (ccoth <z>)
;; hyperbolic cotangent
#? (carcsinh <z>)
;; hyperbolic arcsine
#? (carccosh <z>)
;; hyperbolic arccos
#? (carctanh <z>)
;; hyperbolic arctangent
#? (carcsech <z>)
;; hyperbolic arcsecant


(dm def-complex-op2 (s op func)
  `(de ,(eval op) (a b)
     ((-idx1- (-double-)) a b)
     (let ((r (double-matrix 2)))
       (,(eval func) (idx-ptr a) (idx-ptr b) (idx-ptr r)) r)))

(dm def-complex-op1 (s op func)
  `(de ,(eval op) (a)
     ((-idx1- (-double-)) a)
     (let ((r (double-matrix 2)))
       (,(eval func) (idx-ptr a) (idx-ptr r)) r)))

(each ((z '((c+ gsl_complex_add)
	    (c- gsl_complex_sub)
	    (c* gsl_complex_mul)
	    (c/ gsl_complex_div)
	    (c** gsl_complex_pow)
	    (clogb gsl_complex_log_b))))
  (def-complex-op2 (car z) (cadr z)))

;; These are predefined in latest versions of lush.
(unlock-symbol carg cabs)

(de carg (z)
  ((-idx1- (-double-)) z)
  (double #{ gsl_complex_arg (*IDX_PTR($z,gsl_complex)) #} ))
(de cabs (z)
  ((-idx1- (-double-)) z)
  (double #{ gsl_complex_abs (*IDX_PTR($z,gsl_complex)) #} ))
(de cabs2 (z)
  ((-idx1- (-double-)) z)
  (double #{ gsl_complex_abs2 (*IDX_PTR($z,gsl_complex)) #} ))
(de clogabs (z)
  ((-idx1- (-double-)) z)
  (double #{ gsl_complex_logabs (*IDX_PTR($z,gsl_complex)) #} ))

(each ((z '((cconjugate gsl_complex_conjugate)
	    (cinverse gsl_complex_inverse)
	    (cnegative gsl_complex_negative)
	    (csqrt gsl_complex_sqrt)
	    (cexp gsl_complex_exp)
	    (clog gsl_complex_log)
	    (clog10 gsl_complex_log10)
	    (csin gsl_complex_sin)
	    (ccos gsl_complex_cos)
	    (ctan gsl_complex_tan)
	    (csec gsl_complex_sec)
	    (ccsc gsl_complex_csc)
	    (ccot gsl_complex_cot)
	    (carcsin gsl_complex_arcsin)
	    (carccos gsl_complex_arccos)
	    (carctan gsl_complex_arctan)
	    (csinh gsl_complex_sinh)
	    (ccosh gsl_complex_cosh)
	    (ctanh gsl_complex_tanh)
	    (csech gsl_complex_sech)
	    (ccsch gsl_complex_csch)
	    (ccoth gsl_complex_coth)
	    (carcsinh gsl_complex_arcsinh)
	    (carccosh gsl_complex_arccosh)
	    (carctanh gsl_complex_arctanh)
	    (carcsech gsl_complex_arcsech)
	    (carccsch gsl_complex_arccsch)
	    (carccoth gsl_complex_arccoth)
	    )))
  (def-complex-op1 (car z) (cadr z)))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; complex vectors and matrices.
;; under construction

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(dhc-make ()
 dummy_libcomplex
 c+ c- c* c/ c** clogb
 carg cabs cabs2 clogabs

 cconjugate cinverse cnegative csqrt cexp clog clog10 csin ccos ctan
 csec ccsc ccot carcsin carccos carctan csinh ccosh ctanh csech ccsch
 ccoth carcsinh carccosh carctanh carcsech
)
