(Last Mod: 27 November 2010 21:38:37 )
Understand how to use functions from the standard math and other libraries.
Be able to create new functions that perform tasks not covered by functions in the math library.
Understand basic iterative techniques for solving non-algebraic problems.
Just as the standard header file stdio.h gives us access to a variety of useful function for performing input and output tasks, other standard header files give us access to a wide range of other functions that perform a variety of useful tasks. The C language standard presently requires a minimum set of twenty-four standard header files and most C implementations contain others that provide vendor-specific extensions.
Some of the more commonly used standard header files are listed below:
| Header File | Purpose |
| stdio.h | Input and output functions. |
| stdlib.h | General utility functions such as memory management functions. |
| stddef.h | A variety of symbolic constants and type definitions. |
| math.h | Mathematical functions. |
| limits.h | Integer limitations. |
| float.h | floating point limitations. |
| ctype.h | character functions. |
| string.h | string manipulation functions. |
| time.h | time functions. |
In this lesson we will use some of the math.h functions as examples of how to access and work with the functions in the standard libraries. We will also create a few functions of our own that perform tasks not addressed by the standard functions.
Let's suppose that we wanted to generate a table of trigonometric values similar to the following:
| DEG | SINE | COSINE | TANGENT |
=======================================
| 0.0 | 0.00000 | 1.00000 | 0.00000 |
| 0.5 | 0.00873 | 0.99996 | 0.00873 |
| 1.0 | 0.01745 | 0.99985 | 0.01746 |
| 1.5 | 0.02618 | 0.99966 | 0.02619 |
....
| 89.0 | 0.99985 | 0.01745 | 57.28996 |
| 89.5 | 0.99996 | 0.00873 | >= 100.0 |
| 90.0 | 1.00000 | 0.00000 | >= 100.0 |
The pseudocode for generating this table might be similar to the following:
TASK: Generate the table header
TASK: Generate the table itself
SET: Angle = 0.0
WHILE: Angle <= 90.0
TASK: Output Angle with column delimiters
TASK: Output sin(Angle) with column delimiters
TASK: Output cos(Angle) with column delimiters
TASK: Output tan(Angle) with column delimiters
SET: Angle = Angle + 0.5
While we have already seen how to write a function that can display a floating point value to a specific number of decimal places, this is a good time to introduce the standard library function from stdio.h that performs a similar task. That function, printf(), is very versatile but requires a bit of practice to get comfortable with. For our specific situation, it is enough to note the following:
To print a value of type double to 3 decimal digits beyond the decimal point using a total of 6 spaces in the output, we can use the following statement:
printf("%6.3f", x);
Where x is either a variable of type double or an expression that evaluates to a numerical value. Without the use of additional flag characters in the format string (the first argument to the printf() function), the value will be right-justified within the 6 spaces used and any unused space on the left will be padded with blank spaces. It's important to note that if the value requires more than 6 spaces to display, it will use more than 6 spaces. This "width" specifier only sets a minimum width on the output field. The number of digits to the right of the decimal point is referred to as the "precision" of the format specifier.
Looking at the "DEG" column in the table, we see that the largest value that gets output is 90.0, which requires a field width of four spaces (don't forget to count the decimal point) with one decimal place after the value. Hence the output required by Task #2.1 can be performed with the following statement:
printf("%4.1f", Angle);
While the above statement will print out the value, it will not print out the additional space nor the delimiters that are contained in the table. This can be accomplished by adding two additional printf() statements:
printf("| ");
printf("%4.1f", Angle);
printf(" |");
In general, the format contains one of two types of information - characters that are to be output directly and format instructions that describe how the value passed as the subsequent argument are to be output. The function can tell the difference by looking for the '%' character which marks the beginning of a "format specifier". All of the characters from that point until the end of the specifier tell the printf() function how to display the value. The end of the format specifier is marked by one of a short list of characters. In the case above, the 'f' character marks the end of the specifier and tell the printf() function that the value that is passed as the second argument is of type double.
Since the printf() function has the ability to parse the format string and identify format specifiers that are embedded within it, the above three functions can be combined into a single statement:
printf("| %4.1f |", Angle);
Let's say that we have three variables of type double called sine, cosine, and tangent. After setting these variables to the appropriate values for the current value of Angle, we could then perform the next three tasks in a manner very similar to the above using the following statements:
printf("| %4.1f |", Angle);
printf(" %7.5f |", sine);
printf(" %7.5f |", cosine);
printf(" %8.5f |\n", tangent);
Notice that we increased the precision of the final three statements to match out desired output results and that we increased the total width correspondingly including an additional width of one for the tangent functions to accommodate tangent values that have two digits to the left of the decimal point. Also, notice that we only printed the ending delimiter in the last three statements since the prior column's ending delimiter is also the next column's starting delimiter. Hence we want to print one or the other but not both. Second, notice the '\n' character at the end of the last statement. Since these format strings are string literals, the same escape sequences are used here as in all string literals.
At this point we are very close to being completely finished. Using the printf() statement to generate our table headings, our nearly complete program is:
/* TASK: Generate the table header */
printf("| DEG | SINE | COSINE | TANGENT |\n");
printf("=======================================\n");
/* TASK: Generate the table itself */
Angle = 0.0;
while (Angle <= 90.0)
{
printf("| %4.1f |", Angle);
printf(" %7.5f |", sine);
printf(" %7.5f |", cosine);
printf(" %8.5f |\n", tangent);
Angle += 0.5;
}
The only two items left are getting the actual values for the trig functions and dealing with overly large values of the tangent function. The first of these can be accomplished by using the trig functions provided with the compiler and accessible through function prototypes contained in the math.h header file. These functions are:
double sin(double x);
double cos(double x);
double tan(double x);
Each function takes as its single argument an angle measured in radians and returns the appropriate trig value for that angle. Since our table is in degrees, we need to convert from degrees to radians before we call the function. A rather brute force way to accomplish this would be as follows:
#include <math.h> /* sin(), cos(), tan() */
...
while (Angle <= 90.0)
{
radians = (3.14159/180.0)*Angle;
sine = sin(radians);
cosine = cos(radians);
tangent = tan(radians);
printf("| %4.1f |", Angle);
printf(" %7.5f |", sine);
printf(" %7.5f |", cosine);
printf(" %8.5f |\n", tangent);
Angle += 0.5;
}
We can clean this up by recalling that the printf() function's subsequent arguments can be expressions as well as single variables. We can also use a couple of #define macros to our advantage here:
#include <math.h> /* sin(), cos(), tan() */
#define PI (3.141592653589793238)
#define Deg2Rad(theta) ( (theta)*PI/180.0 )
....
while (Angle <= 90.0)
{
printf("| %4.1f |", Angle);
printf(" %7.5f |", sin(Deg2Rad(Angle)));
printf(" %7.5f |", cos(Deg2Rad(Angle)));
printf(" %8.5f |\n", tan(Deg2Rad(Angle)));
Angle += 0.5;
}
The only remaining task is to deal with large values of tangent. To tackle this, we'll further expand the pseudocode for that task:
TASK: Output tan(Angle) with column delimiters
IF: tan(Angle) >= 100
OUTPUT: " >= 100.0 "
ELSE:
OUTPUT: tan(Angle)
With this, the modification that we need is very clear and our final, complete program can now be written:
#include <stdio.h> /* printf() */
#include <math.h> /* sin(), cos(), tan() */
#define PI (3.141592653589793238)
#define Deg2Rad(theta) ( (theta)*PI/180.0 )
int main(void)
{
double Angle;
/* TASK: Generate the table header */
printf("| DEG | SINE | COSINE | TANGENT |\n");
printf("=======================================\n");
/* TASK: Generate the table itself */
Angle = 0.0;
while (Angle <= 90.0)
{
printf("| %4.1f |", Angle);
printf(" %7.5f |", sin(Deg2Rad(Angle)));
printf(" %7.5f |", cos(Deg2Rad(Angle)));
if ( 100.0 <= tan(Deg2Rad(Angle)) )
printf(" >= 100.0 |\n");
else
printf(" %8.5f |\n", tan(Deg2Rad(Angle)));
Angle += 0.5;
}
return 0;
}
NOTE: The value of PI in the #define statement is pi to 19 decimal places. This matches the level of accuracy for a variable of long double on the Borland TurboC/C++ v4.5 compiler.
We might well discover that, for our needs, it would be very convenient to have a set of trig functions that take arguments in degrees instead of radians. No such library functions exists. But we are free to define our own functions that perform the necessary conversions as part of the function so that we don't have to deal with them at the level that we use the functions. Since we can't use the same function names, we need to come up with something different. While we can call them whatever we want, we should consider something reasonably descriptive such as sin_d(). Using this nomenclature, our functions might be:
double sin_d(double degrees)
{
return sin(Deg2Rad(degrees));
}
double cos_d(double degrees)
{
return cos(Deg2Rad(degrees));
}
double tan_d(double degrees)
{
return tan(Deg2Rad(degrees));
}
Of course, the above code requires that the math.h header file be included as well as out two macros.
We probably also want a set of inverse functions that yield results in degrees. As you should already suspect, a set of inverse trig functions (which happen to be asin(), acos(), and atan()) is made available by math.h but the values returned are angles in radians. Just as we used a macro to convert a value from degrees to radians, we define another to convert a value from radians to degrees. With this, we can write our degree-based inverse trig functions as follows:
#define Rad2Deg(radians) ((radians)*(180.8/PI))
double asin_d(double x)
{
return Rad2Deg(asin(x));
}
double acos_d(double x)
{
return Rad2Deg(acos(x));
}
double atan_d(double x)
{
return Rad2Deg(atan(x));
}
Occasionally is it useful to have a single function that has several different behaviors that are chosen based on one of the parameters that is passed to it. The function printf() would clearly fall into this category. To illustrate how we might implement such a function, consider the following problem statement:
Write a function called Sin() that performs both forward (sin()) and inverse (arcsin()) computations using either radians or degrees. The first argument is the basic value that is being operated on and the second value (the mode) is an integer that controls the function's behavior according to the following table:
| mode | behavior |
| 0 | sin(x) in radians |
| 1 | sin(x) in degrees |
| 2 | arcsin(x) in radians |
| 3 | arcsin(x) in degrees |
The behavior of the function for any value of mode not in the table is unspecified.
While there are many ways to implement such a function, the switch() statement is commonly a natural candidate for such tasks. Leveraging our prior macros and functions, this function is very straight-forward:
double Sin(double x, int mode)
{
switch (mode)
{
case 0: return sin(x);
case 1: return sin(Deg2Rad(x));
case 2: return asin(x);
case 3: return Rad2Deg(asin(x));
}
return sin(x); /* reasonable "unspecified" behavior */
}
The standard math library provides functions for the purpose of computing logs and anti-logs in both base 10 and base e. The prototype for the functions are as follows:
double log(double x); /* natural log of x */
double exp(double x); /* e^x */
double log10(double x); /* common log of x */
double pow10(double x); /* 10^x */
Let's say that we want to create a pair of functions that compute the log and anti-log in an arbitrary base b. Clearly, our functions will need to take a second argument, namely the value of the base. To use a nomenclature that is consistent with the standard libraries, we will declare our functions as follows:
double log_b(double x, double b); /* base b log of x */
double pow_b(double x, double b); /* b^x */
By recalling the fundamental definition of a logarithm, we can quickly develop the code for these functions:
The value of logb(x) is given by the relationship:
y = logb(x) iff x = by
Taking the log (in any base) of both sides yields:
log(x) = y*log(b)
Solving for y yields:
y = logb(x) = log(x)/log(b)
Hence our first function is extremely simple:
double log_b(double x, double b)
{
return log(x)/log(b);
}
The second function is almost as simple. Once again going back to the fundamental relationship we are working with, we want the value y such that:
y = bx
Taking the log of both sides, we have:
log(y) = x * log(b)
Taking the anti-log of both sides, we recover the value y that we are interested in:
y = exp( x * log(b) )
double pow_b(double x, double b)
{
return exp( x * log(b) );
}
An examination of the functions available in the math library reveals the function pow() which takes two arguments and returns the value obtained by raising the first argument to the power given by the second argument.
double pow(double x, double y); /* x^y */
With this knowledge, would could rewrite our pow_b() function to simply be:
double pow_b(double x, double b)
{
return pow(b, x);
}
Up to this point, we have dealt with functions that we could express directly in terms of functions provided by the standard math library. This is not always the case. For instance, let's say that we want a pair of functions that perform x*log(x) and its inverse. The forward function is as straight-forward as any of the other functions we have thus worked with:
double xlogx(double x)
{
return x * log(x);
}
But what about the inverse function? Namely, given a value y, find the value of x that satisfies the relationship:
y = x*log(x)
There is no algebraic solution to this problem. Indeed, most of the functions we have been using in this lesson, such as the trigonometric, exponential, and logarithm functions, are non-algebraic. There are two basic approaches to solving such a problem - develop an infinite series expansion for the expression and evaluate the number of terms necessary to achieve the desired level of accuracy. This is essentially the method used by the standard library functions.
A more generic way is to iteratively converge to a solution, assuming one exists, through some sort of systematic search. In general, there is no guarantee that such a search will actually find a solution even if one exists, but there are techniques that enjoy a reasonably high probability of success that they are worth understanding.
If the function is monotonic (always increasing or decreasing for increasing values of x), then we can use a very simply search procedure to find the answer.
Consider the function x*log(x). This function is defined for all values of x that are strictly positive. The reason for this is that the log(x) function itself is undefined for values less than or equal to zero. For values of x greater than 1, this is clearly a monotonic function since both x and log(x) increase monotonically. The situation is not so clear cut for values of x between zero and one, although it is clear that all such values are negative since x is positive but log(x) is negative for values less than one..
Limiting the discussion to positive values of y, meaning that the corresponding value of x will be greater than or equal to one, how can we find the value of x such that y = x*log(x)?
Instead of exploiting the behavior of this specific function, we will limit ourselves to techniques that will work, albeit perhaps not efficiently, on any monotonic function.
A good place to start with nearly any problem is to consider how we would do it by hand if forced to. In this case, we would probably start with an initial value of x and see if it works. If it didn't, we would consider it's relationship to the desired value and pick a new guess for x accordingly. So let's do that.
Find the value of x such that x*log(x) is equal to 100.
If we start with a trial value of x = 1, we get a result of 0. If we try a trial value of 2 we get a result of 1.386. We are moving in the right direction, but have a long way to go. At this point, let's start doubling the value of x until we get a result greater than 100:
| x | x*log(x) |
| 1 | 0 |
| 2 | 1.386294 |
| 4 | 5.545177 |
| 8 | 16.635532 |
| 16 | 44.361420 |
| 32 | 110.903549 |
We now know that the value is between 16 and 32. So let's split the difference and evaluate the function at 24. If the value is still greater than 100, we know that the value is between 24 and 32. Otherwise, we know the value is between 16 and 24. We will split the difference again going in the correct direction. We will repeat this process until we get as close to the correct value of x as we desire.
| x | x*log(x) |
| 1 | 0 |
| 2 | 1.386294 |
| 4 | 5.545177 |
| 8 | 16.635532 |
| 16 | 44.361420 |
| 32 | 110.903549 |
| 24 | 76.273292 |
| 28 | 93.301726 |
| 30 | 102.035921 |
| 29 | 97.651579 |
| 29.5 | 99.839513 |
| 29.75 | 100.936667 |
| 29.625 | 100.387826 |
| 16 iterations | |
| 29.536599 | 100.000000 |
We were able to find a value of x that yielded the target value of y to six decimal places (eight significant figures) with a total of 29 iterations. The first five find a region where our goal was within a a known range. We then subdivided that range in half and examined the midpoint of the appropriate half each time. Performing 16 iterations, even by hand (with a calculator), is not a significant burden and, in most circumstances, a computer program could easily crank through a thousand or more iterations if necessary.
How do we know when we are finished? There are a few reasonable ways to answer this question. If what we are looking for is a certain number of significant figures in the answer, then we could keep going until the ratio of the answer to the amount that we are changing it by indicates that we have achieved the desired level. For instance, if the answer is 10 and our present increment amount is 0.01, then the ratio is 1000 and indicates that we have achieved roughly three significant figures of accuracy. Notice that this is the base ten log of the ratio. Therefore, we could stop when the base ten log of the ratio exceeds the number of sig figs we are striving for.
So, assuming that this level of performance is acceptable, let's turn out hand-example into an algorithm:
TASK: Find the appropriate general range
SET: x = 1
WHILE: y > x*log(x)
SET: x = 2x
TASK: Zero in on the actual value
inc = 0.5x
WHILE: log10(x / inc) < SIGFIGS
IF: y > x*log(x)
SET: x = x + inc
ELSE
SET: x = x - inc
SET: inc = 0.5 * inc
The code for this function, directly from the above pseudocode, is then simply:
#define SIGFIGS (16)
....
double arc_xlogx(double y)
{
double x, inc;
/* TASK: Find the appropriate general range */
x = 1.0;
while (y > x*log(x))
x *= 2.0;
/* TASK: Zero in on the actual value */
inc = 0.5*x;
while ( log10(x/inc) < SIGFIGS )
{
if (y > x*log(x))
x += inc;
else
x -= inc;
inc *= 0.5;
}
return x;
}
What if we wanted to reduce the number of iterations that are required? Is it possible?
To answer this question we once again look to our hand example and ask how we might do it by hand. Once we have found a pair of bounding value for x - i.e., a value of x that we know is too small and a value of x that we know is too large, our next guess was simply the midpoint between those two values. While this gave us a simple recipe to follow, we are clearly throwing away useful information. For instance, consider the very first time we found a value of x such that x*log(x) was greater than our target value of 100.
| x | x*log(x) |
| 1 | 0 |
| 2 | 1.386294 |
| 4 | 5.545177 |
| 8 | 16.635532 |
| 16 | 44.361420 |
| 32 | 110.903549 |
Just looking at the values of x*log(x) for our two bounding values of x, we can surmise that the value of x we are looking for is closer to 32 than it is to 16. What if the function we were using was a simple straight line that passed through these two bounding points? Would it be possible to calculate the exact value of x that we want in that case? Clearly the answer is yes. While this won't give us an exact value when our functional relationship is not linear, it will generally give us a better estimate than just taking the midpoint. In fact, as the difference between our bounding values shrink, our linear estimation will become increasingly more accurate because we know that, for any continuous function, the function can be approximated by a straight line between any two sufficiently close limits - that is, after all, part of the fundamental underpinnings of both differential and integral calculus.
The method involved is known as "linear interpolation" and proceeds as follows:
Given two pairs of values <x1, y1> and <x2, y2>, find the value of x that corresponds to a value of y provided that y1 < y < y2 and assuming a linear relationship between x and y.
A linear relationship implies y = mx+b
Since we want the value of x that corresponds to a known value of y, we rewrite this as:
x = (y-b)/m
Now we just use our two known endpoints to compute the values of m and b and we are finished (more or less).
Recall that m is the slope of the line and is therefore equal to "rise over run":
m = (y2 - y1) / (x2 - x1)
We can now solve for the value of b using either of the two known endpoints. For instance:
b = y - mx = y2 - m*x2
Now recall that b is the y-intercept, which is the value of y when x is equal to zero:
Our algorithm is very similar to our previous one. The major difference is that before we simply used the midpoint and then used a shortcut in order to pick the upper or lower half of the range by checking the value at the midpoint against the target value and then either moving halfway into the lower range or halfway into the upper range as appropriate. Here the concept is nearly identical but we don't have as clear a shortcut.
Now we have a range, x1 to x2, and we compute a new estimate, x, somewhere in between. Our actual value is either between x1 and x (the lower portion) or between x and x2 (the upper portion). We make the decision the same way as before, if our target value is greater than x*log(x), then we know that our actual value lies in the upper half. Otherwise it lies in the lower half. If it lies in the upper half, then our new range has the same value of x2, but x1 must be adjusted to be equal to x. Similarly, if our value lies in the lower half, our new range has the same value of x1 but x2 must be set equal to x.
| x | x*log(x) |
| 1 | 0 |
| 2 | 1.386294 |
| 4 | 5.545177 |
| 8 | 16.635532 |
| 16 | 44.361420 |
| 32 | 110.903549 |
| 29.378251 | 99.305968 |
| 29.535144 | 99.993617 |
| 29.536586 | 99.999941 |
| 29.536599 | 99.999999 |
Whereas our prior algorithm, which was simpler, required 24 iterations after finding a pair of bounding values, this algorithm converges to the same result after only four iterations. So we have traded a bit of additional complexity at each iteration in exchange for a significant reduction in the total number of iterations that must be performed. In most cases, this is a favorable tradeoff.
Our new version, with no optimizations, is:
#define SIGFIGS (16)
....
double arc_xlogx(double y)
{
double x1, y1;
double x2, y2;
double xnew, ynew;
double m, b;
/* TASK: Find the appropriate general range */
x2 = 1.0;
y2 = 0.0;
do
{
x1 = x2;
y1 = y2;
x2 *= 2.0;
y2 = x2*log(x2);
} while (y > y2);
/* TASK: Zero in on the actual value */
do
{
m = (y2 - y1) / (x2 - x1);
b = y2 - m*x2;
xnew = (y-b)/m;
ynew = xnew*log(xnew);
if (y > ynew)
{
x1 = xnew;
y1 = ynew;
}
else
{
x2 = xnew;
y2 = ynew;
}
} while ( log10(xnew/(x2-x1)) < SIGFIGS );
return xnew;
}
Although there are a number of optimizations and refinements that could be made to the above function, there is two major flaws that could conceivably result in runtime errors. It is a good practice to examine any place where division is performed and consider the potential for division-by-zero errors. Likewise, all function calls should be examined for the potential of passing arguments that result in problems.
In the code above, we have a problem if (x2-x1) is too small or if (y2-y1) is too small. Also notice that if y is equal to b, that xnew will be zero which will cause problems with log(xnew). Since, in the first pass, we know that x2 will be twice the value of x1, and since we know that the function is monotonically increasing, we won't have a problem on the first pass.
The first problem spot is therefore the test for the while() loop. If (x2-x1) is zero, then we will have a divide-by-zero error. But if (x2-x1) is zero, it means that we have hit the the desired value of x exactly and therefore we are done. To guard against this error, we will utilize guaranteed short-circuiting and fail the loop if (x2-x1) is equal to zero or if the number of sig figs has been achieved. The other potential trouble spots are not a problem given the nature of the function we are using (a monotonically increasing function) and the fact that we have ensured that x2 is strictly greater than x1.
#define SIGFIGS (16)
....
double arc_xlogx(double y)
{
double x1, y1;
double x2, y2;
double xnew, ynew;
double m, b;
/* TASK: Find the appropriate general range */
x2 = 1.0;
y2 = 0.0;
do
{
x1 = x2;
y1 = y2;
x2 *= 2.0;
y2 = x2*log(x2);
} while (y > y2);
/* TASK: Zero in on the actual value */
do
{
m = (y2 - y1) / (x2 - x1);
b = y2 - m*x2;
xnew = (y-b)/m;
ynew = xnew*log(xnew);
if (y > ynew)
{
x1 = xnew;
y1 = ynew;
}
else
{
x2 = xnew;
y2 = ynew;
}
} while ( (x1 < x2) && ( log10(xnew/(x2-x1)) < SIGFIGS ) );
return xnew;
}