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) {
	   /* The matrix is SINGULAR ! */
	   return;
	}

	/* 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];
}