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