/*---------------------------------------------------------------------- File : poly.c Contents: evaluate a polynomial at a given point x Author : Christian Borgelt History : 13.12.1997 file created ----------------------------------------------------------------------*/ #include #include #ifndef VERSION #define VERSION 3 /* default: compile third version */ #endif /*---------------------------------------------------------------------- Preprocessor Definitions ----------------------------------------------------------------------*/ #define MAXDEGREE 64 /* maximal degree of polynomial */ /*---------------------------------------------------------------------- Functions ----------------------------------------------------------------------*/ #if (VERSION == 1) /* if to compile first version */ double power (double x, int n) { /* --- compute x^n, n >= 1 */ if (n <= 1) return x; /* x^1 can be computed directly */ if (n & 1) return power(x *x, n >> 1) *x; else return power(x *x, n >> 1); } /* power() */ /* otherwise apply recursive formula */ /*--------------------------------------------------------------------*/ double poly (int n, double c[], double x) { /* --- eval. the degree n polynomial */ /* with coefficients c at point x */ double y = 0; /* result (value of polynomial at x) */ while (n > 0) { /* while not all coeff. processed, */ y += c[n] *power(x, n); /* add next term of the polynomial */ n--; /* and reduce the exponent by one */ } /* return the computed value plus */ return y +c[0]; /* the last coefficient * x^0 */ } /* poly() */ #endif /*---------------------------------------------------------------------- The above function computes the value of the polynomial by successively summing the terms c_k * x^k of the polynomial, where x^k is computed using a `power' function. (See the program `power.c' for an explanation of the `power' function.) It is obvious that the loop has to be executed n times to compute the result. In each execution the power function is called, which has the time complexity O(log(n)). Hence the overall time complexity of this function is O(n log(n)). Note that it would be even worse, if the `power' function used a naive approach, i.e. computed x^n as follows: double power (double x, int n) { double y = x; while (--n > 0) y *= x; return y; } This implementation of the `power' function obviously has a time complexity of O(n) and thus the overall time complexity of the above 'poly' function would worsen to O(n^2). ----------------------------------------------------------------------*/ #if (VERSION == 2) /* if to compile second version */ double poly (int n, double c[], double x) { /* --- eval. the degree n polynomial */ /* with coefficients c at point x */ int k = 0; /* loop variable, coefficient index */ double x_k = 1; /* to traverse the values x^k, k >= 0 */ double y; /* result (value of polynomial at x) */ y = c[0]; /* start with last coefficient */ while (++k <= n) { /* while not all coeff. processed, */ x_k *= x; /* compute x^k from x^(k-1) and */ y += c[k] *x_k; /* add next term of the polynomial */ } /* to the result computed up to now */ return y; /* return the computed value */ } /* poly() */ #endif /*---------------------------------------------------------------------- This function also computes the value of the polynomial by successively summing the terms c_k * x^k of the polynomial, but it does so without using a `power' function. It computes the values of the x^k by successive multiplications by x. This brings the time complexity down to O(n), since the costs of the loop body are constant now. ----------------------------------------------------------------------*/ #if (VERSION == 3) /* if to compile third version */ double poly (int n, double c[], double x) { /* --- eval. the degree n polynomial */ /* with coefficients c at point x */ double y; /* result (value of polynomial at x) */ y = c[n]; /* start with highest coefficient */ while (--n >= 0) /* while not all coeff. processed, */ y = y *x +c[n]; /* compute step of Horner's algorithm */ return y; /* return the computed value */ } /* poly() */ #endif /*---------------------------------------------------------------------- This function uses Horner's algorithm to compute the value of the polynomial. It also has a time complexity of O(n), since the loop is executed n times and the cost of the loop body is constant. Nevertheless, this function is faster than the second version since in this version only one multiplication has to be carried out in the loop body in contrast to two multiplications in the second version. ----------------------------------------------------------------------*/ int main (int argc, char *argv[]) { /* --- main function */ int n; /* degree of polynomial */ double c[MAXDEGREE+1]; /* coefficients of polynomial */ double x; /* point to eval. the polynomial 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 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 */ printf("%g\n", poly(n-1, c, x)); /* compute and print the result */ return 0; /* return 'ok' */ } /* main() */ /*---------------------------------------------------------------------- See the program `mps.c' for an explanation of how to read input from a file. We chose to read the coefficients of the polynomial from a file in order to free the user from having to type in the coefficients over and over again, if he/she needs the value of the polynomial for several values of x. ----------------------------------------------------------------------*/