; f below compute metric tensor.
(define f
  (lambda (input) (let z ((r (cdr input)))
   (if (null? r) '() (cons
                  (let y ((s (cdar r))(v (cdar input)))
   (if (null? s) '() (cons (/ (- (+ (caar r) (car v)) (car s)) 2)
     (y (cdr s)(cdr v)))))
     (z (cdr r)))))))

(define exit (let* ((s (cons 0 0))
   (z (call-with-current-continuation (lambda (e) (cons s e)))))
      (if (and (pair? z) (eq? (car z) s)) (cdr z) exit)))

(define mer (lambda (x) (let ((r (call-with-current-continuation
  (lambda(d) (cons #f (lambda(y) (d (cons #t y))))))))
   (if (car r) (begin (write (list "Exception:" x (cdr r))) (newline) (exit 0))
   (cdr r)))))

(define cc (mer "Singular:"))

(define (trnxx zer zer?) (lambda (x) (if (null? x) '()
  (let ((z ((trnxx zer zer?) (cdr x)))) (let m ((u (car x))(v z))
    (if (null? u) (map (lambda (q) (if (null? q) '() (cons zer q))) v)
     (if (null? v) (let y ((q u)) (if (null? q) '() (let ((p (y (cdr q))))
         (if (and (null? p) (zer? (car q))) '()
           (cons (if (zer? (car q)) '() (list (car q))) p)))))
       (cons (cons (car u)(car v)) (m (cdr u)(cdr v))))))))))

(define matpak (lambda (zero zer? one fa fs fm fi) (letrec (
   (tf (lambda (x zp) (if (null? x) (zp (list one))
     (if (zer? (caar x)) (let ((z (tf (cdr x) zp)))
     (cons (car z) (cons (car x) (cdr z)))) x))))
   (tr (trnxx zero zer?))
   (sm (lambda (s a)(map (lambda (w) (fm s w)) a)))
   (fms (lambda (a s b) (if (null? a) (sm (fs zero s) b) (if (null? b) a
     (cons (fs (car a) (fm s (car b))) (fms (cdr a) s (cdr b)))))))
   (deter (lambda (a) (letrec ((p #f)
     (tf (lambda (x) (if (null? x) (list one)
     (if (zer? (caar x)) (let ((z (tf (cdr x) zp)))
     (set! p (not p))
     (cons (car z) (cons (car x) (cdr z)))) x)))))
     (let inx ((d one)(a a))
          (if (null? a) (if p (fs zero d) d) (let* (
            (A (tf a))
            (i (fi (caar A)))
            (b (map (lambda (z) (fm z i)) (cdar A))))
         (inx (fm (caar A) d) (map (lambda (x w) (fms x w b))
          (map cdr (cdr A)) (map car (cdr A))))))))))
   (inv (lambda (a nxp) (let ol ((ut (let inx (
      (a (let pad ((x a)(e (list one)))
         (if (null? x) '() (cons (let ap ((z (car x))(ln a))
           (if (null? ln) e (if (null? z) (cons zero (ap z (cdr ln)))
             (cons (car z)(ap (cdr z)(cdr ln))))))
          (pad (cdr x) (cons zero e))))))
      (np nxp))
          (if (null? a) '() (let* (
            (A (tf a np))
            (i (fi (caar A)))
            (b (map (lambda (z) (fm z i)) (cdar A))))
         (cons b (inx (map (lambda (x w) (fms x w b))
          (map cdr (cdr A)) (map car (cdr A)))
          (lambda (w) (np (cons (fs zero (ip w b)) w))))))))))
      (if (null? ut) '() (cons
        (let eg ((top (car ut))(bod (cdr ut))) (if (null? bod) top
          (eg (fms (cdr top) (car top) (car bod))(cdr bod))))
      (ol (cdr ut)))))))
   (ip (lambda (x y)(if (or (null? x) (null? y)) zero (fa (fm (car x)(car y))
        (ip (cdr x)(cdr y))))))
   (mp (lambda (a b)(let ((b (tr b)))
       (map (lambda(x) (map (lambda (y) (ip x y)) b)) a)))))

(list mp inv ip tr deter))))

(define gr (lambda (f) (lambda z (let w ((z z))
   (if (null? (car z)) '() (cons (apply f (map car z)) (w (map cdr z))))))))
(define v+ (gr +))
(define v- (gr -))

(define (sum x)(if (null? x) 0 (+ (car x)(sum (cdr x)))))
(define vsum (gr sum))
(define pi 3.1415926535897932385)

(define (pl x)(if (pair? x)(begin (display (car x)) (newline) (pl (cdr x)))))

(define wk (matpak 0 zero? 1 + - * (lambda (x) (/ x))))
(define (inv x) ((cadr wk) x cc))

(define lt (lambda (f) (lambda (x) (let z ((x x))
   (if (null? x) '() (let w ((v (cdr x))) (if (null? v) (z (cdr x))
      (cons (f (car x) (car v)) (w (cdr v))))))))))

(define (ip x y)(let lx ((xr x)(g g))(if (null? xr) 0 (+ 
                (let ly ((yr y)(g (car g))) (if (null? yr) 0 (+
                (* (car g)(car yr)(car xr))(ly (cdr yr)(cdr g)))))
            (lx (cdr xr)(cdr g))))))

(define (nrm x) (let* ((y (ip x x))(z (/ (sqrt y))))
   (let q ((x x)) (if (null? x) '() (cons (* z (car x)) (q (cdr x)))))))

(define (fact n)(if (zero? n) 1 (if (= n -1/2) (sqrt pi) (* n (- n 1)))))
(define (cosAngle a b) (/ (ip a b) (sqrt (* (ip a a) (ip b b)))))

(define (rg n)(if (= n 1) '((0)) (let ((sv (rg (- n 1))))
   (cons (cons 0 (cons 1 (cdar sv)))
         (let x ((q sv)) (if (null? q) '() (cons (cons 1 (car q)) (x (cdr q)))))))))

;(define reg '((0 1 1 1) (1 0 2 2) (1 2 0 2) (1 2 2 0))) ; x, y, z > 0 and x+y+z < 1
(define reg '((0 1 2 2) (1 0 1 1) (2 1 0 2) (2 1 2 0))) ; x, y, z > 0 and x+y+z < 1
;(define reg '((0 1 1) (1 0 1) (1 1 0))) ; equilateral triangle
;(define reg '((0 1 1 1) (1 0 1 1) (1 1 0 1) (1 1 1 0))) ; Regular tetrahedron
;(define reg '((0 1 1 1 1) (1 0 1 1 1) (1 1 0 1 1) (1 1 1 0 1) (1 1 1 1 0))) ; regular 4D simplex
;(define reg (rg 24)) ; regular 24-simplex
(define g (f reg))
(define gi (inv g))
(define normals (cons (v- (vsum gi)) gi))
(define unitnrms ((gr nrm) normals))

(define trih (let ((a (car unitnrms))(b (cadr unitnrms))(c (caddr unitnrms)))
   (let ((A (acos (cosAngle b c))) (B (acos (cosAngle c a))) (C (acos (cosAngle a b))))
     (- (* 2 pi) (+ A B C)))))

(pl (list 
(cons "Simplex edge lengths squared:" reg)
(cons "Volume: " (/ (sqrt ((car (cddddr wk)) g)) (fact (length g))))
(cons "Covariant Metric tensor: " g)
(cons "Contravariant metric tensor: " gi)
(cons "Cosines of dihedral angles: " (v- ((lt cosAngle) normals)))
(cons "First trihedral angle: " trih)))