(use-modules (oop goops))
(use-modules (srfi srfi-1))

;; A matrix has only one slot, containing a list
;; of rows, where each row is a list of elements.

(define-class <matrix> ()
  (rows #:getter matrix-rows #:init-keyword #:rows))

(define (matrix rows)
  (make <matrix> #:rows rows))

(define (matrix-ref a i j)
  (list-ref (list-ref (matrix-rows a)
		      i)
	    j))

(define (display-matrix a port)
  (newline port)
  (for-each (lambda (row)
	      (display row port)
	      (newline port))
	    (matrix-rows a)))

(define (matrix-cols a)
  (apply map list (matrix-rows a)))

(define (matrix-map f a)
  (matrix (map (lambda (row) (map f row))
	       (matrix-rows a))))

(define (scale-matrix s a)
  (cond ((= s 1) a)
	(else (matrix-map (lambda (elem) (* s elem)) a))))

(define (zero-matrix n m)
  (matrix (make-list n (make-list m 0))))

(define (identity-matrix n)
  (matrix
   (list-tabulate n (lambda (i)
		      (append (make-list i 0)
			      (list 1)
			      (make-list (- n i 1) 0))))))

(define (row-vect . elems)  (matrix (list elems)))
(define (col-vect . elems)  (matrix (map list elems)))

(define (diagonal-matrix diag)
  (let ((n (length diag)))
    (matrix (map (lambda (i x)
		   (append (make-list i 0)
			   (list x)
			   (make-list (- n i 1) 0)))
		 (iota n) diag))))

(define (matrix-mult a b)
  (let ((a-rows (matrix-rows a))
	(b-cols (matrix-cols b)))
    (if (not (= (length (first a-rows))
		(length (first b-cols))))
	(throw 'matrix-mult a b))
    (matrix (map (lambda (a-row)
		   (map (lambda (b-col)
			  (apply + (map * a-row b-col)))
			b-cols))
		 a-rows))))

(define-method (display (a <matrix>) port)  (display-matrix a port))
(define-method (write   (a <matrix>) port)  (display-matrix a port))

(define-method (* (s <number>) (a <matrix>))  (scale-matrix s a))
(define-method (* (a <matrix>) (s <number>))  (scale-matrix s a))
(define-method (* (a <matrix>) (b <matrix>))  (matrix-mult a b))

;;
;; Fast fibonacci using matrix exponentiation
;;
;; [0 1] [a]   [ b ]
;; [1 1] [b] = [a+b]
;;
;; As shown above, the matrix above transforms a column vector (a b)
;; into (b a+b).  This column vector is the state of an iterative
;; algorithm to compute fibonacci, and the matrix is the
;; transformation done at each step.
;;
;; However, instead of multiplying the state vector times the matrix
;; at each step, we exponentiate the matrix to a power, where the
;; expononent N equals the number of steps we want to take.  This
;; yields a matrix that advances the state vector by N steps.  The
;; trick is to use the same fast exponentiation algorithm used on
;; integers.  The number of matrix multiplications performed is
;; proportional to the logarithm of the exponent.
;;
;; With the help of GOOPS, we use expt to exponentiate the matrix.  We
;; don't even have to define an expt method for matrices.  All that we
;; need is multiplication.  expt calls integer-expt, since the
;; exponent is an integer.  Thus this depends on integer-expt
;; accepting a base that is not a number.
;; 
(define (fib n)
  (let* ((m1 (matrix '((0 1)
		       (1 1))))
	 (m (expt m1 n))
	 (v (* m (col-vect 0 1))))
    (matrix-ref v 0 0)))