aptdec/reg.c

/* ---------------------------------------------------------------------------

 Polynomial regression, freely adapted from :
 
 NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995
 Algorithm translated to  C  by: Dr. Norman Fahrer
 NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
 Prentice Hall, International Editions:   ISBN 0-13-625047-5
 This free software is compliments of the author.
 E-mail address:       in%"mathews@fullerton.edu"
*/  
  
#include<math.h>
  
#define DMAX  5    /* Maximum degree of polynomial */
#define NMAX  10    /* Maximum number of points     */

static void FactPiv (int N, double A[DMAX][DMAX], double B[], double Cf[]);


void
polyreg (const int M, const int N, const double X[], const double Y[],
	 double C[]) 
{
  
int R, K, J;			/* Loop counters     */
  
double A[DMAX][DMAX];	/* A                 */
  
double B[DMAX];
  
double P[2 * DMAX + 1];
  
double x, y;
  
double p;
  
    /* Zero the array */ 
    
for (R = 0; R < M + 1; R++)
    B[R] = 0;
  
    /* Compute the column vector */ 
    
for (K = 0; K < N; K++)
    
 {
    
y = Y[K];
    
x = X[K];
    
p = 1.0;
    

for (R = 0; R < M + 1; R++)
      
 {
      
B[R] += y * p;
      
p = p * x;
      
}
    
}
  
    /* Zero the array */ 
    
for (J = 1; J <= 2 * M; J++)
    P[J] = 0;
  

P[0] = N;
  
    /* Compute the sum of powers of x_(K-1) */ 
    
for (K = 0; K < N; K++)
    
 {
    
x = X[K];
    
p = X[K];
    

for (J = 1; J <= 2 * M; J++)
      
 {
      
P[J] += p;
      
p = p * x;
      
}
    
}
  
    /* Determine the matrix entries */ 
    
for (R = 0; R < M + 1; R++)
    
 {
    
for (K = 0; K < M + 1; K++)
      A[R][K] = P[R + K];
    
}
  
    /* Solve the linear system of M + 1 equations : A*C = B 
       for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */ 
    FactPiv (M + 1, A, B, C);


}				/* end main */


/*--------------------------------------------------------*/ 
  
static void
FactPiv (int N, double A[DMAX][DMAX], double B[], double Cf[]) 
{
  

int K, P, C, J;		/* Loop counters               */
  
int Row[NMAX];		/* Field with row-number       */
  
double X[DMAX], Y[DMAX];
  
double SUM, DET = 1.0;
  
int T;
  
    /* Initialize the pointer vector */ 
    
for (J = 0; J < N; J++)
    Row[J] = J;
  
    /* Start LU factorization */ 
    
for (P = 0; P < N - 1; P++) {
    
      /* Find pivot element */ 
      
for (K = P + 1; K < N; K++) {
      
if (fabs (A[Row[K]][P]) > fabs (A[Row[P]][P])) {
	
	  /* Switch the index for the p-1 th pivot row if necessary */ 
	  T = Row[P];
	
Row[P] = Row[K];
	
Row[K] = T;
	
DET = -DET;
      
}
    

}				/* End of simulated row interchange */
    

if (A[Row[P]][P] == 0) {
      
printf ("The matrix is SINGULAR !\n");
      
printf ("Cannot use algorithm --> exit\n");
      
exit (1);
    
}
    
      /* Multiply the diagonal elements */ 
      
DET = DET * A[Row[P]][P];
    
      /* Form multiplier */ 
      
for (K = P + 1; K < N; K++) {
      
A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P];
      
	/* Eliminate X_(p-1) */ 
	
for (C = P + 1; C < N + 1; C++) {
	
A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C];
      
}
    
}
  

}				/* End of  L*U  factorization routine */
  

DET = DET * A[Row[N - 1]][N - 1];
  
    /* Start the forward substitution */ 
    
for (K = 0; K < N; K++)
    Y[K] = B[K];
  

Y[0] = B[Row[0]];
  
for (K = 1; K < N; K++) {
    
SUM = 0;
    
for (C = 0; C <= K - 1; C++)
      SUM += A[Row[K]][C] * Y[C];
    
Y[K] = B[Row[K]] - SUM;
  
}
  
    /* Start the back substitution */ 
    
X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1];
  

for (K = N - 2; K >= 0; K--) {
    
SUM = 0;
    
for (C = K + 1; C < N; C++) {
      
SUM += A[Row[K]][C] * X[C];
    
}
    

X[K] = (Y[K] - SUM) / A[Row[K]][K];
  

}				/* End of back substitution */
  
    /* Output */ 
    for (K = 0; K < N; K++)
    
Cf[K] = X[K];

}