;; SAMPLE USAGE
;; instantiating and filling a matrix
(setf foo (make-mat '(2 2)))
(setf (aref (mat-rc foo) 0 0) 2)
(setf (aref (mat-rc foo) 0 1) 1)
(setf (aref (mat-rc foo) 1 0) 0)
(setf (aref (mat-rc foo) 1 1) 1)
;; see the 2D array
(mat-rc foo)
;=> #2A((2 1) (0 1))
;; find eigenvalues
(mat-eigenvalues foo)
;=> 2.0
;  -1.0
;; find determinnt
(mat-det foo)
;=> -2
;; create another
(setf bar (make-mat '(2 2)))
(setf (aref (mat-rc bar) 0 0) 38)
(setf (aref (mat-rc bar) 0 1) -20)
(setf (aref (mat-rc bar) 1 0) 3)
(setf (aref (mat-rc bar) 1 1) 91)
(add-or-sub-mats foo bar #'-)
;=>#<a COMMON-LISP-USER::MAT>
;; subtract and add them
(mat-rc (add-or-sub-mats foo bar #'-))
;=>#2A((-36 21) (-3 -92))
(mat-rc (add-or-sub-mats foo bar #'+))
;=>#2A((40 -19) (3 90))

;; The Work...
(defclass mat ()
  ((rc :initform (make-array '(1)) ;fixme
       :type 'array
       :initarg :dimensions ;fixme
       ;:initarg :the-array
       :accessor mat-rc)))

;; override print-object to avoid hideous object printer if there's time
;; (defmethod print-object ((m mat) *standard-output*)
;;   (let ((rows (car (array-dimensions (mat-rc m))))
;; 	(cols (cadr (array-dimensions (mat-rc m)))))
;;     (do ((i 0 (incf i)))
;; 	((>= i rows))
;;       (format t "~%| ")
;;       (do ((j 0 (incf j)))
;; 	  ((>= j cols))
;; 	(format t "~A " (aref (mat-rc m) i j)))
;;       (format t " |~%"))))

(defun make-mat (dimensions)
  (if (= (length dimensions) 2)
      (let ((m (make-instance 'mat)))
	(setf (mat-rc m)
	      (make-array dimensions :initial-element 0))
	m)
      (format t "matrices are only in 2 dimensions, ")
      (format t "not ~A~%" (length dimensions))))

(defmethod add-or-sub-mats ((m1 mat) (m2 mat) op)
  (if (same-dimensions-p m1 m2)
      (let ((em (make-mat (array-dimensions (mat-rc m1))))
	    (num-rows (car (array-dimensions (mat-rc m1))))
	    (num-cols (cadr (array-dimensions (mat-rc m1)))))
	(do ((i 0 (incf i)))
	    ((>= i num-rows)
	     em)
	  (do ((j 0 (incf j)))
	      ((>= j num-cols)
	       em);probablynecessary? 
	    (setf (aref (mat-rc em) i j)
		  (funcall op
			   (aref (mat-rc m1) i j)
			   (aref (mat-rc m2) i j))))	  ))
      (format t "To add/subtract matrices, ")
      (format t "they must be of the same dimension~%")))

(defmethod scalar-mul-mat ((m1 mat) scalar)
  (let ((em (make-mat (array-dimensions (mat-rc m1))))
	(num-rows (car (array-dimensions (mat-rc m1))))
	(num-cols (cadr (array-dimensions (mat-rc m1)))))
    (do ((i 0 (incf i)))
	((>= i num-rows)
	 em)
      (do ((j 0 (incf j)))
	  ((>= j num-cols))
	(setf (aref (mat-rc em) i j)
	      (* (aref (mat-rc m1) i j) scalar))))))

(defmethod mul-mat ((m1 mat) (m2 mat))
  (if (mat-mul-allowed m1 m2)
      (let* ((m1-rows (car (array-dimensions (mat-rc m1))))
	     (m1-cols (cadr (array-dimensions (mat-rc m1))))
	     (m2-rows (car (array-dimensions (mat-rc m2))))
	     (m2-cols (cadr (array-dimensions (mat-rc m2))))
	     (em (make-mat `(,m1-rows ,m2-cols))))
	(do ((i 0 (incf i)))
	    ((>= i m1-rows))
	  (do ((j 0 (incf j)))
	      ((>= j m2-cols))
	    (do ((k 0 (incf k)))
		((>= k m1-cols))
	      (setf (aref (mat-rc em) i j)
		    (+ (aref (mat-rc em) i j)
		       (* (aref (mat-rc m1) i k)
			  (aref (mat-rc m2) k j)))))))
	em)))

(defmethod mat-det ((m mat))
  (cond ((mat-is-2-by-2 m)
	 (- (* (aref (mat-rc m) 0 0)
	       (aref (mat-rc m) 1 1))
	    (* (aref (mat-rc m) 0 1)
	       (aref (mat-rc m) 1 0))))
	((mat-is-3-by-3 m)
	 (- (+ (* (aref (mat-rc m) 0 0)
		  (aref (mat-rc m) 1 1)
		  (aref (mat-rc m) 2 2))
	       (* (aref (mat-rc m) 0 1)
		  (aref (mat-rc m) 1 2)
		  (aref (mat-rc m) 2 0))
	       (* (aref (mat-rc m) 0 2)
		  (aref (mat-rc m) 1 0)
		  (aref (mat-rc m) 2 0)))
	    (* (aref (mat-rc m) 0 2)
	       (aref (mat-rc m) 1 1)
	       (aref (mat-rc m) 2 0))
	    (* (aref (mat-rc m) 0 1)
	       (aref (mat-rc m) 1 0)
	       (aref (mat-rc m) 2 2))
	    (* (aref (mat-rc m) 0 0)
	       (aref (mat-rc m) 1 2)
	       (aref (mat-rc m) 2 1))))))

(defmethod mat-eigenvalues ((m mat))
  (if (mat-is-2-by-2 m)
      (let* ((a (aref (mat-rc m) 0 0))
	     (b (aref (mat-rc m) 0 1))
	     (c (aref (mat-rc m) 1 0))
	     (d (aref (mat-rc m) 1 1))
	     (discrim (- (square (- (- d) a)) (* 4 (- (* a d) (* b c))))))
	(if (>= discrim 0)
	    (values (/ (+ (+ a d)
			  (sqrt discrim))
		       2.0)
		    (/ (- (+ a d)
			  (sqrt discrim))
		       2.0))
	    (format t "profound mathematical reason about null space~%")))
      (format t "To find eigenvalue, matrix must be 2x2 only.~%"xs)))

(defun square (x)
  (* x x))

(defmethod mat-is-2-by-2 ((m mat))
  (and (= (car (array-dimensions (mat-rc m))) 2)
       (= (cadr (array-dimensions (mat-rc m))) 2)))

(defmethod mat-is-3-by-3 ((m mat))
  (and (= (car (array-dimensions (mat-rc m))) 3)
       (= (cadr (array-dimensions (mat-rc m))) 3)))

(defmethod same-dimensions-p ((m1 mat) (m2 mat))
  (let ((d1 (array-dimensions (mat-rc m1)))
	(d2 (array-dimensions (mat-rc m2))))
    (and (= (car d1) (car d2))
	 (= (cadr d1) (cadr d2)))))

(defmethod mat-mul-allowed ((m1 mat) (m2 mat))
  (let ((d1 (array-dimensions (mat-rc m1)))
	(d2 (array-dimensions (mat-rc m2))))
    (= (cadr d1) (car d2))))

(defmethod mat-is-square ((m mat))
  (let ((d (array-dimensions (mat-rc m))))
    (= (car d) (cadr d))))