/**************************************************************** Treibergs Feb. 24, 2006 HW 3.1 Function subprograms mycos and fastcos ****************************************************************/ /* Maclaurin's series for cosine is given by S = a_0 + a_1 + a_2 + a_3 + a_4 + ... where a_0 = 1 a_1 = -x^2/2 a_2 = x^4/(2*3*4) = -x^2*[-x^2/2]/(3*4) a_3 = -x^6/(2*3*4*5*6) = -x^2*[x^4/(2*3*4)]/(5*6) a_4 = x^8/(2*3*4*5*6*7*8) = -x^2*[-x^6/(2*3*4*5*6)]/(7*8) ... Let y = - x*x. Thus a_i satisfies the recursion a_1 = y/2 and a_i = -x^2*a_{i-2}/(i*(i+1)) for i=3,5,7,.... The partial sums S_n = a_1 + ... + a_n corresponds to the Taylor polynomias p_{2n} = p_{2n+1} = S_n of degrees 2n and 2n+1 since the series consists only of even terms. Thus the remainder |cos(x) - S_n| = |cos(x) - p_{2n+1}| = |R_{2n+1}|. For a general function f(x), the error term (the Lagrange formula for the remainder) is f(x) - p_{2n+1}(x) = R_{2n+1} = f^{2n+2}(c) x^{2n+2}/(2n+2)! where c is a number strictly between 0 and x. But if f(x)=cos(x), |f^{2n+2}(c)|=|cos(c)| <= 1 so that |R_{2n+1}| <= |x|^{2n+2}/(2n+2)! = |p_{2n+2}| = |a_{n+1}|. Thus in this instance, the error has the magnitude of the next term neglected. The basic way to make the series sum quickly is to sum the series over a limited range of x. Knowing that 0< deg <= 360 we know that 0< x <= 2 pi. So we evaluate the sum when |x| <= pi / 2, but if pi / 2 < x < 3 pi / 2 we do cos x = - cos ( pi - x ) instead because | pi - x | < pi / 2. Similarly if 3 pi / 2 <= x ( <= 5 pi / 2 ) we do cos x = cos ( 2 pi - x ) instead because | 2 pi - x | <= pi / 2. There is still the matter of how many terms are needed to get the desired accuracy. In the subroutine mycos, since the next term is an estimate of the current error, we loop through the sum as long as the next term |a| >= .0000005. When the desired accuracy for the sum is reached, then mycos exits and returns the sum. We have reduced the computation of the sum for the range |x| < pi/2 < 1.58 only. For this range |R_{11}| <= 5.053e-07 but |R_{9}| <= 2.072e-05. Hence, for every possible x in this range we need at most five therms (up to x^{10}). We can loop over five terms, or what is equivalent, evaluate the polynomial p_{10}(x). Now, the polynomial can be evaluated even faster if we store the constants c_i first, and regroup the computation as follows: Let s = x*x. p_5 = 1 + c2 * x^2 + c4 * x^4 + c6 * x^6 (16 mult and additions) = s*(s*(s*c6 + c4) + c2) + 1 (7 additioms and multiplications: MUCH BETTER!) This is coded in fastcos. The actual cosine routine in a calculator or canned in C is similar: the input x is reduced to a small range, and so a fixed number of terms is all that is needed to reach the fixed accuracy. A fixed degree Maclaurin polynomial is evaluated and returned. Actually, there is an improvement of Taylor polynomials, called Chebychev polynomials which does an even better job! The remaining thing to comment on is that occasionally the functions give different values. This cannot be avoided, because, irregardless of accuracy, sometimes the total computed from mycos is rounded up or down, even though the actual value would have been rounded the other way. Observe that fastcos agrees with machine cos much better than mycos, since fastcos often has better accuracy. */ # include # include # include # define C2 -0.5 # define C4 4.166666666666666e-02 # define C6 -1.388888888888889e-03 # define C8 2.480158730158730e-05 # define C10 -2.755731922398589e-07 # define PIOT 1.570796326794897 # define TP 6.283185307179586 double mycos( double x ); double fastcos( double x ); int main(void) { int deg, degeq, n = 1; double x, factor; printf ("\n Degree Mycos Fastcos Machine cos\n" ); for ( deg = 1; deg <= 360 ; deg = deg + n ) { x = deg * M_PI / 180.; printf ( "%7d %10f %10f %10f\n", deg, mycos(x), fastcos(x), cos(x) ); } return EXIT_SUCCESS; } /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ double mycos( double x ) { double sum = 1.0, y, a, s; int i = 3; if( fabs ( x ) <= PIOT ) { y = - x * x; a = 0.5 * y; do { sum = sum + a; a = y * a / ( i * ( i + 1 ) ); i = i + 2; } while ( fabs ( a ) >= .0000005 ); return sum; } else { s = x - M_PI; if( fabs ( s ) < PIOT ) return -mycos ( s ); else return mycos ( TP - x ); } } /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ double fastcos( double x ) { double s; if ( fabs ( x ) <= PIOT ) { s = x*x; return s*(s*(s*(s*(s*C10 + C8) + C6) + C4) + C2) + 1.0; } else { s = x - M_PI; if ( fabs ( s ) < PIOT ) return -fastcos ( s ); else return fastcos ( TP - x ); } } /***************************************************************************** Treibergs Feb. 24, 2006 Solution to Problem 3.2, to write a fuunction subprogram that return a p-periodic function of x. The main program will ask for the endpoints of an interval and print out 31 equally spaced values from a to b. ******************************************************************************/ # include # include double /* Function prototype. */ perfun ( double x, double p ); int main ( void ) { double a, b, c, d, p, x; int i; printf ( "Periodic Function Via Recursive Function Calls\n" ); printf ( " Enter period : " ); scanf ( "%lf", &p ); printf ( " Enter two endpoint values : " ); scanf ( "%lf %lf", &a, &b ); if ( b < a) { /* Swap a and b if b < a. */ c = a; a = b; b = c; } d = ( b - a ) / 30.0; printf ( "\n Listing of function values from %f to %f\n", a, b ); printf ( "\t\t x\t\t\t f(x)\n" ); for ( i = 0; i <= 30 ; i++ ) { x = a + i * d; printf ( "%25.15f %25.15f\n", x, perfun ( x, p ) ); } return EXIT_SUCCESS; } /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ double perfun ( double x, double p ) { /* perfun grows linearly from 0 to p/4, then droops. */ double result; if ( x < 0.0 ) result = perfun ( x + p, p ); else if ( x >= p ) result = perfun ( x - p, p ); else { if( x <= 0.25*p) result = 16.0 * x / p; else result = 1.0 / ( 0.25 + x / p ) - 1.0; } return result; }