29 April 2007

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))))

5 comments:

Brad Lucier said...

This "Pure Scheme" code does not at all do the same thing as the C code used in the FFI. It would not be hard to take SO31.c and rewrite it using the primitives of Gambit-C or of the proposed R6RS and get just about the same speed and the same lack of flexibility and safety as you get from the C code. (It would not be hard, but it would be just as tedious as writing the C code in the first place.)

In fact, the only thing you rely on LAPACK for is dgemm to multiply 4 x 4 matrices. I will speculate that while dgemm is very fast for matrices of moderate size, there's enough overhead in dgemm that hand-writing a 4x4 matrix multiplier (which, given all the other low-level code you're written in SO31.c, would not be that much extra code) would be even *faster* than the LAPACK code.

But speculation is easy, and I decided not to take the time to rewrite SO31.c myself.

Will Farr said...

Hey Brad,

Thanks for the comments. Actually, at the moment I rely on LAPACK for computing matrix inverses, too, but I just realized (after your comment) that I don't need to do that at all. Given a set of parameters for a group element, negating the parameters, and reversing the product order of individual boost and rotation matrices will give the matrix inverse of the transformation matrix associated with those parameters. Something to do in version 1.2, I suppose.

As to your larger point, you're entirely correct---it would be basically just as much work to do this (at the same low level) in Gambit-C or any scheme which does a reasonable job of optimizing the new R6RS floating-point specific operations (but that probably won't be MzScheme). My point was, given that the amount of effort is a wash, at the end of the day, I like the library options which PLT Scheme gives me better than I like those of Gambit-C. So that's why I do things the way I do.

By the way, you're right that the code I posted looks very different from the SO31.plt package code (and, in fact, I didn't post the benchmark code at all---just the library code which would allow you to perform the benchmark). The code I posted is part of a larger collection of code, and provides a much richer library than SO31.plt. I just broke out the bare-bones stuff for SO31.plt because I didn't want people to have to get the entire library of my code in order to do something with Lorentz transformations.

Finally (sorry for the long comment, which you probably won't see at all anyway), I think you're probably right about specialized 4x4 matrix multipliers (definitely LAPACK isn't doing any blocking at 4x4), but I don't really see a need for that kind of speed yet. If I need it eventually, I'll keep it in mind.

Thanks,
Will

Brad Lucier said...

The main point I was trying to make was that the timings are of different algorithms implemented in different languages, so there's too much variation to learn much from them.

And of course I saw your comment! ;-)

Will Farr said...

Brad,

I agree with your point---I don't really think the comparison is very relevant. But, Jens asked to see the pure Scheme code, and I couldn't resist.... Who knows---maybe someone will find the pure-scheme Gaussian elimination matrix inversion routine useful (though perhaps I shouldn't encourage matrix inversion since it's such a crappy way to solve Mx=b).

Thanks for keeping up, not only with the blog, but with comments!

Will

Kathy said...

Hey there, Will! I searched for you and found your blog! I can't say that I understand any of it, but hello :)