#include <stdio.h>

int main(){

  int i,j,k;

  double a[4][3],b[3][4],x[4][4],davidmeburk[3][4];

  void mmult(double a[][3],int ar,int ac,double b[][4],int bc, double x[][4]);


  FILE *fi,*fo;

  fi=fopen("data1.txt","r");

  for(i=0;i<4;i++){

    fscanf(fi,"%lf %lf %lf", &a[i][0],&a[i][1],&a[i][2]);

  }

  fclose(fi);

  fi=fopen("data2.txt","r");

  for(i=0;i<3;i++){

    fscanf(fi,"%lf %lf %lf %lf", &b[i][0],&b[i][1],&b[i][2],&b[i][3]);

  }

  fclose(fi);


  /* compute product using mmult function */

  mmult(a,4,3,b,4,x);

  
  fo=fopen("out1.txt","w");

  for(i=0;i<4;i++){

    fprintf(fo,"%f %f %f %f \n", x[i][0], x[i][1], x[i][2], x[i][3]);

  }

  fclose(fo);


  mmult(b,3,4,x,4,x,davidmeburk);

  fo=fopen("out2.txt","w");

  for(i=0;i<3;i++){

    fprintf(fo,"%f %f %f %f \n", davidmeburk[i][0], davidmeburk[i][1], davidmeburk[i][2], davidmeburk[i][3]);

  }

  fclose(fo);


  return 0;

}


void mmult(double a[][3],int ar,int ac,double b[][4],int bc, double x[][4]){

/* product of matrices */

/* a is ar x ac
   b is br x bc, where ac=br
*/

  int i,j,k;

  for(i=0;i<ar;i++){
    for(j=0;j<bc;j++){
      x[i][j]=0.0;
      for(k=0;k<ac;k++){

	x[i][j]=x[i][j]+a[i][k]*b[k][j];

      }
    }
    }
}