#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void make_rnd(double a[][3]);
void make_rotx(double a[][3], double x);
void make_rotz(double a[][3], double z);
void transpose_mat(double x[][3], double a[][3]);
void inverse_mat(double x[][3]);
void mat_mul(double x[][3], double a[][3], double b[][3]);
void print_mul(char *s, double a[][3]);

int main(void)
{
  double x, z;
  char buf[128];
  
  double a[3][3], b[3][3], u[3][3], t[3][3], s[3][3];
  
  printf("x(rad) = ");
  scanf("%lf", &x);
  printf("z(rad) = ");
  scanf("%lf", &z);
  gets(buf);
  
  make_rotx(a, x);
  make_rotz(b, z);

  mat_mul(u, a, b);
  print_mul("直交行列", u);
  
  transpose_mat(t, u);
  print_mul("転置行列", t);
  
  mat_mul(s, u, t);
  print_mul("直交行列と転置行列との積（＝単位行列）", s);
  
  inverse_mat(u);
  print_mul("直交行列の逆行列(＝転置行列）", u);
  
  return 0;
}

void make_rnd(double a[][3])
{
  int i, j;
  
  srand((unsigned)time(NULL));
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      a[i][j] = (double)rand() / RAND_MAX;
}

/*
 * (cosX -sinX 0)
 * (sinX cosX 0) 
 * (0 0 1)
 */

void make_rotx(double a[][3], double x)
{
  a[0][0] = a[1][1] = cos(x);
  a[0][1] = -sin(x);
  a[1][0] = sin(x);
  a[0][2] = a[1][2] = a[2][0] = a[2][1] = 0;
  a[2][2] = 1;
}

/*
 * (1 0 0 )
 * (0 cosZ -sinZ)
 * (0 sinZ cosZ)
 */

void make_rotz(double a[][3], double z)
{
  a[1][1] = a[2][2] = cos(z);
  a[1][2] = -sin(z);
  a[2][1] = sin(z);
  a[0][1] = a[0][2] = a[1][0] = a[2][0] = 0;
  a[0][0] = 1;
}

void mat_mul(double x[][3], double a[][3], double b[][3])
{
  int i, j, k;
  double s;
  
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      s = 0;
      for (k = 0; k < 3; k++)
        s += a[i][k] * b[k][j];
      x[i][j] = s;
    }
  }
}

void transpose_mat(double x[][3], double a[][3])
{
  int i, j;
  
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
        x[i][j] = a[j][i];
}

void inverse_mat(double x[][3])
{
  int i, j, k;
  double p, q;
  
  for (k = 0; k < 3; k++) {
    p = x[k][k];
    x[k][k] = 1;
    for (j = 0; j < 3; j++) x[k][j] /= p;
    for (i = 0; i < 3; i++)
      if (i != k) {
        q = x[i][k];
        x[i][k] = 0;
        for (j = 0; j < 3; j++) x[i][j] -= q * x[k][j];
      }
  }
}

void print_mul(char *s, double a[][3])
{
  int i, j;
  
  puts(s);
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) printf("%10.6f ", a[i][j]);
    putchar('\n');
  }
}