;;; Lisplab, store-operators.lisp
;;; Double float and complex double float operators (such as +,-,*, 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.

;;; TODO: 
;;; - change name of the file to something about blas store
;;; - make the file completely independent of the rest of lisplab
;;;
;;; 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.

(in-package :lisplab)

;;; TODO: there must be some easier way to generate the code in this file,
;;;       but I have not the energy to do it. I do, however, think that
;;;       the basic idea of having a layer of ordinary functions is a good one. 

;;; The reason for generating ordinary functions and not using methods,
;;; is that the real and complex stores have the same type. The fortran-compatible
;;; complex arrays are just subsequent real and complex double-floats.

;;; The reason for having both real and complex in the same file is that 
;;; not all operators function on both real and complex arguments. Care must 
;;; be taken. This is also the reason why it's hard to generate more code 
;;; automatically. 

;;; The code below generates ordinary lisp functions 
;;; for elementwise operations on simple double-float arrays. 
;;; They use a naming conventions, which should be pretty easy to 
;;; guess, such as df = double float and cdfa = complex double float array. 
;;;
;;; (The convention for complex should for consistent naming be changed to zdfa, 
;;; but its not really important)
;;; 
;;; I use map-into when its performance is equal or better to the iterations. 
;;; The iterative version for all operations are still in the file, since other lisps  
;;; than sbcl might have a slower map-into, so that they can be needed later.
;;; For real numbers, map-into can be used for all operations, while for complex 
;;; number only + and - (*, / and expt mix the real and complex parts)

;;; The below operations are based on map-into. It is hope that some 
;;; clever lisp machine can be very fast on them (On my SBCL 32 they are exactly 
;;; the same speed as the iterative versions

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun +_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'+ a b))

(defun +_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (+ x b)) a))

(defun +_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (+ a x)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun -_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'- a b))

(defun -_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (- x b)) a))

(defun -_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (- a x)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun *_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'* a b))

(defun *_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (* x b)) a))

(defun *_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (* a x)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun /_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'/ a b))

(defun /_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (/ x b)) a))

(defun /_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (/ a x)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun ^_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'expt a b))

(defun ^_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (expt x b)) a))

(defun ^_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (expt a x)) b))

(defun ^_dfa>=0-dfa (a b c)
  (declare (type (simple-array (double-float 0d0) (*)) a b c))
  (map-into c #'expt a b))

(defun ^_df>=0-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type (double-float 0d0) a))
  (map-into c (lambda (x) (expt a x)) b))

(defun ^_dfa>=0-df (a b c)
  (declare (type (simple-array (double-float 0d0) (*)) a c)
	   (type double-float b))
  (map-into c (lambda (x) (expt x b)) a))

(defun ^_dfa-+2 (a c)
  (declare (type type-blas-store a c))
  (dotimes (i (length a))
    (setf (aref c i) (expt (aref a i) 2)))
  c)

(defun ^_dfa-+3 (a c)
  (declare (type type-blas-store a c))
  (dotimes (i (length a))
    (setf (aref c i) (expt (aref a i) 3)))
  c)

(defun ^_dfa--1 (a c)
  (declare (type type-blas-store a c))
  (dotimes (i (length a))
    (setf (aref c i) (expt (aref a i) -1)))
  c)

(defun ^_dfa--2 (a c)
  (declare (type type-blas-store a c))
  (dotimes (i (length a))
    (setf (aref c i) (expt (aref a i) -2)))
  c)

(defun ^_dfa--3 (a c)
  (declare (type type-blas-store a c))
  (dotimes (i (length a))
    (setf (aref c i) (expt (aref a i) -3)))
  c)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun max_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'max a b))

(defun max_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (max x b)) a))

(defun max_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (max a x)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun min_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (map-into c #'min a b))

(defun min_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (x) (min x b)) a))

(defun min_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (x) (min a x)) b))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun =_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (length a))
    (setf (aref c i) (if (= (aref a i)
			    (aref b i))
			 1.0
			 0.0))))
  
(defun =_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (a) (if (= a b) 1.0 0.0)) a))

(defun =_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (b) (if (= a b) 1.0 0.0)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun /=_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (length a))
    (setf (aref c i) (if (/= (aref a i)
			    (aref b i))
			 1.0
			 0.0))))
  
(defun /=_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (a) (if (/= a b) 1.0 0.0)) a))

(defun /=_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (b) (if (/= a b) 1.0 0.0)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun >_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (length a))
    (setf (aref c i) (if (> (aref a i)
			    (aref b i))
			 1.0
			 0.0))))
  
(defun >_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (a) (if (> a b) 1.0 0.0)) a))

(defun >_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (b) (if (> a b) 1.0 0.0)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun >=_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (length a))
    (setf (aref c i) (if (>= (aref a i)
			    (aref b i))
			 1.0
			 0.0))))
  
(defun >=_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (a) (if (>= a b) 1.0 0.0)) a))

(defun >=_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (b) (if (>= a b) 1.0 0.0)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun <_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (length a))
    (setf (aref c i) (if (< (aref a i)
			    (aref b i))
			 1.0
			 0.0))))
  
(defun <_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (a) (if (< a b) 1.0 0.0)) a))

(defun <_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (b) (if (< a b) 1.0 0.0)) b))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun <=_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (length a))
    (setf (aref c i) (if (<= (aref a i)
			    (aref b i))
			 1.0
			 0.0))))
  
(defun <=_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (map-into c (lambda (a) (if (<= a b) 1.0 0.0)) a))

(defun <=_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (map-into c (lambda (b) (if (<= a b) 1.0 0.0)) b))


;;; Complex map-into-operations

(defun +_cdfa-cdfa (a b c)
  (+_dfa-dfa a b c))

(defun -_cdfa-cdfa (a b c)
  (-_dfa-dfa a b c))




;;; Now the complex operators


;;;  Array and number

(define-constant +operators-cdfa-cdf-map+ 
    '((+ . +_cdfa-cdf) ; iterative version
      (- . -_cdfa-cdf) ; iterative version
      (* . *_cdfa-cdf)
      (/ . /_cdfa-cdf)
      (expt . ^_cdfa-cdf) 
      ))
	   
(defmacro defun-cdfa-cdf (name op)
  (let ((a (gensym))
	(b (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,b ,out) 
       (declare (type type-blas-store ,a ,out)
		(type (complex double-float) ,b))
       (dotimes (,i (floor (length ,a) 2))
	 (setf (vref-blas-complex-store ,out ,i) 
	       (coerce (,op (vref-blas-complex-store ,a ,i) ,b) '(complex double-float))))
       ,out)))

(defmacro expand-operators-cdfa-cdf-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-cdfa-cdf ,(cdr x) ,(car x)))
	      +operators-cdfa-cdf-map+)))

(expand-operators-cdfa-cdf-map)

;;; The three parts of code below contains some common strucutre that could 
;;; in principle be joined, and there is also some clumsy code, 
;;; but I have not the energy to find out how to clean up.

;;; Number and array

(define-constant +operators-cdf-cdfa-map+ 
    '((+ . +_cdf-cdfa)
      (- . -_cdf-cdfa)
      (* . *_cdf-cdfa)
      (/ . /_cdf-cdfa)
      (expt . ^_cdf-cdfa)))

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

(defmacro expand-operators-cdf-cdfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-cdf-cdfa ,(cdr x) ,(car x)))
	      +operators-cdf-cdfa-map+)))

(expand-operators-cdf-cdfa-map)

;;; Array and array

(define-constant +operators-cdfa-cdfa-map+ 
    '(; (+ . +_cdfa-cdfa)
      ; (- . -_cdfa-cdfa)
      (* . *_cdfa-cdfa)
      (/ . /_cdfa-cdfa)
      (expt . ^_cdfa-cdfa)))

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

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

(expand-operators-cdfa-cdfa-map)


;;;; Now, some special cases of real and imaginary mixing
;;; Other cases could be optimized too, but these cases are not so obvious.

(defun complex_dfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (dotimes (i (the type-blas-idx (length a)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (aref a i)
	    (aref c 2i+1) b)))
  c)

(defun complex_df-dfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (dotimes (i (the type-blas-idx (length b)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) a
	    (aref c 2i+1) (aref b i))))
  c)

(defun complex_dfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (the type-blas-idx (length a)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (aref a i)
	    (aref c 2i+1) (aref b i))))
  c)

(defun +_cdfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (dotimes (i (floor (the type-blas-idx (length a)) 2))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (+ (aref a 2i) b)
	    (aref c 2i+1) (aref a 2i+1) )))
  c)

(defun +_cdfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (the type-blas-idx (length b)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (+ (aref a 2i) (aref b i))
	    (aref c 2i+1) (aref a 2i+1) )))
  c)
	
(defun -_cdfa-df (a b c)
  (declare (type type-blas-store a c)
	   (type double-float b))
  (dotimes (i (floor (the type-blas-idx (length a)) 2))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (- (aref a 2i) b)
	    (aref c 2i+1) (aref a 2i+1) )))
  c)

(defun -_cdfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (the type-blas-idx (length b)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (- (aref a 2i) (aref b i))
	    (aref c 2i+1) (aref a 2i+1) )))
  c)

(defun -_df-cdfa (a b c)
  (declare (type type-blas-store b c)
	   (type double-float a))
  (dotimes (i (floor (the type-blas-idx (length b)) 2))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (- a (aref b 2i))
	    (aref c 2i+1) (- (aref b 2i+1) ))))
  c)

(defun -_dfa-cdfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (the type-blas-idx (length a)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (- (aref a i) (aref b 2i))
	    (aref c 2i+1) (- (aref b 2i+1) ))))
  c)

(defun *_cdfa-df (a b c)
  ;; Well, the same as +_dfa-df, but length is twice
  (*_dfa-df a b c))

(defun *_cdfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (the type-blas-idx (length b)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (* (aref a 2i) (aref b i))
	    (aref c 2i+1) (* (aref a 2i+1) (aref b i)))))
  c)


(defun /_cdfa-df (a b c)
  ;; Well, the same as +_dfa-df, but length is twice
  (/_dfa-df a b c))

(defun /_cdfa-dfa (a b c)
  (declare (type type-blas-store a b c))
  (dotimes (i (the type-blas-idx (length b)))
    (let* ((2i  (truly-the type-blas-idx (* 2 i)))
	   (2i+1  (the type-blas-idx (1+ 2i))))
      (declare (type type-blas-idx 2i 2i+1))
      (setf (aref c 2i) (/ (aref a 2i) (aref b i))
	    (aref c 2i+1) (/ (aref a 2i+1) (aref b i)))))
  c)









#| Below is the iterative verions. They might be outdated
;;;  Array and number

(define-constant +operators-dfa-df-map+ 
    '((+ . +_dfa-df)
      (- . -_dfa-df)
      (* . *_dfa-df)
      (/ . /_dfa-df)
      (expt . ^_dfa-df) 
      (max . max_dfa-df) 
      (min . min_dfa-df)
      (log . log_dfa-df) ;; Note the log here
      ))
	   
(defmacro defun-dfa-df (name op)
  (let ((a (gensym))
	(b (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,b ,out) 
       (declare (type type-blas-store ,a ,out)
		(type double-float ,b))
       (dotimes (,i (length ,a))
	 (setf (aref ,out ,i) (,op (aref ,a ,i) ,b)))
       ,out)))

(defmacro expand-operators-dfa-df-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-dfa-df ,(cdr x) ,(car x)))
	      +operators-dfa-df-map+)))

#+nil (expand-operators-dfa-df-map) ; These are the iterative versions 

;;; The three parts of code below contains some common strucutre that could 
;;; in principle be joined, and there is also some clumsy code, 
;;; but I have not the energy to find out how to clean up.

;;; Number and array

(define-constant +operators-df-dfa-map+ 
    '((+ . +_df-dfa)
      (- . -_df-dfa)
      (* . *_df-dfa)
      (/ . /_df-dfa)
      (expt . ^_df-dfa) 
      (max . max_df-dfa) 
      (min . min_df-dfa)
      (log . log_df-dfa)))

(defmacro defun-df-dfa (name op)
  (let ((a (gensym))
	(b (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,b ,out) 
       (declare (type type-blas-store ,b ,out)
		(type double-float ,a))
       (dotimes (,i (length ,b))
	 (setf (aref ,out ,i) (,op ,a (aref ,b ,i))))
       ,out)))

(defmacro expand-operators-df-dfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-df-dfa ,(cdr x) ,(car x)))
	      +operators-df-dfa-map+)))

#+nil (expand-operators-df-dfa-map) ; These are the iterative versions 

;;; Array and array

(define-constant +operators-dfa-dfa-map+ 
    '((+ . +_dfa-dfa)
      (- . -_dfa-dfa)
      (* . *_dfa-dfa)
      (/ . /_dfa-dfa)
      (expt . ^_dfa-dfa) 
      (max . max_dfa-dfa) 
      (min . min_dfa-dfa)
      (log . log_dfa-dfa)))

(defmacro defun-dfa-dfa (name op)
  (let ((a (gensym))
	(b (gensym))
	(out (gensym))
	(i (gensym)))
    `(defun ,name (,a ,b ,out) 
       (declare (type type-blas-store ,a ,b ,out))
       (dotimes (,i (length ,a))
	 (setf (aref ,out ,i) (,op (aref ,a ,i) (aref ,b ,i))))
       ,out)))

(defmacro expand-operators-dfa-dfa-map ()
  (cons 'progn
      (mapcar (lambda (x)
		`(defun-dfa-dfa ,(cdr x) ,(car x)))
	      +operators-dfa-dfa-map+)))

#+nil (expand-operators-dfa-dfa-map) ; These are the iterative versions 


|#