Pure Scheme Code for SO(3,1)

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. 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 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. 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 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 defmethod  mul 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))))

Why Use the FFI for Simple Numerical Operations in MzScheme?

Brad Lucier (of Numerical Partial Differential Equations in Scheme fame) posted a comment on my last post announcing the code for numerical manipulation of Lorentz group elements which reads:

I'm surprised that you decided to use the C FFI to manipulate 4x4 matrices; they are so small that it would seem better just to use mzscheme's numerical facilities which are getting better all the time.

This is an interesting question to answer, so I'm going to devote this post to answering it instead of leaving a short response in the comments section.

Before writing this code (which uses some C to do the manipulations), I wrote some code which used purely Scheme to perform the following sequence of operations, given a set of parameters for an SO(3,1) element:

1. Compute the matrix associated with the parameters, M.
2. Invert that matrix to produce another matrix M-inv.
3. Compute the parameters associated with M-inv, p-inv.

This is a natural sequence of operations to perform on the group---it's just the group inverse function on the parameters. The timings for 10K repeats of this operation are:

Code Style	Total CPU Time	GC Time
Pure Scheme	51 sec	27 sec
Mixed Scheme/C	4 sec	2.8 sec

I think the large GC time for the Scheme/C code is probably due to the fact that each call goes through several contract checks, which probably dominate the run-time and allocation. Clearly, there's still an advantage to doing significant manipulations in C! I wish that it was practical to stay completely within PLT Scheme for doing this sort of work, but I'm afraid that it isn't. I suspect that Brad's implementation of choice (I assume), Gambit, would do much better on this sort of comparison.

To philosophize: the role I see for Scheme in my computational world is as an excellent driver to glue a bunch of simple C or Fortran computational elements together. (This isn't my idea---have a look at this for a similar approach with Python.) PLT Scheme is my implementation of choice for this because of its excellent FFI, it's extensive libraries (particularly its support of software contracts), and it's relatively easy C functions for memory allocation and Scheme object manipulation. The contracts, in particular, are a big help, because they let me avoid doing any argument checking on the C side (where I want to minimize coding). Memory allocation is nice, because with a combination of mzc --xform and scheme_malloc[_atomic], I never again have to mess with malloc and free in the C that I write. I've written this code with an eye toward including it in my thesis work, which will also be written in this style, so it seemed natural to do the manipulations in C.

This is, of course, my view. People like Brad, Marc, Manuel, and Jeff have a different vision, and have done tremendous work in optimizing Scheme performance on numeric codes. I've tried that approach, and always find myself back with PLT's scheme and C.

SO(3,1) (Lorentz Group) Manipulations from PLT Scheme

I just released a PlaneT package which manipulates elements of the Lorentz group (here). Get it using

(require (planet "SO31.ss" ("wmfarr" "SO31.plt" 1 0)))

At the moment, it only runs on Mac OS 10.3 or 10.4, ppc, because it uses some C code I wrote (for efficiency), and I haven't had the patience to understand all the hoops I need to jump through to get PLT Scheme to compile this to a library on every system.

If someone on Windows or Linux (or an x86 mac) really wants to use the package, it's pretty easy to compile manually (after downloading and installing from PlaneT): just compile the SO31.c file in the SO31c subdirectory of the package to a dynamic library in the same directory named "libSO31.", linking with BLAS and LAPACK.

I hope people find all this linear algebra stuff useful (it certainly will be for my thesis code).

Linear Algebra in PLT Scheme

Last night, I posted a library which uses BLAS and LAPACK from within PLT Scheme to perform simple numerical linear algebra.

In doing so, I've learned a lot about the FFI for PLT Scheme. The conclusion: it's great! In particular, the "custom _fun types", discussed here are tremendously helpful! Particularly when interfacing to Fortran libraries, where every argument to a function must be a pointer. For example, here's how I make a call to dgesv from LAPACK to solve linear systems:

  (define matrix-solve-many
    (get-ffi-obj 'dgesv_ *lapack*
                 (_fun (m b) ::
                       ((_ptr i _int) = (matrix-rows m))
                       ((_ptr i _int) = (matrix-cols b))
                       (_matrix = (matrix-copy m))
                       ((_ptr i _int) = (matrix-rows m))
                       (_u32vector o (matrix-rows m))
                       (x : _matrix = (matrix-copy b))
                       ((_ptr i _int) = (matrix-rows b))
                       (_ptr o _int) ->
                       _void ->
                       x)))

all the (_ptr o ...) calls allocate the appropriate memory blocks and pass pointers to them for dgesv, and then handle the translation back to the scheme world automatically. I didn't have to write any C code to construct the library.

Up next: C code and PLT FFI for fast manipulation of SO(3,1) group elements (i.e. Lorentz transformations). Stay tuned.

Pearls Before Breakfast

Nothing to do with Scheme, but a really neat article here in the Washington Post. I'm really busy now, but I'm glad I took the time out to read this article. Found via the fascinating blog of Terrence Tao.