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