// This presents a mathmatical mystery. It says such things as ABCA'DB'C'D' = I
// where only BC and AD commute. (A' is the inverse of A)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define d 4
#define dep 8
typedef double r;
r c, s, th = 1;
char hx[dep][2];
static void to(char * q, r m[d][d]){int j=d; while(j--){int k=j+1; while(k--){
  int n=d; r s = 0;
  while(n--) s += m[j][n]*m[k][n]; {r e = s - (j==k);
    if(fabs(e) > 1.e-14) printf("Grock %s %d %d %e\n", q, j, k, e);}}}}
int cnt=0;
r o[d][d]; // the main matrix
r ID[d][d]; // comparand
static void id(r w[d][d]){int j=d; while(j--){int k=d; while(k--) w[j][k] = j==k;}}
static r sq(r x){return x*x;}
static void twist(int f, int g){int n=d; while(n--){
  r t = o[f][n], b = o[g][n];
  o[f][n] = c*t+s*b; o[g][n]=-s*t+c*b;} ++cnt;}
static void pt(char * m, r t[d][d]){int j, k; for(j=0; j<d; ++j) {
  for(k=0; k<d; ++k) printf("%10.6f ", t[j][k]); printf("%s\n", m);}}
static void rec(int j, int k, int mn, int D){ // j!=k && mn has just two bits.
  {int jk = (1<<j)|(1<<k);
  if((jk&mn) && (jk != mn) || (jk<mn)) {r hs[2][d];
    hx[D-1][0] = j; hx[D-1][1] = k;
    {int x=d; while(x--){hs[0][x] = o[j][x]; hs[1][x] = o[k][x];}}
    twist(j, k);
    {int n=16; r s=0; while(n--) {s += sq(o[0][n] - ID[0][n]);
     if(s > .046) goto xx;}
     pt("cls", o); to("cl", o); printf("s = %e c=%d\n", s, cnt);
     {int j=dep; while(j-->=D) printf("%d:%d, ", hx[j][0], hx[j][1]); printf("\n\n");} xx:;}
    {int a=d; while(a--){int b=d; while(b--) if(a!=b && D) rec(a, b, jk, D-1);}}
    {int x=d; while(x--){o[j][x] = hs[0][x]; o[k][x] = hs[1][x];}}}}}

int main(){s = sin(th); c = cos(th); id(o); id(ID);
  {int j=d; while(j--){int k=d; while(k--) if(j!=k) rec(j, k, 3<<d, dep);}}
printf("count=%d\n", cnt); return 0;}