Linear Algebra and the C Language/a0jz

/* ------------------------------------ */
/*  Save as :   c00f.c                  */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RA  R5
#define   CA  C5
#define   CXY C2
/* ------------------------------------ */
int main(void)
{
double tA[RA*CA]={
/*  x**4   x**3   x**2    x**1    x**0     */
  +1,     +1,     +1,     +1,     +1,       
 +16,     +8,     +4,     +2,     +1,        
 +81,    +27,     +9,     +3,     +1,        
+256,    +64,    +16,     +4,     +1,         
+625,   +125,    +25,     +5,     +1,       
};

double tb[RA*C1]={
/*    y  */
   -2,
   -2,
    3,
   -9,
    4,
};

double xy[RA*CXY] ={
  1,   -2,
  2,   -2,
  3,    3,
  4,   -9,
  5,    4,      };
  
double **XY   = ca_A_mR(xy,i_mR(RA,CXY));
double **A    = ca_A_mR(tA,i_mR(RA,CA));
double **b    = ca_A_mR(tb,i_mR(RA,C1));
double **Inv  =            i_mR(CA,RA);          
double **Invb =            i_mR(CA,C1);          

  clrscrn();
  printf(" Fitting a linear Curve to Data :\n\n");
  printf("    x     y");
  p_mR(XY,S5,P0,C6);
  printf(" A :\n  x**4    x**3    x**2    x**1    x**0");
  p_mR(A,S7,P2,C7);
  printf(" b :\n   y ");
  p_mR(b,S7,P2,C7);
  stop();
  
  clrscrn();   
  printf(" Inv : ");
  invgj_mR(A,Inv); 
  pE_mR(Inv,S12,P4,C10); 
  
  printf(" x = Inv * b ");   
  mul_mR(Inv,b,Invb); 
  p_mR(Invb,S10,P4,C10);
  printf(" The Quartic equation Curve to Data : \n\n"
         "  s = %+.3ft**4 %+.3ft**3 %+.3f*t**2 %+.3ft %+.3f\n\n"
            ,Invb[R1][C1],Invb[R2][C1],Invb[R3][C1],
             Invb[R4][C1],Invb[R5][C1]);        
  stop();  

  f_mR(XY);
  f_mR(b);  
  f_mR(A); 
  f_mR(Inv);
  f_mR(Invb); 

  return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */

  Let's calculate the coefficients of a polynomial.
 
              y =  ax**4 + bx**3 + cx**2 + dx + e        
  
         Which passes through these five points.    
          
       x[1],  y[1] 
       x[2],  y[2] 
       x[3],  y[3] 
       x[4],  y[4] 
       x[5],  y[5] 
       
  Using the points we obtain the matrix:

      x**4         x**3        x**2      x**1      x**0     y

     x[1]**4     x[1]**3     x[1]**2   x[1]**1   x[1]**0   y[1]
     x[2]**4     x[2]**3     x[2]**2   x[2]**1   x[2]**0   y[2]
     x[3]**4     x[3]**3     x[3]**2   x[3]**1   x[3]**0   y[3]
     x[4]**4     x[4]**3     x[4]**2   x[4]**1   x[4]**0   y[4]
     x[5]**4     x[5]**3     x[5]**2   x[5]**1   x[5]**0   y[5]

  That we can write:

       x**4          x**3        x**2     x     1    y

      x[1]**4       x[1]**3    x[1]**2   x[1]   1   y[1]
      x[2]**4       x[2]**3    x[2]**2   x[2]   1   y[2]
      x[3]**4       x[3]**3    x[3]**2   x[3]   1   y[3]
      x[4]**4       x[4]**3    x[4]**2   x[4]   1   y[4]
      x[5]**4       x[5]**3    x[5]**2   x[5]   1   y[5]   

     Let's use the invgj_mR() function to solve
     the system that will give us the coefficients a, b, c, d, e

 Fitting a linear Curve to Data :

    x     y
   +1    -2 
   +2    -2 
   +3    +3 
   +4    -9 
   +5    +4 

 A :
  x**4    x**3    x**2    x**1    x**0
  +1.00   +1.00   +1.00   +1.00   +1.00 
 +16.00   +8.00   +4.00   +2.00   +1.00 
 +81.00  +27.00   +9.00   +3.00   +1.00 
+256.00  +64.00  +16.00   +4.00   +1.00 
+625.00 +125.00  +25.00   +5.00   +1.00 

 b :
   y 
  -2.00 
  -2.00 
  +3.00 
  -9.00 
  +4.00 

 Press return to continue. 


 Inv : 
 +4.1667e-02  -1.6667e-01  +2.5000e-01  -1.6667e-01  +4.1667e-02 
 -5.8333e-01  +2.1667e+00  -3.0000e+00  +1.8333e+00  -4.1667e-01 
 +2.9583e+00  -9.8333e+00  +1.2250e+01  -6.8333e+00  +1.4583e+00 
 -6.4167e+00  +1.7833e+01  -1.9500e+01  +1.0167e+01  -2.0833e+00 
 +5.0000e+00  -1.0000e+01  +1.0000e+01  -5.0000e+00  +1.0000e+00 

 x = Inv * b 
   +2.6667 
  -30.3333 
 +117.8333 
 -181.1667 
  +89.0000 

 The Quartic equation Curve to Data : 

  s = +2.667t**4 -30.333t**3 +117.833*t**2 -181.167t +89.000

 Press return to continue.