;;; Lisplab, store-ordinary-functions.lisp
;;; Double float and complex double float ordinary functions (such as sin, log, etc) on
;;; simple arrays. Used by the matrix implementations.
;;; 

;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; 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.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

;;; This file contains manipulations of simple double-float arrays 
;;; and should be called by the spesialized matrix methods. 
;;; The purpose of this layer is that it can be used by 
;;; many classes such as matrix-base-dge and matrix-base-ddi, etc. 
;;; 
;;; The content of this file must be highly optimized 
;;; and should not depend anything exept Common Lisp itself.
;;;
;;; I use map-into for the real functions (impossible for the complex functions)
;;; but keep the iterative versions still in the file since future versions 
;;; of lisplab might use them.

;;; Generate more real-to-real functions. With some kind of input these will 
;;; fail and give complex output but for speed it can be ok to have them

(in-package :lisplab)

;;; Now the ordinary functions 

;;; Double-float to double float

(define-constant +ordinary-functions-dfa-to-dfa-map+ 
    ;; Real to real functions. Note that not all of the 
    ;; functions are safe to use, such as sqrt. The caller 
    ;; must make sure that the input give real output also.
    '((sin . sin_dfa-to-dfa) 
      (cos . cos_dfa-to-dfa) 
      (tan . tan_dfa-to-dfa)
      (asin . asin_dfa-to-dfa)
      (acos . acos_dfa-to-dfa)
      (atan . atan_dfa-to-dfa) 
      (sinh . sinh_dfa-to-dfa) 
      (cosh . cosh_dfa-to-dfa) 
      (tanh . tanh_dfa-to-dfa) 
      (asinh . asinh_dfa-to-dfa) 
      (acosh . acosh_dfa-to-dfa) 
      (atanh . atanh_dfa-to-dfa)
      (sqrt . sqrt_dfa-to-dfa)
      (exp . exp_dfa-to-dfa) 
      (log . log_dfa-to-dfa)
      (abs . abs_dfa-to-dfa)
      (signum . signum_dfa-to-dfa)))

;; the iterative version 
#+lisplab-iterative (defmacro defun-dfa-to-dfa-iterative (name op)
  (let ((a (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,out) 
       (declare (type type-blas-store ,a ,out))
       (dotimes (,i (length ,a))
	 (setf (aref ,out ,i) (,op (aref ,a ,i))))
       ,out)))

;; the iterative version 
#+lisplab-iterative (defmacro expand-ordinary-functions-dfa-to-dfa-map-iterative ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-dfa-to-dfa-iterative ,(cdr x) ,(car x)))
	      +ordinary-functions-dfa-to-dfa-map+)))

;; the iterative version 
#+lisplab-iterative (expand-ordinary-functions-dfa-to-dfa-map-iterative) 

;;; Non-iterative version
#-lisplab-iterative 
(defmacro defun-dfa-to-dfa (name op)
  (let ((a (gensym))
	(out (gensym)))
    `(defun ,name (,a ,out) 
       (declare (type type-blas-store ,a ,out))
       (map-into ,out #',op ,a)
       ,out)))

;;; Non-iterative version
#-lisplab-iterative 
(defmacro expand-ordinary-functions-dfa-to-dfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-dfa-to-dfa ,(cdr x) ,(car x)))
	      +ordinary-functions-dfa-to-dfa-map+)))

;;; Non-iterative version
#-lisplab-iterative 
(expand-ordinary-functions-dfa-to-dfa-map) ; the iterative version 

(define-constant +ordinary-functions-dfa-to-cdfa-map+ 
    ;; double float to complex double float
    ;; Hm the list should include most functions.  
    '((asin . asin_dfa-to-cdfa) 
      (acos . acos_dfa-to-cdfa) 
      (atanh . atanh_dfa-to-cdfa)
      (log . log_dfa-to-cdfa) 
      (sqrt . sqrt_dfa-to-cdfa)
      (acosh . acosh_dfa-to-cdfa)
      ))

(defmacro defun-dfa-to-cdfa (name op)
  (let ((a (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,out) 
       (declare (type type-blas-store ,a ,out))
       (dotimes (,i (length ,a))
	 (setf (vref-blas-complex-store ,out ,i) 
	       (coerce (,op (aref ,a ,i)) '(complex double-float))))
       ,out)))

(defmacro expand-ordinary-functions-dfa-to-cdfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-dfa-to-cdfa ,(cdr x) ,(car x)))
	      +ordinary-functions-dfa-to-cdfa-map+)))

(expand-ordinary-functions-dfa-to-cdfa-map)

;;; Complex double float to double float

(defun abs_cdfa (a out) 
  (declare (type type-blas-store a out))
  (dotimes (i (length out))
    (setf (aref out i) 
	  (abs (vref-blas-complex-store a i))))
  out)

(defun realpart_cdfa (a out)
  (declare (type type-blas-store a out))
  (dotimes (i (length out))
    (setf (aref out i) 
	  (aref a (truly-the type-blas-idx (* 2 i)))))
  out)

(defun imagpart_cdfa (a out)
  (declare (type type-blas-store a out))
  (dotimes (i (length out))
    (setf (aref out i) 
	  (aref a (truly-the type-blas-idx (1+ (* 2 i))))))
  out)

;;; Complex double float to complex double float

(define-constant +ordinary-functions-cdfa-to-cdfa-map+ 
    '((sin . sin_cdfa-to-cdfa) 
      (cos . cos_cdfa-to-cdfa) 
      (tan . tan_cdfa-to-cdfa) 
      (atan . atan_cdfa-to-cdfa)
      (sinh . sinh_cdfa-to-cdfa) 
      (cosh . cosh_cdfa-to-cdfa) 
      (tanh . tanh_cdfa-to-cdfa) 
      (asinh . asinh_cdfa-to-cdfa) 
      (acosh . acosh_cdfa-to-cdfa) 
      (exp . exp_cdfa-to-cdfa) 
      (log . log_cdfa-to-cdfa) 
      (sqrt . sqrt_cdfa-to-cdfa) 
      #+nil (conjugate . conjugate_cdfa) ;; separate implementation!
      (asin . asin_cdfa-to-cdfa) 
      (acos . acos_cdfa-to-cdfa) 
      (atanh . atanh_cdfa-to-cdfa)))

(defmacro defun-cdfa-to-cdfa (name op)
  (let ((a (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,out) 
       (declare (type type-blas-store ,a ,out))
       (dotimes (,i (floor (length ,a) 2))
	 (setf (vref-blas-complex-store ,out ,i) (,op (vref-blas-complex-store ,a ,i))))
       ,out)))

(defmacro expand-ordinary-functions-cdfa-to-cdfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-cdfa-to-cdfa ,(cdr x) ,(car x)))
	      +ordinary-functions-cdfa-to-cdfa-map+)))

(expand-ordinary-functions-cdfa-to-cdfa-map)

;; Conjugate

(defun conjugate_cdfa-to-cdfa (a out)
  (declare (type type-blas-store a out))
  (dotimes (i (floor (length a) 2))
    (let* ((2i (* 2 i))
	   (2i+1 (1+ 2i)))
      (declare (type type-blas-idx i 2i 2i+1))    
      (setf (aref out 2i) (aref a 2i)
	    (aref out 2i+1) (- (aref a 2i+1)))))
  out) 

;;; Copies contents from real to complex
(defun copy_dfa-to-cdfa (a out)
  (declare (type type-blas-store a out))
  (dotimes (i (length a))
    (let* ((2i (* 2 i)))
      (declare (type type-blas-idx i 2i))    
      (setf (aref out 2i) (aref a i))))
  out)