行列演算

ちょいと3Dを勉強しなおしてるわけですが計算用の道具でよさげなのが手元にないのです。
maximaがいい感じなのですが検算用に使いづらい。一応ループや条件分岐などのロジックは書けるようなのですが、中途半端にlispっぽさを引きずってて馴染めそうに無いのです。

一応C++で書いた行列演算などのコードはあるのですがちょっと値をかえるだけでリコンパイルとかやってられんので（入力用のI／Fを作るのはもっとやってられなかった）手軽に実行できる言語とかでないかなぁ、と探した結果、lispに白羽の矢がたちました。

で、行列関連のコードを探したわけですが処理系依存してたりインストールが面倒そうだったり・・・

実際のところ、和、差、積、転置くらいが出来ればよく、逆行列とかもいらなかったので勉強かねて自分で実装してみました。
ほんとに出来るのは検算用の計算くらいです。

(defpackage :math.geometry
  (:use :cl)
  (:export :transpose
	   :add
	   :add-scalar
	   :sub
	   :sub-scalar
	   :multiply
	   :multiply-scalar
	   :column
	   :row
	   :element
	   :make-matrix
	   :make-scalar-matrix))

(in-package :math.geometry)

(defun make-matrix (row column)
  (loop repeat row collect (loop repeat column collect 0)))

#+test (tree-equal (make-matrix 3 4)
		   '((0 0 0 0)
		     (0 0 0 0)
		     (0 0 0 0)))

(defun make-scalar-matrix (size value)
  (loop for i below size
     collect (loop for j below size collect (if (eq i j) value 0))))

#+test (tree-equal (make-scalar-matrix 4 2)
		   '((2 0 0 0)
		     (0 2 0 0)
		     (0 0 2 0)
		     (0 0 0 2)))

(defun row (mat)
  (length mat))

#+test (and (equal (row '((0 0) (0 0))) 2)
	    (equal (row '((0 0 0) (0 0 0))) 2))

(defun column (mat)
  (let ((line (car mat)))
    (if (listp line)
	(length line)
	1)))

#+test (and (equal (column '((0 0) (0 0))) 2)
	    (equal (column '((0 0 0 0) (0 0 0 0))) 4))

(defun element (mat row column)
  (nth column (nth row mat)))

#+test (let ((mat '((1 2) (3 4))))
	 (and (equal (element mat 1 1) 4)
	      (equal (element mat 1 0) 3)
	      (equal (element mat 0 1) 2)
	      (equal (element mat 0 0) 1)))

(defun set-element (mat row column value)
  (setf (nth column (nth row mat)) value))

(defsetf element set-element)

#+test (let ((mat '((1 2 3) (4 5 6) (7 8 9))))
	 (progn
	   (setf (element mat 0 0) 11)
	   (setf (element mat 2 1) 16)
	   (setf (element mat 2 2) 19)
	   (tree-equal mat '((11 2 3)
			     (4 5 6)
			     (7 16 19)))))

(defun transpose (mat)
  (let* ((column (column mat))
	 (row (row mat))
	 (result (make-matrix column row)))
    (loop for i below row
       do (loop for j below column
	       do (setf (element result j i)
		      (element mat i j))))
    result))

#+test (and
	(tree-equal (transpose '((1 2) (3 4)))
		    '((1 3) (2 4)))
	(tree-equal (transpose '((1 2)))
		    '((1) (2))))

(defun add (lhs rhs)
  (unless 
      (and (eq (column lhs) (column rhs))
	   (eq (row lhs) (row rhs)))
    (error "not equals lhs and rhs column or/and row"))
  (let* ((column (column lhs))
	 (row (row lhs))
	 (result (make-matrix row column)))
    (loop for i below column
       do (loop for j below row
	     do (setf (element result j i)
		      (+ (element lhs j i)
			 (element rhs j i)))))
    result))

#+test (tree-equal (add '((1 2 3) (4 5 6)) '((7 8 9) (10 11 12)))
		   '((8 10 12) (14 16 18)))

(defun square-matrixp (mat)
  (eq (column mat) (row mat)))

#+test (and (square-matrixp (make-scalar-matrix 2 1))
	    (not (square-matrixp '((1 2 3) (4 5 6)))))

(defun add-scalar (lhs rhs)
  (unless (square-matrixp lhs)
    (error "not equals lhs's column and row"))
  (let ((rhs (make-scalar-matrix (column lhs) rhs)))
    (add lhs rhs)))

#+test (tree-equal 
	(add-scalar '((1 2) (3 4)) 5)
	'((6 2) (3 9)))

(defun sub (lhs rhs)
  (unless 
      (and (eq (column lhs) (column rhs))
	   (eq (row lhs) (row rhs)))
    (error "not equals lhs and rhs column or/and row"))
  (let* ((column (column lhs))
	 (row (row lhs))
	 (result (make-matrix row column)))
    (loop for i below column
       do (loop for j below row
	     do (setf (element result j i)
		      (- (element lhs j i)
			 (element rhs j i)))))
    result))

#+test (tree-equal (sub '((1 2 3) (4 5 6)) '((7 8 9) (10 11 12)))
		   '((-6 -6 -6) (-6 -6 -6)))

(defun sub-scalar (lhs rhs)
  (unless (square-matrixp lhs)
    (error "not equals lhs's column and row"))
  (let ((rhs (make-scalar-matrix (column lhs) rhs)))
    (sub lhs rhs)))

#+test (tree-equal 
	(sub-scalar '((1 2) (3 4)) 5)
	'((-4 2) (3 -1)))

(defun multiply (lhs rhs)
  (unless (eq (column lhs) (row rhs))
    (error "not equal lhs's column and rhs's row"))
  (let* ((column (column rhs))
	 (row (row lhs))
	 (max-k (column lhs))
	 (result (make-matrix row column)))
    (loop for j below column do
	 (loop for i below row do
	      (setf (element result i j)
		    (loop for k below max-k sum
			 (* (element lhs i k)
			    (element rhs k j))))))
    result))

#+test (and
	(tree-equal (multiply '((1 2) (3 4)) '((5 6) (7 8)))
		    '((19 22) (43 50)))
	(tree-equal (multiply '((5 6) (7 8)) '((1 2) (3 4)))
		    '((23 34) (31 46)))
	(tree-equal (multiply '((1 2) (3 4)) '((5) (6)))
		    '((17) (39))))

(defun multiply-scalar (lhs rhs)
  (unless (square-matrixp lhs)
    (error "not square matrix for lhs."))
  (multiply lhs (make-scalar-matrix (column lhs) rhs)))

#+test (tree-equal (multiply-scalar '((1 2) (3 4)) 5)
		   '((5 10) (15 20)))

・・・短けぇ!!
多少のエラーチェックを含んだ行列演算のコードがたった100行ちょいかぁ・・・
あ、でもC++で組んだ同レベルのコードも200行ちょいくらいだな。
それでも実装にかかった時間は断然lispの方が短いや。

で、こんな感じで検算用に使ってます。

(use-package :math.geometry)

(defun projection (left right top bottom near far)
  (let ((width_distance (- right left))
	(width_center (+ right left))
	(height_distance (- top bottom))
	(height_center (+ top bottom))
	(depth_distance (- far near))
	(depth_center (+ far near)))
    (let ((width_ratio (/ width_center width_distance))
	  (height_ratio (/ height_center height_distance))
	  (depth_ratio (/ depth_center depth_distance)))
      `((,(/ (* 2 near) width_distance) 0 0 0)
	(0 ,(/ (* 2 near) height_distance) 0 0)
	(,width_ratio ,height_ratio ,depth_ratio -1)
	(0 0 ,(/ (* -2 far near) depth_distance) 0)))))

(let ((proj (projection -1 1 1 -1 -1 -100)))
  (loop for v from -1 downto -100 collect
       (multiply `((1 0 ,v 1)) proj)))

とまあ、3D用の透視投影の検算したかっただけなんですけどね。

※ 行と列の扱いがおかしいところがいくつか散見されたので直したつもり。あくまでつもり(笑)

※ うがー、やっぱりダメな場所があった・・・elementの行と列が逆になってる・・・
俺みたいなおバカには関連メソッドは一ヶ所に強制的に書かせるクラス系の表記のほうがいいのかもしらん・・・

※ テストコードつきで上げ直した。これで一応動くはず・・・