/*---------------------------------------------------------------------- File : horner.c Contents: extended Horner method to evaluate a polynomial Author : Christian Borgelt History : 20.01.1998 file created 21.01.1998 comments added ----------------------------------------------------------------------*/ #include #include #ifndef VERSION #define VERSION 2 /* default: define second version */ #endif /*---------------------------------------------------------------------- Preprocessor Definitions ----------------------------------------------------------------------*/ #define MAXDEGREE 64 /* maximal degree of polynomial */ /*---------------------------------------------------------------------- Functions ----------------------------------------------------------------------*/ #if (VERSION == 1) /* if to compile first version */ void horner (int n, double c[], double x, double y[]) { /* --- extended Horner method 1 */ int i, k; /* loop variables */ for (k = n+1; --k >= 0; ) /* initialize the result vector to */ y[k] = c[n]; /* the coeff. of the highest exponent */ for (k = n; --k >= 0; ) { /* traverse the columns of the scheme */ y[0] = y[0] *x +c[k]; /* compute value in first line */ for (i = 1; i <= k; i++) /* traverse the lines of one column */ y[i] = y[i] *x +y[i-1]; /* and compute next column's values */ } i = 1; /* initialize i to 1! = 1 */ for (k = 2; k <= n; k++) /* traverse derivations 2 to n */ y[k] *= i *= k; /* and multiply them by i = k! */ } /* horner() */ #endif /*---------------------------------------------------------------------- This version of the extended Horner method to evaluate a polynomial of degree n with coefficients c[] and its derivations at a given point x computes the numbers in the evaluation scheme column by column. It starts by initializing the result vector y to the coefficient of the highest exponent in the polynomial. This yields the first column of the scheme. Then it computes the remaining columns one by one, each from top to bottom. To go right one column in a given line, one has to multiply by x and to add the value from the line of the new column that lies above the current one. For the first line, this is the line containing the coefficients. For any other line, this is the value computed in the step preceding the current one. Finally, after all columns are computed, the results in the vector y are multiplied by the factorials of the degrees of the derivations, as required to obtain the correct values. ----------------------------------------------------------------------*/ #if (VERSION == 2) /* if to compile second version */ void horner (int n, double c[], double x, double y[]) { /* --- extended Horner method 2 */ int i, k; /* loop variables */ for (k = n+1; --k >= 0; ) /* copy the coefficients */ y[k] = c[k]; /* to the resulkt vector */ for (k = 0; k < n; k++) /* traverse the lines of the scheme */ for (i = n; --i >= k; ) /* traverse the columns of one line */ y[i] += y[i+1] *x; /* and compute next column's value */ i = 1; /* initialize i to 1! = 1 */ for (k = 2; k <= n; k++) /* traverse derivations 2 to n */ y[k] *= i *= k; /* and multiply them by i = k! */ } /* horner() */ #endif /*---------------------------------------------------------------------- This version of the extended Horner method, which is slightly more elegant than the first, computes the numbers in the evaluation scheme line by line. In order not to destroy the coefficients, we first copy them to the result vector. We thus obtain the first line of the scheme. (Note that we could do entirely without the vector y, if we computed the results in the vector c.) Then we traverse the lines of the scheme. The first execution of the outer for-loop simply applies the normal Horner method, which yields the value of the polynomial. The second execution computes the value of the first derivation, and so on. Finally, after these nested loops, we multiply the results in the vector y by the factorials of the degrees of the derivations, just as in the first version. ----------------------------------------------------------------------*/ int main (int argc, char *argv[]) { /* --- main function */ int i; /* loop variable */ int n; /* number of coefficients (degree +1) */ double c[MAXDEGREE+1]; /* coefficients of polynomial */ double y[MAXDEGREE+1]; /* result vector for y^(i)(x) */ double x; /* point to compute values at */ FILE *file = stdin; /* file to read polynomial from */ if (argc != 3) { /* if wrong number of arguments given */ printf("usage: %s filename x\n", argv[0]); printf("evaluate a polynomial and its derivations " "at the point x\n"); printf("The coefficients of the polynomial are " "read from the file 'filename';\n"); printf("if 'filename' starts with a '-', " "input is read from stdin.\n"); return 0; /* print a usage message */ } /* and abort the function */ x = atof(argv[2]); /* convert x to a number */ if (argv[1][0] != '-') { /* if a file name is given, */ file = fopen(argv[1], "r"); /* open the file for reading */ if (!file) { /* if the file cannot be opened */ printf("%s: cannot open file %s\n", argv[0], argv[1]); return -1; /* print an error message */ } /* and abort the function */ } for (n = 0; n <= MAXDEGREE; n++) { if (fscanf(file, "%lf", c+n) != 1) break; /* read the coefficients of the */ } /* polynomial from the input file */ if (n < 1) { /* check degree of polynomial */ printf("%s: degree must be >= 0", argv[0]); return -1; } if (file != stdin) /* if input has been read */ fclose(file); /* from a file, close the file */ horner(n-1, c, x, y); /* call evaluation function and */ for (i = 0; i < n; i++) /* print the results from vector y */ printf("p^(%i)(%g) = %g\n", i, x, y[i]); return 0; /* return 'ok' */ } /* main() */ /*---------------------------------------------------------------------- For comments on the main function see the program `poly.c'. ----------------------------------------------------------------------*/