// So far this code computes a random orthogonal matrix.
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define M 4
typedef double r;
r m[M][M];
const int sig = 0; // 0 for Euclid, 1 for Minkowski
static void twist(int i, int j, r theta){
  void tw(r a, r b, r c){int k=M;
    while(k--){r t = a*m[k][i] + b*m[k][j]; m[k][j] = c*m[k][i] + a*m[k][j]; m[k][i] = t;}}
  if (((sig>>i)^(sig>>j))&1) {r s=sinh(theta), c=cosh(theta); tw(c, s, s);}
  else {r s=sin(theta), c=cos(theta); tw(c, s, -s);}}
  
static void to(char * cl){int j=M; while(j--){int k=j+1; while(k--){int n=M; r s = 0;
  while(n--) s += ((sig>>n)&1?-1.:1.)*m[j][n]*m[k][n];
  {r e = s - ((sig>>j)&1?-1.:1.)*(j==k);
    if(fabs(e) > 1e-12) printf("err: %s %d %d %e %e\n", cl, j, k, e, s);}}}}

static void pm(char * y){int j; for(j=0;j<M;++j){int k;
  for(k=0;k<M;++k) printf(" %18.15f", m[j][k]); printf("\n");} printf("%s\n", y);}

static void stir(){int j = 1024; while(j--){int k=random()%M, l = random()%M;
  if(k!=l) {r v = 1?(0.125 - (float)random()*0x1.p-33):acosh(1.1); twist(k, l, v);}}}

int main(){{int j=M; while(j--){int k=M; while(k--) m[j][k] = (k==j);}}
to("gli"); stir(); to("glo"); pm ("om");
return 0;}

// AηAtη = I
// (define eta '((-1 0 0 0) (0 1 0 0) (0 0 1 0) (0 0 0 1)))

//(define (z x)(let ((x (quad x))) (mm (mm x eta) (mm (transpose x) eta))))
//