Computing the cosine function
-----------------------------
You know how to compute the cosine for several special
angles (0, 30, 45, 60, 90) and how to make use of symmetries
in the cosine function to find equivalent angles in the
other quadrants.
You can use these in conjunction with trig identities to
compute other angles. For example, you can use a half-angle
identity to compute cos(15 degrees).
How do you compute the cosine for an arbitrary angle, for
example 14.78 degrees?
The answer is to use the Maclaurin series expansion for cosine.
cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! - ......
In a large simulation (like a weather simulation or computing
the flow around an aircraft), a function like this might be called
several billion times.
It is critical that this be both fast (in terms of the number of
operations required) and accurate (no more than 1 bit of error in
the last place).
Version 1:
- Straightforward implementation of series.
- uses factorial function we wrote earlier in the semester to compute
the factorials.
- uses a while loop because we don't know how many terms
it will take for the series to converge.
- the idea is to keep adding terms until the sum is no longer changing.
---------------------------------------------
x = some angle in radians;
oldsum = 1;
newsum = 0;
n = 1;
flag = 0;
while(flag == 0)
numerator = (-1)^n * x^(2*n);
denominator = nfact(2*n);
newsum = oldsum + numerator/denominator;
if(newsum == oldsum)
flag = 1;
end
oldsum = newsum;
n = n + 1;
end
----------------------------------------------
Questions:
Q1) Suppose x = 10000? What happens if we try to compute
cos(10000)?
Answer: Result overflows. In exact arithmetic, the
series converges for any x because the (2*n)!
denominator eventually dominates the
x^(2*n) numerator. In finite arithmetic, the
denominator overflows before this can happen.
Solution: a) Use periodic nature of cosine function. For
any angle x, there is an equivalent angle in
the interval [0, 2*pi], so set
x = mod(x,2*pi)
b) Use symmetry properties of cosine function.
In Quadrant 2, cos(x) = -cos(pi - x)
In Quadrant 3, cos(x) = -cos(x - pi)
In Quadrant 4, cos(x) = cos(2*pi-x)
c) Using these, the angle can always be put into
the interval [0, pi/2], so overflow can't
occur.
Q2) Can the amount of work be reduced?
Answer: a) Yes. We are recomputing (2*n)! from scratch each pass through
the loop even though we have most of this available from
the previous calculation. For example, if we need to compute
the 6! in the denominator of the 3rd term, we are not using
the fact that we have 4! available from the previous
pass through the loop.
b) Something similar occurs with the x^(2*n) numerator. In computing
the 3rd term, we are computing x^6 from scratch even though
we have x^4 available from the previous pass through the loop.
c) The key is to efficiently generate the numerator and denominator
of the new term using information from the previous term.
Here, it can be seen that
new_numerator = -x^2 * old_numerator
new_denominator = n*(n+1) * old_denominator
d) In addition, since x only appears as x^2 in the formulas
above, we can replace this with xsq = x^2
new_numerator = -xsq * old_numerator
new_denominator = n*(n+1) * old_denominator
Version 2: More efficient
----------------------------------------------------------
x = some angle in radians;
x = (shift angle into equivalent angle in [0,pi/2]);
xsq = x*x;
oldsum = 1;
newsum = 0;
numerator = 1;
denominator = 1;
n = 1;
flag = 0;
while(flag = 0)
numerator = -xsq*numerator;
denominator = n*(n+1)*denominator;
newsum = oldsum + numerator/denominator;
if(newsum == oldsum)
flag = 1;
end
oldsum = newsum;
n = n + 2;
end
--------------------------------------------------------
More improvements:
1) Run for hundreds of angles and see how many terms it takes
the series to converge (ie, how many times does the loop
execute?). This will vary depending on the value of x.
Once the maximum number of terms is known for an average x,
the uncounted loop can be replaced with a counting loop.
This eliminates the IF test in the middle of the loop.
IF tests are a critical component to programming, but they
are kind of like speed bumps in that they slow the program
down. Anytime you can eliminate one, especially one that is
in a loop, you will gain speed (even at the cost of a few
operations).
2) Since you now have a counting loop, you can define a vector
of values that contains precomputed factorials.
3) If we also have a routine for computing the sine, we can use
the identity cos(x) = sin(pi/2 - x) and reduce the angle
even further to [0, pi/4]. The series will converge very
rapidly since x < 1 and the x^(2*n) terms will decay very
quickly.
Final version of cosine routine (see mycos.m file)
Other notes:
We have focused on speed but have not addressed accuracy.
1) gcc 4.5 code for cosine is on the website. This routine
uses nterms = 7 and also contains other accuracy tweaks.
This code uses a common trick to gain more digits. It takes
a number
z
and writes it as
z = x + y
where x contains the high digits and y contains the low digits,
for example, if
z = 0.1234567891234567
then
x = 0.12345678d0 y = 0.91234567d-8
x and y now have room for 8 additional digits in subsequent calculations.
Using an identity for
cos(x+y)
and some results from the remainder term for the Taylor series, the
last 2 or 3 digits have improved accuracy over our algorithm.
Other transcendental functions
1) Computing exp(x)
As with the cosine, exp(x) is computed using the Maclaurin series
exp(x) = 1 + x + x^2/2! + x^3/3! + x^4/4! + .....
This series also converges rapidly if x is small. In this case, we
don't have any periodic or symmetry properties to exploit, so overflow
is a big concern. This is addressed by using the relation
exp(x) = exp(x/2^n)^(2^n)
Since the computer uses binary, dividing and multiplyig by 2^n is
just a change in the decimal point (ie, no work is required). By
choosing n correctly for a given x, you can ensure that 0 <= x/2^n < 1.
The series will converge very rapidly and the same ideas as from the
cosine discussion above can be used.
2) Computing tan(x)
tan(x) = x + x^3/3 + 2*x^5/15 + ....
Unfortunately, the Maclaurin series for tangent converges too slowly to
be useful. To compute the tangent, one technique is to use a variation
on our cosine algorithm that also computes a sine for (almost) free via
the Pythagorean identity cos^2(x) + sin^2(x) = 1.
Then the identity
tan(x) = sin(x)/cos(x)
is used to get the tangent.
3) Computing ln(x)
ln(x-1) = x - x^2/2 + x^3/3 - x^4/4 + .......
The Taylor series for the natural log function has the same
problem as the tangent; it doesn't converge fast enough.
To obtain a rapid (and accurate) series for ln is not easy.
One of the most common techniques uses a complicated function
known as an elliptic integral.