Jens Axel Soegaard wrote a comment on my last post asking for the Scheme code for comparison, so here it is. Note that the code provides a (moderately efficient, but conses too much) matrix inversion routine in pure Scheme (you want to look at the matrix module, for matrix-inverse)---some people may find that useful independently. Sorry it's so long, but there's a lot of background stuff to get through (it's part of a larger code base, too, but I think this is enough to be self-contained). Enjoy!

Here's the code for the generics that I use at the end of the matrix code:

#| base.ss: Generic operations and base language. Copyright (C) 2007 Will M. FarrThis program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |# (module base (lib "swindle.ss" "swindle") (provide (all-from-except (lib "swindle.ss" "swindle") + - * / add negate) (rename + mz:+) (rename - mz:-) (rename * mz:*) (rename / mz:/) (rename my-+ +) (rename my-- -) (rename my-* *) (rename my-/ /) (rename my-add add) sub mul div invert (rename my-negate negate)) (define my-+ (case-lambda (() 0) ((x) x) ((x y) (my-add x y)) ((x y . ys) (apply my-+ (my-add x y) ys)))) (define my-- (case-lambda ((x) (my-negate x)) ((x y) (sub x y)) ((x . ys) (sub x (apply my-+ ys))))) (define my-* (case-lambda (() 1) ((x) x) ((x y) (mul x y)) ((x y . ys) (apply my-* (mul x y) ys)))) (define my-/ (case-lambda ((x) (invert x)) ((x y) (div x y)) ((x . ys) (div x (apply my-* ys))))) (defgeneric (my-add x y)) (defmethod (my-add (x ) (y )) (+ x y)) (defgeneric (my-negate x)) (defmethod (my-negate (x )) (- x)) (defgeneric (sub x y)) (defmethod (sub (x ) (y )) (- x y)) (defgeneric (mul x y)) (defmethod (mul (x ) (y )) (* x y)) (defgeneric (invert x)) (defmethod (invert (x )) (/ x)) (defgeneric (div x y)) (defmethod (div (x ) (y )) (/ x y)))

Here's the matrix module (before we get to the Lorentz group):

#| matrix.ss: Basic matrix utilities. Copyright (C) Will M. FarrThis program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. |# (module matrix mzscheme (require (lib "42.ss" "srfi") (lib "contract.ss") (only "base.ss" struct-type->class add sub mul div invert negate defmethod )) (define mz:+ +) (define mz:- -) (define mz:* *) (define mz:/ /) (define-struct matrix (rows cols elts) #f) (define (struct-type->class struct:matrix)) (provide struct:matrix :matrix matrix-ec *small* ) (define (valid-matrix-row/c m) (and/c natural-number/c (between/c 0 (mz:- (matrix-rows m) 1)))) (define (valid-matrix-col/c m) (and/c natural-number/c (between/c 0 (mz:- (matrix-cols m) 1)))) (define (same-matrix-dimensions/c m) (flat-named-contract "same matrix dimensions" (let ((mc (matrix-cols m)) (mr (matrix-rows m))) (lambda (m2) (and (= (matrix-cols m2) mc) (= (matrix-rows m2) mr)))))) (define square-matrix/c (flat-named-contract "square matrix" (lambda (m) (= (matrix-rows m) (matrix-cols m))))) (provide/contract (rename my-make-matrix make-matrix (-> natural-number/c natural-number/c any/c matrix?)) (rename my-matrix matrix (->r ((r natural-number/c) (c natural-number/c)) elts (and/c list? (lambda (els) (= (length els) (mz:* r c)))) matrix?)) (matrix? (-> any/c boolean?)) (matrix-ref (->r ((m matrix?) (i (valid-matrix-row/c m)) (j (valid-matrix-col/c m))) any)) (matrix-set! (->r ((m matrix?) (i (valid-matrix-row/c m)) (j (valid-matrix-col/c m)) (x any/c)) any)) (matrix-mul (->r ((m1 matrix?) (m2 (and/c matrix? (flat-named-contract "compatible matrix dimensions" (lambda (m2) (= (matrix-cols m1) (matrix-rows m2))))))) matrix?)) (matrix-transpose (-> matrix? matrix?)) (matrix-vector-mul (->r ((m matrix?) (v (and/c vector? (flat-named-contract "compatible vector dimension" (lambda (v) (= (vector-length v) (matrix-cols m))))))) vector?)) (vector-matrix-mul (->r ((v vector?) (m (and/c matrix? (flat-named-contract "compatible vector dimension" (lambda (m) (= (vector-length v) (matrix-rows m))))))) vector?)) (matrix-add (->r ((m1 matrix?) (m2 (and/c matrix? (same-matrix-dimensions/c m1)))) matrix?)) (matrix-sub (->r ((m1 matrix?) (m2 (and/c matrix? (same-matrix-dimensions/c m1)))) matrix?)) (matrix-scale (->r ((s (or/c number? matrix?)) (m (if (number? s) matrix? number?))) matrix?)) (matrix-norm (-> (and/c integer? (>/c 0)) (-> matrix? number?))) (matrix-infinity-norm (-> matrix? number?)) (matrix-inverse (-> square-matrix/c square-matrix/c)) (matrix-identity (-> natural-number/c matrix?)) (matrix-copy (-> matrix? matrix?))) (define (my-make-matrix rows cols elt) (make-matrix rows cols (make-vector (mz:* rows cols) elt))) (define (matrix-ref m i j) (vector-ref (matrix-elts m) (mz:+ (mz:* i (matrix-cols m)) j))) (define (matrix-set! m i j x) (vector-set! (matrix-elts m) (mz:+ (mz:* i (matrix-cols m)) j) x)) (define (my-matrix rows cols . elts) (make-matrix rows cols (list->vector elts))) (define-syntax matrix-ec (syntax-rules () ((matrix-ec rrows ccols etc ...) (let ((rows rrows) (cols ccols)) (make-matrix rows cols (vector-of-length-ec (mz:* rows cols) etc ...)))))) (define-syntax :matrix (syntax-rules (index) ((:matrix cc var arg) (:vector cc var (matrix-elts arg))) ((:matrix cc var (index i j) arg) (:do cc (let ((m arg) (rows #f) (cols #f)) (set! rows (matrix-rows m)) (set! cols (matrix-cols m))) ((i 0) (j 0)) (< i rows) (let ((i+1 (mz:+ i 1)) (j+1 (mz:+ j 1)) (wrapping? #f)) (set! wrapping? (>= j+1 cols))) #t ((if wrapping? i+1 i) (if wrapping? 0 j+1)))))) (define (matrix-mul m1 m2) (let ((rows (matrix-rows m1)) (cols (matrix-cols m2)) (n (matrix-cols m1))) (matrix-ec rows cols (:range i rows) (:range j cols) (sum-ec (:range k n) (* (matrix-ref m1 i k) (matrix-ref m2 k j)))))) (define (matrix-transpose m) (matrix-ec (matrix-cols m) (matrix-rows m) (:range i (matrix-cols m)) (:range j (matrix-rows m)) (matrix-ref m j i))) (define (matrix-vector-mul m v) (vector-of-length-ec (matrix-rows m) (:range i (matrix-rows m)) (sum-ec (:range j (matrix-cols m)) (* (matrix-ref m i j) (vector-ref v j))))) (define (vector-matrix-mul v m) (vector-of-length-ec (matrix-cols m) (:range j (matrix-cols m)) (sum-ec (:range i (matrix-rows m)) (* (vector-ref v i) (matrix-ref m i j))))) (define (matrix-add m1 m2) (matrix-ec (matrix-rows m1) (matrix-cols m1) (:parallel (:matrix x1 m1) (:matrix x2 m2)) (+ x1 x2))) (define (matrix-sub m1 m2) (matrix-ec (matrix-rows m1) (matrix-cols m1) (:parallel (:matrix x1 m1) (:matrix x2 m2)) (- x1 x2))) (define (matrix-scale m s) (when (number? m) (let ((temp m)) (set! m s) (set! s temp))) (matrix-ec (matrix-rows m) (matrix-cols m) (:matrix x m) (* x s))) (define ((matrix-norm n) m) (expt (sum-ec (:matrix x m) (expt (abs x) n)) (/ n))) (define (matrix-infinity-norm m) (max-ec (:matrix x m) (abs x))) (define (scale-matrix-row! m i s) (do-ec (:range j (matrix-cols m)) (matrix-set! m i j (* (matrix-ref m i j) s)))) (define swap-matrix-rows! (case-lambda ((m i j pvec) (let ((oldi (vector-ref pvec i))) (vector-set! pvec i (vector-ref pvec j)) (vector-set! pvec j oldi)) (swap-matrix-rows! m i j)) ((m i j) (do-ec (:range k (matrix-cols m)) (let ((oldik (matrix-ref m i k))) (matrix-set! m i k (matrix-ref m j k)) (matrix-set! m j k oldik)))))) (define (add-scaled-matrix-row-to-row! m s i j) (do-ec (:range k (matrix-cols m)) (matrix-set! m j k (+ (matrix-ref m j k) (* s (matrix-ref m i k)))))) (define (pivot-row m j) (fold3-ec 'no-max (:range i j (matrix-rows m)) i (lambda (x) x) (lambda (i iold) (if (> (abs (matrix-ref m i j)) (abs (matrix-ref m iold j))) i iold)))) (define (matrix-pref m p i j) (matrix-ref m (vector-ref p i) j)) (define (matrix-identity n) (matrix-ec n n (:range i n) (:range j n) (if (= i j) 1 0))) (define (matrix-copy m) (matrix-ec (matrix-rows m) (matrix-cols m) (:matrix x m) x)) (define *small* (make-parameter 1.4901161193847656e-08 (lambda (sm) (if (not (>= sm 0)) (error '*small* "Cannot set to negative number: ~a" sm) sm)))) (define (matrix-inverse m-in) (let ((rows (matrix-rows m-in)) (cols (matrix-cols m-in))) (let ((m (matrix-copy m-in)) (id (matrix-identity rows)) (pvec (vector-of-length-ec rows (:range i rows) i))) (do-ec (:range i rows) (begin (let ((pi (pivot-row m i))) (swap-matrix-rows! m i pi pvec) (swap-matrix-rows! id i pi)) (let ((pelt (matrix-ref m i i))) (when (< (abs pelt) (*small*)) (raise-mismatch-error 'matrix-inverse "Singular or nearly singular (by *small* parameter) matrix " (cons pelt m-in))) (scale-matrix-row! m i (/ pelt)) (scale-matrix-row! id i (/ pelt)) (do-ec (:range j rows) (not (= j i)) (let ((s (- (matrix-ref m j i)))) (add-scaled-matrix-row-to-row! m s i j) (add-scaled-matrix-row-to-row! id s i j)))))) id))) (define (relative-size x norm) (abs (/ x norm))) (defmethod (add (m1 ) (m2 )) (matrix-add m1 m2)) (defmethod (sub (m1 ) (m2 )) (matrix-sub m1 m2)) (defmethod (mul (m ) (s )) (matrix-scale m s)) (defmethod (mul (s ) (m )) (matrix-scale m s)) (defmethod (mul (m ) (v )) (matrix-vector-mul m v)) (defmethod (mul (v ) (m )) (vector-matrix-mul v m)) (defmethod (mul (m1 ) (m2 )) (matrix-mul m1 m2)) (defmethod (negate (m )) (matrix-scale m -1)) (defmethod (invert (m )) (matrix-inverse m)) (defmethod (div (v ) (m )) (matrix-vector-mul (matrix-inverse m) v)) (defmethod (div (m ) (s )) (matrix-scale m (mz:/ s))))

And here's (finally), the code for SO(3,1)---see the documentation for SO31.plt for an explanation of how we parameterize the group if this seems a little opaque :).

#| lorentz-transform.ss: SO(3,1) matrix group. Copyright (C) 2007 Will M. Farr This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. |# (module lorentz-transform mzscheme (require "matrix.ss" (lib "42.ss" "srfi") (lib "contract.ss") (lib "math.ss") "simplex.ss" "simplex-map-utils.ss" (only "base.ss" struct-type->class defmethodmul div invert) (all-except (planet "table.ss" ("soegaard" "galore.plt" 3)) map ormap andmap for-each values equal? union)) (define-struct lorentz-transform-field (table) #f) (define-struct lorentz-transform-element (params inverse-params matrix inverse-matrix) #f) (define (struct-type->class struct:lorentz-transform-field)) (define (struct-type->class struct:lorentz-transform-element)) (provide eta lt-params/c struct:lorentz-transform-field struct:lorentz-transform-element ) (define lt-params/c (and/c vector? (flat-named-contract "6-element vector" (lambda (v) (= (vector-length v) 6))))) (define zero-simplex/c (and/c simplex? (flat-named-contract "simplex-dimension = 0" (lambda (s) (= (simplex-dimension s) 0))))) (define ltf-alist/c (listof (cons/c zero-simplex/c (or/c matrix? lt-params/c lorentz-transform-element?)))) (provide/contract (matrix-eta-trace (-> matrix? number?)) (lorentz-matrix (-> number? number? number? number? number? number? matrix?)) (make-lorentz-matrix (-> lt-params/c matrix?)) (lorentz-matrix->params (-> matrix? lt-params/c)) (Dlorentz-matrix (-> lt-params/c (vectorof matrix?))) (lorentz-transform-field? (-> any/c boolean?)) (lorentz-transform-element? (-> any/c boolean?)) (params->lorentz-transform-element (-> lt-params/c lorentz-transform-element?)) (matrix->lorentz-transform-element (-> matrix? lorentz-transform-element?)) (lorentz-transform-element-params (-> lorentz-transform-element? lt-params/c)) (lorentz-transform-element-inverse-params (-> lorentz-transform-element? lt-params/c)) (lorentz-transform-element-matrix (-> lorentz-transform-element? matrix?)) (lorentz-transform-element-inverse-matrix (-> lorentz-transform-element? matrix?)) (lorentz-transform-element-compose (-> lorentz-transform-element? lorentz-transform-element? lorentz-transform-element?)) (lorentz-transform-element-invert (-> lorentz-transform-element? lorentz-transform-element?)) (lorentz-transform-element-transform-vector (-> lorentz-transform-element? vector? vector?)) (rename my-make-lorentz-transform-field make-lorentz-transform-field (-> ltf-alist/c lorentz-transform-field?)) (rename my-lorentz-transform-field lorentz-transform-field (->* () ltf-alist/c (lorentz-transform-field?))) (lorentz-transform-field-element (-> lorentz-transform-field? zero-simplex/c lorentz-transform-element?)) (lorentz-transform-field-compose (-> lorentz-transform-field? lorentz-transform-field? lorentz-transform-field?)) (lorentz-transform-field-fold (-> (-> zero-simplex/c lorentz-transform-element? any/c any) any/c lorentz-transform-field? any))) (define (matrix-eta-trace M) (+ (- (matrix-ref M 0 0)) (matrix-ref M 1 1) (matrix-ref M 2 2) (matrix-ref M 3 3))) (define Rx-gen (matrix 4 4 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 1 0)) (define Ry-gen (matrix 4 4 0 0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0)) (define Rz-gen (matrix 4 4 0 0 0 0 0 0 -1 0 0 1 0 0 0 0 0 0)) (define Bx-gen (matrix 4 4 0 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0)) (define By-gen (matrix 4 4 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0 0 0)) (define Bz-gen (matrix 4 4 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0)) (define eta (matrix 4 4 -1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1)) (define (Rx x) (matrix 4 4 1 0 0 0 0 1 0 0 0 0 (cos x) (- (sin x)) 0 0 (sin x) (cos x))) (define (Ry y) (matrix 4 4 1 0 0 0 0 (cos y) 0 (- (sin y)) 0 0 1 0 0 (sin y) 0 (cos y))) (define (Rz z) (matrix 4 4 1 0 0 0 0 (cos z) (- (sin z)) 0 0 (sin z) (cos z) 0 0 0 0 1)) (define (Bx x) (matrix 4 4 (cosh x) (- (sinh x)) 0 0 (- (sinh x)) (cosh x) 0 0 0 0 1 0 0 0 0 1)) (define (By y) (matrix 4 4 (cosh y) 0 (- (sinh y)) 0 0 1 0 0 (- (sinh y)) 0 (cosh y) 0 0 0 0 1)) (define (Bz z) (matrix 4 4 (cosh z) 0 0 (- (sinh z)) 0 1 0 0 0 0 1 0 (- (sinh z)) 0 0 (cosh z))) (define gens (list Bx-gen By-gen Bz-gen Rx-gen Ry-gen Rz-gen)) (define mat-funs (list Bx By Bz Rx Ry Rz)) (define rgens (reverse gens)) (define rmat-funs (reverse mat-funs)) (define (make-lorentz-matrix param-v) (let ((! vector-ref) (* matrix-mul) (p param-v)) (* (Bx (! p 0)) (* (By (! p 1)) (* (Bz (! p 2)) (* (Rx (! p 3)) (* (Ry (! p 4)) (Rz (! p 5))))))))) (define (lorentz-matrix vx vy vz x y z) (make-lorentz-matrix (vector vx vy vz x y z))) (define (Dlorentz-matrix params) (let* ((lparams (vector->list params)) (mats (map (lambda (f x) (f x)) mat-funs lparams)) (rmats (reverse mats))) (vector-of-length-ec 6 (:range i 6) (fold-ec (matrix-identity 4) (:parallel (:list M rmats) (:list g rgens) (:range j 5 -1 -1)) (if (= i j) (matrix-mul M g) M) matrix-mul)))) (define (asinh z) (log (+ z (sqrt (+ (sqr z) 1))))) (define (acosh z) (log (+ z (* (sqrt (- z 1)) (sqrt (+ z 1)))))) (define (lorentz-matrix->params M) (let* ((sinh-vz (- (matrix-ref M 3 0))) (vz (asinh sinh-vz)) (cosh-vz (cosh vz)) (sinh-vy (- (/ (matrix-ref M 2 0) cosh-vz))) (vy (asinh sinh-vy)) (cosh-vy (cosh vy)) (sinh-vx (- (/ (matrix-ref M 1 0) (* cosh-vy cosh-vz)))) (vx (asinh sinh-vx))) (let ((new-M (matrix-mul (Bz (- vz)) (matrix-mul (By (- vy)) (matrix-mul (Bx (- vx)) M))))) (let* ((sin-y (- (matrix-ref new-M 1 3))) (y (asin sin-y)) (cos-y (cos y)) (sin-z (- (/ (matrix-ref new-M 1 2) cos-y))) (z (asin sin-z)) (cos-z (cos z)) (cos-x (/ (matrix-ref new-M 3 3) cos-y)) (sin-x (/ (+ (matrix-ref new-M 3 2) (* sin-y sin-z cos-x)) cos-z)) (x (asin sin-x))) (vector vx vy vz x y z))))) (define (my-lorentz-transform-field . maps) (my-make-lorentz-transform-field maps)) (define (params->lorentz-transform-element p) (let* ((M (make-lorentz-matrix p)) (M-inv (matrix-inverse M)) (p-inv (lorentz-matrix->params M-inv))) (make-lorentz-transform-element p p-inv M M-inv))) (define (matrix->lorentz-transform-element M) (let* ((M-inv (matrix-inverse M)) (p (lorentz-matrix->params M)) (p-inv (lorentz-matrix->params M-inv))) (make-lorentz-transform-element p p-inv M M-inv))) (define (my-make-lorentz-transform-field ltf-alist) (let ((elt-alist (map (lambda (s-p-or-M) (let ((p-or-M (cdr s-p-or-M))) (cons (car s-p-or-M) (cond ((vector? p-or-M) (params->lorentz-transform-element p-or-M)) ((matrix? p-or-M) (matrix->lorentz-transform-element p-or-M)) (else p-or-M))))) ltf-alist))) (make-lorentz-transform-field (alist->ordered simplex-compare elt-alist)))) (define (lorentz-transform-field-element ltf s) (lookup s (lorentz-transform-field-table ltf) (lambda () (make-lorentz-transform-element (make-vector 6 0) (make-vector 6 0) (matrix-identity 4) (matrix-identity 4))) values)) (define (lorentz-transform-element-compose elt1 elt2) (let* ((M1 (lorentz-transform-element-matrix elt1)) (M2 (lorentz-transform-element-matrix elt2)) (M (matrix-mul M1 M2)) (p (lorentz-matrix->params M)) (M-inv (matrix-inverse M)) (p-inv (lorentz-matrix->params M-inv))) (make-lorentz-transform-element p p-inv M M-inv))) (defmethod (mul (e1 ) (e2 )) (lorentz-transform-element-compose e1 e2)) (define (lorentz-transform-element-transform-vector elt v) (matrix-vector-mul (lorentz-transform-element-matrix elt) v)) (defmethod (mul (e ) (v )) (lorentz-transform-element-transform-vector e v)) (define (lorentz-transform-element-invert elt) (make-lorentz-transform-element (lorentz-transform-element-inverse-params elt) (lorentz-transform-element-params elt) (lorentz-transform-element-inverse-matrix elt) (lorentz-transform-element-matrix elt))) (defmethod (invert (e )) (lorentz-transform-element-invert e)) (defmethod (div (v ) (e )) (matrix-vector-mul (lorentz-transform-element-inverse-matrix e) v)) (define lorentz-transform-field-compose (let ((table-compose (simplex-map-lift-combiners (lambda () (make-ordered simplex-compare)) insert (lambda (sl el sr er accu) (insert sl (lorentz-transform-element-compose el er) accu)) insert))) (lambda (ltf1 ltf2) (make-lorentz-transform-field (table-compose (lorentz-transform-field-table ltf1) (lorentz-transform-field-table ltf2)))))) (defmethod (mul (ltf1 ) (ltf2 )) (lorentz-transform-field-compose ltf1 ltf2)) (define (lorentz-transform-field-fold f start ltf) (fold f start (lorentz-transform-field-table ltf))))