;;;;;;;;;;;;;;;;;;;;;; ftensor.lisp ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                                                               ;
;    for Franz Lisp Opus 38.92 on 4.3BSD Quasijarus0c           ;
;    numerical tensor operation                                 ;
;    (c) Kenrou Adachi, 2007/02/10                              ;
;                                                               ; 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;; Utility (from the book "ANSI Common Lisp" by Paul Graham) ;;

(defun copylist (lst)
  (if (atom lst)
      lst
    (cons (car lst) (copylist (cdr lst)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;; type of tensor ;;;

(defun rectyp (typ tmp)
  (progn
    (setf tmp (car tmp))
    (if (listp tmp)
      (progn
        (setf typ (cons (length tmp) typ))
        (rectyp typ tmp))
      typ)))


(defun tensortype (a)
  (let* ((ttmp a) (ttyp (cons (length ttmp) nil)))
    (reverse (rectyp ttyp ttmp))))


;;; rank of tensor ;;;

(defun tensorrank (f)
  (length (tensortype f)))


;;; make a zero tensor ;;;

(defun zerolist (m)
  (if (= m 0)
      nil
    (cons 0 (zerolist (1- m)))))

(defun rec (m rz)
  (if (= m 0)
      nil
    (cons rz (rec (1- m) rz))))

(defun zerotensor (ttype)
  (let*
      ((n (length ttype))
       (rtyp (reverse ttype))
       (rlst (car rtyp))
       (zerot nil))

    (do ((i 1 (1+ i))
	 (gstop n))
	((> i gstop))
	(if (= i 1)
	    (progn
	      (setf zerot (zerolist rlst))
	      (setf rtyp (cdr rtyp)))
	  (progn
	    (setf rlst (car rtyp))
	    (setf zerot (rec rlst zerot))
	    (setf rtyp (cdr rtyp)))))
    
    zerot))
            
;;; tensor product ;;;

(defun sp (x) (* k x))

(defun mt (x) (mtimes k x))

(defun mtimes (k b)
  (if (not (listp (car b)))
    (mapcar #'sp b)
    (mapcar #'mt b)))


(defun sproduct (k b)
  (mtimes k b))


(defun mp (x) (mtimes x b))

(defun tp (x) (tproduct x b))


(defun tproduct (a b)
  (if (not (listp (car a)))
    (mapcar #'mp a)
    (mapcar #'tp a)))


;;; exchange tensor type ;;;

(defun endcons (x a)
  (reverse (cons x (reverse a))))

;;;

(defun exchinext (a)
  (let
    ((ni  (length a))
     (nip (length (car a)))
     (b nil))

    (do ((j 1 (1+ j))
	 (gstop nip))
	((> j gstop))
	(let ((c nil))
	  (do ((k 1 (1+ k))
	       (gstop2 ni))
	      ((> k gstop2))
	      (setf c (endcons (nth (1- j) (nth (1- k) a)) c)))
	  (setf b (endcons c b))))
    b))


(defun exchnext (a i)
  (let ((n (tensorrank a)))
    (if (or (> i (1- n)) (< i 1))
      nil
      (if (= i 1)
        (exchinext a)
        (mapcar #'(lambda (x) (exchnext x (1- i))) a)))))


(defun exchttype (a i j)
  (if (= i j)
      nil
    (if (> i j)
        (exchttype a j i)
      (if (= i (1- j))
          (exchnext a i)
        (let ((b (copylist a)))
	  (do ((k i (1+ k))
	       (gstop (1- j)))
	      ((> k gstop))
	      (setf b (exchnext b k)))
	  (do ((k i (1+ k))
	       (gstop2 (- j 2)))
	      ((> k gstop2))
	      (setf b (exchnext b (- (+ (- j 2) i) k))))
          b)))))


;;; contraction ;;;


(defun sqmatp (a)
  (if (not (numberp (car (car a))))
      nil
    (if (not (= (length a) (length (car a))))
        nil
      t)))


(defun mattrace (a)
  (if (not (sqmatp a))
      nil
    (let ((n (length a)) (tra 0))
      (do ((i 1 (1+ i))
	   (gstop n))
	  ((> i gstop))
	  (setf tra (+ tra (nth (1- i) (nth (1- i) a)))))
      tra)))


(defun cntrlst2 (a)
  (if (sqmatp a)
      (mattrace a)
    (mapcar #'(lambda (x) (cntrlst2 x)) a)))


(defun contract (a i j)
  (let ((n (tensorrank a)))
    (if (or (= i j) (< i 1) (> i n) (< j 1) (> j n))
        nil
      (if (> i j)
          (contract a j i)
        (if (= i (1- n))
            (cntrlst2 a)
          (let ((b (copylist a)))
            (if (not (= j n))
		(do ((p j (1+ p))
		     (gstop (1- n)))
		    ((> p gstop))
		    (setf b (exchnext b p)))
              nil)
            (if (not (= i (1- n)))
		(do ((p i (1+ p))
		     (gstop2 (- n 2)))
		    ((> p gstop2))
		    (setf b (exchnext b p)))
              nil)
            (setf b (cntrlst2 b))
            b))))))


;;; plus function of two tensors ;;;


(defun tpls (a b)
  (let
      ((ta (copylist a))
       (tb (copylist b)))
    (if (numberp (car ta))
        (mapcar #'+ ta tb)
      (mapcar #'tpls ta tb))))


(defun tplus (a b)
  (if (not (equal (tensortype a) (tensortype b)))
      nil
    (tpls a b)))

;;; end of ftensor.lisp ;;;