/*
   % The computations were worked out with the Symbolic
   % toolbox of Matlab:
   syms a00 a01 a02 a03 a10 a11 a12 a13 a20 a21 a22 a23 a30 a31 a32 a33
   syms b00 b01 b02 b03 b10 b11 b12 b13 b20 b21 b22 b23 b30 b31 b32 b33

   A = [a00 a01 a02 a03; a10 a11 a12 a13; a20 a21 a22 a23; a30 a31 a32 a33];
   B = [b00 b01 b02 b03; b10 b11 b12 b13; b20 b21 b22 b23; b30 b31 b32 b33];

   det(A)
   inv(A)*det(A)
   A*B
*/

/* Determinant of a 4x4 matrix */
double det44(double a[4][4])
{
    double dt;
    dt = a[0][0]*( a[1][1]*( a[2][2]*a[3][3]-a[2][3]*a[3][2])
                  +a[2][1]*(-a[1][2]*a[3][3]+a[1][3]*a[3][2])
                  +a[3][1]*( a[1][2]*a[2][3]-a[1][3]*a[2][2]))
        +a[1][0]*( a[0][1]*(-a[2][2]*a[3][3]+a[2][3]*a[3][2])
                  +a[2][1]*( a[0][2]*a[3][3]-a[0][3]*a[3][2])
                  +a[3][1]*(-a[0][2]*a[2][3]+a[0][3]*a[2][2]))
        +a[2][0]*( a[0][1]*( a[1][2]*a[3][3]-a[1][3]*a[3][2])
                  +a[1][1]*(-a[0][2]*a[3][3]+a[0][3]*a[3][2])
                  +a[3][1]*( a[0][2]*a[1][3]-a[0][3]*a[1][2]))
        +a[3][0]*( a[0][1]*(-a[1][2]*a[2][3]+a[1][3]*a[2][2])
                  +a[1][1]*( a[0][2]*a[2][3]-a[0][3]*a[2][2])
                  +a[2][1]*(-a[0][2]*a[1][3]+a[0][3]*a[1][2]));
    return(dt);
}

/* B = inv(A) */
void inv44(double a[4][4], double b[4][4])
{
    double idet;
    idet    = 1.0/det44(a);

    b[0][0] = idet*(a[1][1]*( a[2][2]*a[3][3]-a[2][3]*a[3][2])
                   +a[2][1]*(-a[1][2]*a[3][3]+a[1][3]*a[3][2])
                   +a[3][1]*( a[1][2]*a[2][3]-a[1][3]*a[2][2]));
    b[0][1] = idet*(a[0][1]*(-a[2][2]*a[3][3]+a[2][3]*a[3][2])
                   +a[2][1]*( a[0][2]*a[3][3]-a[0][3]*a[3][2])
                   +a[3][1]*(-a[0][2]*a[2][3]+a[0][3]*a[2][2]));
    b[0][2] = idet*(a[0][1]*( a[1][2]*a[3][3]-a[1][3]*a[3][2])
                   +a[1][1]*(-a[0][2]*a[3][3]+a[0][3]*a[3][2])
                   +a[3][1]*( a[0][2]*a[1][3]-a[0][3]*a[1][2]));
    b[0][3] = idet*(a[0][1]*(-a[1][2]*a[2][3]+a[1][3]*a[2][2])
                   +a[1][1]*( a[0][2]*a[2][3]-a[0][3]*a[2][2])
                   +a[2][1]*(-a[0][2]*a[1][3]+a[0][3]*a[1][2]));
    b[1][0] = idet*(a[1][0]*(-a[2][2]*a[3][3]+a[2][3]*a[3][2])
                   +a[2][0]*( a[1][2]*a[3][3]-a[1][3]*a[3][2])
                   +a[3][0]*(-a[1][2]*a[2][3]+a[1][3]*a[2][2]));
    b[1][1] = idet*(a[0][0]*( a[2][2]*a[3][3]-a[2][3]*a[3][2])
                   +a[2][0]*(-a[0][2]*a[3][3]+a[0][3]*a[3][2])
                   +a[3][0]*( a[0][2]*a[2][3]-a[0][3]*a[2][2]));
    b[1][2] = idet*(a[0][0]*(-a[1][2]*a[3][3]+a[1][3]*a[3][2])
                   +a[1][0]*( a[0][2]*a[3][3]-a[0][3]*a[3][2])
                   +a[3][0]*(-a[0][2]*a[1][3]+a[0][3]*a[1][2]));
    b[1][3] = idet*(a[0][0]*( a[1][2]*a[2][3]-a[1][3]*a[2][2])
                   +a[1][0]*(-a[0][2]*a[2][3]+a[0][3]*a[2][2])
                   +a[2][0]*( a[0][2]*a[1][3]-a[0][3]*a[1][2]));
    b[2][0] = idet*(a[1][0]*( a[2][1]*a[3][3]-a[2][3]*a[3][1])
                   +a[2][0]*(-a[1][1]*a[3][3]+a[1][3]*a[3][1])
                   +a[3][0]*( a[1][1]*a[2][3]-a[1][3]*a[2][1]));
    b[2][1] = idet*(a[0][0]*(-a[2][1]*a[3][3]+a[2][3]*a[3][1])
                   +a[2][0]*( a[0][1]*a[3][3]-a[0][3]*a[3][1])
                   +a[3][0]*(-a[0][1]*a[2][3]+a[0][3]*a[2][1]));
    b[2][2] = idet*(a[0][0]*( a[1][1]*a[3][3]-a[1][3]*a[3][1])
                   +a[1][0]*(-a[0][1]*a[3][3]+a[0][3]*a[3][1])
                   +a[3][0]*( a[0][1]*a[1][3]-a[0][3]*a[1][1]));
    b[2][3] = idet*(a[0][0]*(-a[1][1]*a[2][3]+a[1][3]*a[2][1])
                   +a[1][0]*( a[0][1]*a[2][3]-a[0][3]*a[2][1])
                   +a[2][0]*(-a[0][1]*a[1][3]+a[0][3]*a[1][1]));
    b[3][0] = idet*(a[1][0]*(-a[2][1]*a[3][2]+a[2][2]*a[3][1])
                   +a[2][0]*( a[1][1]*a[3][2]-a[1][2]*a[3][1])
                   +a[3][0]*(-a[1][1]*a[2][2]+a[1][2]*a[2][1]));
    b[3][1] = idet*(a[0][0]*( a[2][1]*a[3][2]-a[2][2]*a[3][1])
                   +a[2][0]*(-a[0][1]*a[3][2]+a[0][2]*a[3][1])
                   +a[3][0]*( a[0][1]*a[2][2]-a[0][2]*a[2][1]));
    b[3][2] = idet*(a[0][0]*(-a[1][1]*a[3][2]+a[1][2]*a[3][1])
                   +a[1][0]*( a[0][1]*a[3][2]-a[0][2]*a[3][1])
                   +a[3][0]*(-a[0][1]*a[1][2]+a[0][2]*a[1][1]));
    b[3][3] = idet*(a[0][0]*( a[1][1]*a[2][2]-a[1][2]*a[2][1])
                   +a[1][0]*(-a[0][1]*a[2][2]+a[0][2]*a[2][1])
                   +a[2][0]*( a[0][1]*a[1][2]-a[0][2]*a[1][1]));
}

/* C = A*B */
void mul44(double a[4][4], double b[4][4], double c[4][4])
{
    c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0]+a[0][3]*b[3][0];
    c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1]+a[0][3]*b[3][1];
    c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2]+a[0][3]*b[3][2];
    c[0][3] = a[0][0]*b[0][3]+a[0][1]*b[1][3]+a[0][2]*b[2][3]+a[0][3]*b[3][3];
    c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0]+a[1][3]*b[3][0];
    c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1]+a[1][3]*b[3][1];
    c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2]+a[1][3]*b[3][2];
    c[1][3] = a[1][0]*b[0][3]+a[1][1]*b[1][3]+a[1][2]*b[2][3]+a[1][3]*b[3][3];
    c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0]+a[2][3]*b[3][0];
    c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1]+a[2][3]*b[3][1];
    c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2]+a[2][3]*b[3][2];
    c[2][3] = a[2][0]*b[0][3]+a[2][1]*b[1][3]+a[2][2]*b[2][3]+a[2][3]*b[3][3];
    c[3][0] = a[3][0]*b[0][0]+a[3][1]*b[1][0]+a[3][2]*b[2][0]+a[3][3]*b[3][0];
    c[3][1] = a[3][0]*b[0][1]+a[3][1]*b[1][1]+a[3][2]*b[2][1]+a[3][3]*b[3][1];
    c[3][2] = a[3][0]*b[0][2]+a[3][1]*b[1][2]+a[3][2]*b[2][2]+a[3][3]*b[3][2];
    c[3][3] = a[3][0]*b[0][3]+a[3][1]*b[1][3]+a[3][2]*b[2][3]+a[3][3]*b[3][3];
}

#include <stdio.h>
#include <strings.h>
void printmat44(double a[4][4])
{
    int i, j;
    (void)printf("\nmatrix = [\n");
    for(i=0; i<4; i++)
        (void)printf("%12.8g %12.8g %12.8g %12.8g\n", a[i][0],a[i][1],a[i][2],a[i][3]);
    (void)printf("]\n\n");
}

int scanmat44(char c[], double a[4][4])
{
    return(sscanf(c,"%lg,%lg,%lg,%lg;%lg,%lg,%lg,%lg;%lg,%lg,%lg,%lg;%lg,%lg,%lg,%lg",
        &a[0][0],&a[0][1],&a[0][2],&a[0][3],
        &a[1][0],&a[1][1],&a[1][2],&a[1][3],
        &a[2][0],&a[2][1],&a[2][2],&a[2][3],
        &a[3][0],&a[3][1],&a[3][2],&a[3][3]));
}
main(int argc, char *argv[])
{
    double a[4][4], b[4][4], c[4][4];
    if (argc==1)
    {
        char Am[] = "3.0,-2.1,0.1,-1.0;0.9,1.3,2.0,-1.0;1.9,-0.9,1.1,-0.1;0.8,-0.1,0.1,1.2";
        char Bm[] = "3.2, 2.0,0.3,1.2;-0.9,2.3,1.0,-1.2;1.2,-0.8,1.4,-0.3;0.1,-0.7,0.3,1.1";
        (void)fprintf(stderr,"%s\n%s '%s' '%s'\n\n%s\nA = [%s];\nB=[%s];\n%s\n%s\n%s\n%s\n\n",
            "Try the following (for example):",
            argv[0],Am, Bm,
            "and compare with this in Matlab:",
            Am, Bm,"disp(A);","disp(inv(A));","disp(B);","disp(A*B);");

    }
    if (argc>1)
    {
        scanmat44(argv[1],a);
        printmat44(a);
        inv44(a,c);
        printmat44(c);
    }
    if (argc>2)
    {
        scanmat44(argv[2],b);
        printmat44(b);
        mul44(a,b,c);
        printmat44(c);
    }
    return(0);
}