導入
逆行列を作るプログラムを書きます。数だけでなく変数を含む行列にも対応します。
プログラム
Lisp ぽく nil を定義（個人的嗜好）。
(define nil '())
次はスタックを扱う手続き。スタックは記憶方法の一つで、最後に格納した値を最初に取り出せる仕組みのことです(LIFO: Last In First Out)。コンピュータにおいては関数を呼び出した後、元の呼び出し元に戻る場所を一時的に覚えるのに使われたりします。
(define *stack* nil)
(define level-ratio 100)
(define (last-pushed)
(if (empty?)
nil
(car *stack*)))
(define (after-popped)
(if (empty?)
nil
(cdr *stack*)))
(define (init)
(set! *stack* nil))
(define (empty?)
(null? *stack*))
(define (push x)
(set! *stack* (cons x *stack*)))
(define (pop)
(if (empty?)
nil
(let ((result (last-pushed)))
(set! *stack* (after-popped))
result)))
(define (pop-all level)
(define (iter operators)
(if (or (empty?) (< (priority (last-pushed)) (* level level-ratio)))
operators
(iter (append operators (cons (operator (pop)) nil)))))
(iter nil))
次は逆ポーランド記法を処理するための手続き。逆ポーランド記法は、計算式を表現するための方式の一つで後置記法とも言います。コンピュータで計算式を扱う場合にいろいろ都合が良いのでよく使われます。
(define (make-priority op level)
(let ((priority
(cond ((eq? op '**) 30)
((eq? op '*) 20)
((eq? op '/) 20)
((eq? op '+) 10)
((eq? op '-) 10)
(else 0))))
(cons op (+ priority (* level level-ratio)))))
(define (operator p)
(car p))
(define (priority p)
(if (null? p)
0
(cdr p)))
(define (operator? x)
(or (eq? x '+) (eq? x '-) (eq? x '*) (eq? x '/) (eq? x '**)))
(define (exchange-operator op level)
(define (iter ops)
(if (or (empty?) (< (priority (last-pushed)) (priority (make-priority op level))))
ops
(iter (append ops (cons (operator (pop)) nil)))))
(iter nil))
(define (exp->rpn exp)
(define (iter e r level)
(if (null? e)
(append r (pop-all level))
(if (list? (car e))
(iter (cdr e) (iter (car e ) r (+ level 1)) level)
(cond ((operator? (car e))
(let ((ops (exchange-operator (car e) level)))
(push (make-priority (car e) level))
(iter (cdr e) (append r ops) level)))
(else (iter (cdr e) (append r (list (car e))) level))))))
(init)
(iter exp nil 0))
(define (rpn->s-exp rpn)
(define (iter r)
(if (null? r)
(pop)
(let ((e (car r)))
(cond ((eq? e '+) (let ((addend (pop)))
(push (make-sum (pop) addend))))
((eq? e '-) (let ((minuend (pop)))
(push (make-sub (pop) minuend))))
((eq? e '*) (let ((multiplier (pop)))
(push (make-product (pop) multiplier))))
((eq? e '/) (let ((divisor (pop)))
(if (and (number? divisor) (= divisor 0))
(error "attempt to calculate a division by zero")
(push (make-division (pop) divisor)))))
((eq? e '**) (let ((exponent (pop)))
(push (make-exponentiation (pop) exponent))))
(else (push e)))
(iter (cdr r)))))
(init)
(iter rpn))
微分手続きとそれに関連する手続き。微分とは与えられた関数の導関数を求めることです。導関数は関数の変化率を表しています。ちなみに、世界初のコンピュータである ENIAC は弾道計算（微分を使います）を行うために製作されました。
(define (deriv exp var)
(define (deriv-iter exp var)
(cond ((number? exp) 0)
((variable? exp)
(if (same-variable? exp var) 1 0))
((sum? exp)
(make-sum (deriv-iter (addend exp) var)
(deriv-iter (augend exp) var)))
((sub? exp)
(make-sub (deriv-iter (minuend exp) var)
(deriv-iter (subtrahend exp) var)))
((product? exp)
(make-sum
(make-product (multiplicand exp)
(deriv-iter (multiplier exp) var))
(make-product (deriv-iter (multiplicand exp) var)
(multiplier exp))))
((exponentiation? exp)
(make-product
(make-product (exponent exp)
(make-exponentiation (base exp) (make-sum (exponent exp) -1)))
(deriv-iter (base exp) var)))
(else
(error "unknown expression type -- DERIV" exp))))
(deriv-iter (rpn->s-exp (exp->rpn exp)) var))
(define (variable? x) (symbol? x))
(define (same-variable? v1 v2)
(and (variable? v1) (variable? v2) (eq? v1 v2)))
(define (=number? exp num)
(and (number? exp) (= exp num)))
(define (sort-sum-values lis)
(let ((num (apply + (filter number? lis)))
(vals (filter symbol? lis)))
(if (= num 0)
(sort vals)
(append `(,num) (sort vals)))))
(define (sort-product-values lis)
(let ((num (apply * (filter number? lis)))
(vals (filter symbol? lis)))
(if (= num 1)
(sort vals)
(append `(,num) (sort vals)))))
(define (opti-sum a1 a2)
(cond ((and (sum? a1) (sum? a2)) (append '(+) (sort-sum-values (append (cdr a1) (cdr a2)))))
((sum? a1) (append '(+) (sort-sum-values (append (cdr a1) `(,a2)))))
((sum? a2) (append '(+) (sort-sum-values (append `(,a1) (cdr a2)))))
(else (list '+ a1 a2))))
(define (opti-product a1 a2)
(cond ((and (product? a1) (product? a2)) (append '(*) (sort-product-values (append (cdr a1) (cdr a2)))))
((product? a1) (append '(*) (sort-product-values (append (cdr a1) `(,a2)))))
((product? a2) (append '(*) (sort-product-values (append `(,a1) (cdr a2)))))
(else (list '* a1 a2))))
(define (make-sum a1 a2)
(cond ((=number? a1 0) a2)
((=number? a2 0) a1)
((and (number? a1) (number? a2)) (+ a1 a2))
(else (opti-sum a1 a2))))
(define (make-product m1 m2)
(cond ((or (=number? m1 0) (=number? m2 0)) 0)
((=number? m1 1) m2)
((=number? m2 1) m1)
((and (number? m1) (number? m2)) (* m1 m2))
(else (opti-product m1 m2))))
(define (make-sub a1 a2)
(cond ((=number? a1 0) (* -1 a2))
((=number? a2 0) a1)
((and (number? a1) (number? a2)) (- a1 a2))
(else (list '- a1 a2))))
(define (make-division m1 m2)
(cond ((=number? m1 0) 0)
((=number? m2 0) (error "attempt to calculate a division by zero"))
((=number? m2 1) m1)
((and (number? m1) (number? m2)) (/ m1 m2))
(else (list '/ m1 m2))))
(define (make-exponentiation e1 e2)
(cond ((=number? e2 0) 1)
((=number? e2 1) e1)
(else (list '** e1 e2))))
(define (sum? x)
(and (pair? x) (eq? (car x) '+)))
;; Japanese style
;; (+ augend addend)
(define (augend s) (cadr s))
(define (addend s) (caddr s))
(define (sub? x)
(and (pair? x) (eq? (car x) '-)))
;; Japanese style
;; (- minuend subtrahend)
(define (minuend s) (cadr s))
(define (subtrahend s) (caddr s))
(define (product? x)
(and (pair? x) (eq? (car x) '*)))
;; Japanese style
;; (* multiplicand multiplier)
(define (multiplicand p) (cadr p))
(define (multiplier p) (caddr p))
(define (exponentiation? x)
(and (pair? x) (eq? (car x) '**)))
;; Japanese style
;; (** base exponent)
(define (base e) (cadr e))
(define (exponent e) (caddr e))
最後は逆行列を扱う手続きです。ピボット選択に対応しています。逆行列の計算方法を紹介しているサイトや動画をみて勉強しながらコードを書きました。
(define (make-key r c)
`(,r ,c))
(define (square-matrix? matrix)
(let ((row (length matrix)))
(equal? (make-list row row) (map length matrix))))
(define (make-unit-matrix matrix)
(define (iter i lst mat)
(if (= i 0)
mat
(iter (- i 1) (append (cdr lst) (cons 0 nil)) (cons lst mat))))
(let ((len (length matrix)))
(iter len (append (make-list (- len 1) 0) '(1)) nil)))
(define (matrix->hash-table matrix)
(define ht (make-hash-table 'equal?))
(define (col-iter r c cols)
(if (null? cols)
'done
(begin
(hash-table-set! ht (make-key r c) (car cols))
(col-iter r (+ c 1) (cdr cols)))))
(define (row-iter r rows)
(if (null? rows)
ht
(begin
(col-iter r 1 (car rows))
(row-iter (+ r 1) (cdr rows)))))
(row-iter 1 matrix))
(define (hash-table->row ht size row)
(define (iter c row-list)
(if (> c size)
row-list
(iter (+ c 1) (append row-list (cons (hash-table-get ht (make-key row c)) nil)))))
(iter 1 nil))
(define (hash-table->matrix ht size)
(define (r-iter r matrix)
(if (> r size)
matrix
(r-iter (+ r 1) (append matrix (cons (hash-table->row ht size r) nil)))))
(r-iter 1 nil))
(define (solve-main ht1 ht2 size)
(define (change-row row1 row2)
(define (iter c r1 r2 r3 r4)
(if (> c size)
'done
(begin
(hash-table-set! ht1 (make-key row1 c) (car r2))
(hash-table-set! ht1 (make-key row2 c) (car r1))
(hash-table-set! ht2 (make-key row1 c) (car r4))
(hash-table-set! ht2 (make-key row2 c) (car r3))
(iter (+ c 1) (cdr r1) (cdr r2) (cdr r3) (cdr r4)))))
(let ((r1 (hash-table->row ht1 size row1))
(r2 (hash-table->row ht1 size row2))
(r3 (hash-table->row ht2 size row1))
(r4 (hash-table->row ht2 size row2)))
(iter 1 r1 r2 r3 r4)))
(define (irekae pivot)
(define (r-iter r max-r max-val)
(if (> r size)
(change-row pivot max-r)
(let ((val (hash-table-get ht1 (make-key r pivot))))
(if (and (number? val) (> (abs val) max-val))
(r-iter (+ r 1) r val)
(r-iter (+ r 1) max-r max-val)))))
(r-iter pivot pivot (hash-table-get ht1 (make-key pivot pivot))))
(define (taikaku pivot)
(define (iter val c)
(if (> c size)
'done
(begin
(hash-table-set! ht1 (make-key pivot c) (make-division (hash-table-get ht1 (make-key pivot c)) val))
(hash-table-set! ht2 (make-key pivot c) (make-division (hash-table-get ht2 (make-key pivot c)) val))
(iter val (+ c 1)))))
(let ((val (hash-table-get ht1 (make-key pivot pivot))))
(and (not (= val 1)) (iter val 1))))
(define (hitaikaku pivot)
(define (c-iter val row col)
(if (> col size)
'done
(begin
(hash-table-set! ht1
(make-key row col)
(make-sum (hash-table-get ht1 (make-key row col))
(make-product (hash-table-get ht1 (make-key pivot col))
val)))
(hash-table-set! ht2
(make-key row col)
(make-sum (hash-table-get ht2 (make-key row col))
(make-product (hash-table-get ht2 (make-key pivot col))
val)))
(c-iter val row (+ col 1)))))
(define (r-iter row)
(if (> row size)
'done
(begin
(let ((val (hash-table-get ht1 (make-key row pivot))))
(and (not (= row pivot)) (c-iter (make-product -1 val) row 1)))
(r-iter (+ row 1)))))
(r-iter 1))
(define (iter pivot)
(if (> pivot size)
ht2
(begin
(irekae pivot)
(taikaku pivot)
(hitaikaku pivot)
(iter (+ pivot 1)))))
(iter 1)
)
(define (solve matrix)
(let ((size (length matrix)))
(hash-table->matrix (solve-main (matrix->hash-table matrix)
(matrix->hash-table (make-unit-matrix matrix))
size)
size)))
(define (inverse-matrix matrix)
(if (not (square-matrix? matrix))
(error "Require square matrix.")
(solve matrix)))
(define ** expt)
実行
逆行列を求めてみます。
gosh> (inverse-matrix '((1 3 2) (-1 0 1) (2 3 0)))
((1 -2 -1) (-2/3 4/3 1) (1 -1 -1))
逆行列の逆行列は元の行列になります。
gosh> (inverse-matrix (inverse-matrix '((1 3 2) (-1 0 1) (2 3 0))))
((1 3 2) (-1 0 1) (2 3 0))
変数混じりの行列ついても扱うことができます。
マイナスの単項演算子がつく変数の表現は、ちょっといまいちです。
gosh> (inverse-matrix '((1 0 0) (a 1 0) (b c 1)))
((1 0 0) ((* -1 a) 1 0) ((+ (* -1 b) (* a c)) (* -1 c) 1))