(define (make-zero-mixer n) (make-mixer n))

(define mat
  (make-procedure-with-setter
   (lambda (m i j)
     (mixer-ref m i j))
   (lambda (m i j x)
     (mixer-set! m i j x))))

(define (mixer-diagonal? m)
  (let ((n (mus-length m)))
    (if (= n 1)
	#t
	(call-with-current-continuation
	 (lambda (return)
	   (do ((i 0 (1+ i)))
	       ((= i n) #t)
	     (do ((j 0 (1+ j)))
		 ((= j n))
	       (if (and (not (= i j))
			(not (= (mat m i j) 0.0)))
		   (return #f)))))))))
	   
(define (mixer-equal? m1 m2) 
  ;; this is less demanding than the built-in function -- it is used in snd-test.scm where we only print out 3 digits
  (let ((len (mus-length m1)))
    (and (= len (mus-length m2))
	 (call-with-current-continuation
	  (lambda (return)
	    (do ((i 0 (1+ i)))
		((= i len) #t)
	      (do ((j 0 (1+ j)))
		  ((= j len))
		(if (> (abs (- (mat m1 i j) (mat m2 i j))) .001)
		    (return #f)))))))))

(define (mixer-symmetric? m1) (mixer-equal? m1 (mixer-transpose m1)))
(define (mixer-hermitian? mx) (mixer-symmetric? mx)) ; assuming reals	   

(define (mixer-transpose mx)
  (let* ((n (mus-length mx))
	 (nmx (make-zero-mixer n)))
    (do ((i 0 (1+ i)))
	((= i n))
      (do ((j 0 (1+ j)))
	  ((= j n))
	(set! (mat nmx j i) (mat mx i j))))
    nmx))

(define (sub-matrix mx row col)
  (let* ((old-n (mus-length mx))
	 (new-n (1- old-n))
	 (nmx (make-zero-mixer new-n)))
    (do ((i 0 (1+ i))
	 (ni 0))
	((= i old-n))
      (if (not (= i row))
	  (begin
	    (do ((j 0 (1+ j))
		 (nj 0))
		((= j old-n))
	      (if (not (= j col))
		  (begin
		    (set! (mat nmx ni nj) (mat mx i j))
		    (set! nj (+ nj 1)))))
	    (set! ni (1+ ni)))))
    nmx))

(define (mixer-determinant mx)
  (let ((n (mus-length mx)))
    (if (= n 1) 
	(mat mx 0 0)
	(if (= n 2)
	    (- (* (mat mx 0 0) (mat mx 1 1))
	       (* (mat mx 0 1) (mat mx 1 0)))
	    (if (= n 3)
		(- (+ (* (mat mx 0 0) (mat mx 1 1) (mat mx 2 2))
		      (* (mat mx 0 1) (mat mx 1 2) (mat mx 2 0))
		      (* (mat mx 0 2) (mat mx 1 0) (mat mx 2 1)))
		   (+ (* (mat mx 0 0) (mat mx 1 2) (mat mx 2 1))
		      (* (mat mx 0 1) (mat mx 1 0) (mat mx 2 2))
		      (* (mat mx 0 2) (mat mx 1 1) (mat mx 2 0))))
		(let ((sum 0.0)
		      (sign 1))
		  (do ((i 0 (1+ i)))
		      ((= i n))
		    (let ((mult (mat mx 0 i)))
		      (if (not (= mult 0.0))
			  (set! sum (+ sum (* sign mult (mixer-determinant (sub-matrix mx 0 i))))))
		      (set! sign (- sign))))
		  sum))))))

(define (mixer-inverse mx)
  (let ((n (mus-length mx))
	(det (mixer-determinant mx)))
    (if (not (= det 0.0))
	(let ((nmx (make-zero-mixer n)))
	  (do ((i 0 (1+ i)))
	      ((= i n))
	    (do ((j 0 (1+ j)))
		((= j n))
	      (mixer-set! nmx j i (* (if (odd? (+ i j)) -1 1) 
				     (/ (mixer-determinant (sub-matrix mx i j)) 
					det)))))
	  nmx)
	#f)))

(define* (mixer-poly mx #:rest coeffs)
  (let* ((n (length coeffs))
	 (nmx (make-scalar-mixer (mus-length mx) (list-ref coeffs (1- n))))
	 (x (mixer-scale mx 1.0)))
    (do ((i (- n 2) (1- i)))
	((< i 0))
      (set! nmx (mixer+ nmx (mixer-scale x (list-ref coeffs i))))
      (set! x (mixer* mx x)))
    nmx))

(define (mixer-normal? mx)
  (mixer-equal? (mixer* mx (mixer-transpose mx)) ; conjugate transpose (hermitian) if complex
		(mixer* (mixer-transpose mx) mx)))

(define (mixer-orthogonal? mx)
  (mixer-equal? (mixer-transpose mx) 
		(mixer-inverse mx)))

(define (mixer-unitary? mx) (mixer-orthogonal? mx)) ; reals


;;; (define (vct-norm v1) (sqrt (dot-product v1 v1)))

(define (mixer-trace mx)
  (let ((sum 0.0)
	(n (mus-length mx)))
    (do ((i 0 (1+ i)))
	((= i n) sum)
      (set! sum (+ sum (mixer-ref mx i i))))))
