Assignment 2
A Small Numerical Library
Prof. Max Dunne
CSE 13S – Fall 2020
1 Introduction
As we know, computers are simple machines that carry out a sequence of very simple steps, albeit very
quickly. Unless you have a special-purpose processor, a computer can only compute addition, subtrac-
tion, multiplication, and division. If you think about it, you will see that the functions that might interest
you when dealing with real or complex numbers can be built up from those four operations. We use
many of these functions in nearly every program that we write, so we ought to understand how they are
created.
If you recall from your calculus class, with some conditions a function f(x) can be represented by its
Taylor series expansion near some point f(a):
f(x)= f(a)+
∞∑
k=1
(x−a)k
k!
f(k)(a).
If you have forgotten (or never taken) calculus, do not despair. Go to a laboratory section for review: the
concepts required for this assignment are just derivatives.
Since we cannot compute an infinite series, we must be content to calculate a finite number of terms.
In general, the more terms that we compute, the more accurate our approximation. For example, if we
expand to 10 terms we get:
f(x)= f(a)+(x−a)f ′(a)+
1
2
(x−a)2f ′′(a)+
1
6
f(3)(a)(x−a)3
+
1
24
f(4)(a)(x−a)4+
1
120
f(5)(a)(x−a)5+
1
720
f(6)(a)(x−a)6
+
f(7)(a)(x−a)7
5040
+
f(8)(a)(x−a)8
40320
+
f(9)(a)(x−a)9
362880
+O((x−a)10) .
Taylor series, named after Brook Taylor, requires that we pick a point a where we will center the
approximation. In the case a= 0, then it is called a Maclaurin series). Often we choose 0, but the closer
to the value of x the better we will approximate the function. For example, let’s consider ex centered
around 0:
ex= 1+x+
x2
2
+
x3
6
+
x4
24
+
x5
120
+
x6
720
+
x7
5040
+
x8
40320
+
x9
362880
+O(x10) .
© 2020 Darrell Long 1
This is one of the simplest series when centered at 0, since e0= 1. Consider the general case:
ex= ea+ea(x−a)+
1
2
ea(x−a)2+
1
6
ea(x−a)3+
1
24
ea(x−a)4+
1
120
ea(x−a)5+
1
720
ea(x−a)6
+
ea(x−a)7
5040
+
ea(x−a)8
40320
+
ea(x−a)9
362880
+
ea(x−a)10
3628800
+O((x−a)11) .
Since ddxe
x = ex the exponential function does not drop out as it does for a = 0, leaving us with our
original problem. If we knew ea for a ≈ x then we could use a small number of terms. However, we do
not know it and so we must use a= 0.
What is theO((x−a)11) term? That is the error term that is “on the order of” the value in parentheses.
This is different from the big-O that we will discuss with regard to algorithm analysis.
2 Your Task
For this assignment, you will be creating a small numerical library. Our goal is for you to have some idea
of what must be done to implement functions that you use all of the time.
You will be writing and implementing sin, cos, tan, and exp using the Taylor series approximations
for exp and Padé approximants for sin, cos, and tan. You will then run these functions and compare
them to the standard library math.h implementations and output the results into a table similar to what
is shown in Figures 1 and 2.
1 x Sin Library Difference
2 - --- ------- ----------
3 -6.2832 0.00000000 0.00000000 -0.0000000000
4 -6.0868 0.19509032 0.19509032 -0.0000000000
Figure 1: Example of program output for sine.
1 x Cos Library Difference
2 - --- ------- ----------
3 -6.2832 1.00000000 1.00000000 0.0000000000
4 -6.0868 0.98078528 0.98078528 0.0000000000
Figure 2: Example of program output for cosine.
From left to right, the columns represent the input number, your program’s cosine value from the
input number, the actual math library’s value from the input number and lastly, the difference between
your value and the library’s value. To ensure easier grading use the format specifier below for printing
the values in your tables.
1 printf("% 6.4lf\t% 10.8lf\t% 10.8lf\t% 12.10lf\n", ...);
You will test sine and cosine from −2pi to 2pi with steps of pi/16. Tangent will be tested from −(pi/3)
to pi/3with steps of pi/16.1 For the exponential function, ex will be from 0 to 9with steps of 0.1.
1The slightly odd range is to make the table look nicer when printed.
© 2020 Darrell Long 2
Each implementation will be a separate function. You must name the functions Exp, Sin, Cos, and
Tan. Since the math library uses exp, sin, cos, and tan, you will not be able to use the same names. To
use math.h in your program you must be sure to link the math library at compile time using the -lm flag.
The following is an example function that implements Newton’s method of computing square roots
that doesn’t conflict with sqrt() found in math.h. Note that the function is named Sqrt(double x).
1 #define EPSILON 0.00001
2 static inline double abs(double x) { return x < 0 ? -x : x; }
3
4 double Sqrt(double x) {
5 double y = 1.0;
6 double old = 0.0;
7 while (abs(y - old) > EPSILON) {
8 old = y;
9 y = 0.5 * (y + x / y);
10 }
11 return y;
12 }
Computing
p
x using Newton’s method.
2.1 Sine and Cosine
The domain of sin and cos is [−2pi,2pi], and so centering them around 0makes sense. Since the domain
is restricted, you should reduce any parameter to [−2pi,2pi] (making your approximation better). The
Taylor series for sin(x) centered about 0 is:
sin(x)= x−
x3
6
+
x5
120
−
x7
5040
+
x9
362880
−
x11
39916800
+
x13
6227020800
+O(x14)
and the corresponding series for cos(x) is:
cos(x)= 1−
x2
2
+
x4
24
−
x6
720
+
x8
40320
−
x10
3628800
+
x12
479001600
+O(x14) .
We can use what is called a Padé Approximant. It’s beyond the scope of this course to go into com-
puting them, but essentially it is the ratio of two polynomials that conform to certain properties. It is
often easier to compute and more accurate than a truncated series. The Padé approximant for a 14 term
series for sin(x) centered around 0 is:
sin(x)≈ −479249x
7+52785432x5−1640635920x3+11511339840x
7(2623x6+453960x4+39702960x2+1644477120)
.
It is a lot easier to square a number than to raise it to a power, so we can simplify it by putting the
formula into Horner normal form, by factoring out x as much as possible:
sin(x)≈ x((x
2 (52785432−479249x2)−1640635920)x2+11511339840)
((18361x2+3177720)x2+277920720)x2+11511339840
.
Why Horner normal form? It has the fewest multiplications.
© 2020 Darrell Long 3
If you want to be clever you can compute x2 once, store it in a variable, using that result again and
save a few instructions. Does this matter? A good compiler will recognize the common subexpression
and do it for you behind the scenes, but numerical code tends to be heavily used so every little bit helps.
Consider the corresponding approximant for cos(x) centered around 0 written in Horner normal
form:
cos(x)≈ (x
2 (1075032−14615x2)−18471600)x2+39251520
((127x2+16632)x2+1154160)x2+39251520
.
Restricting the domain is important. As we move away from the center of the series, our approxima-
tion gets worse. Consider Figure 3, and you can see that when we get much beyond [2pi,2pi] the graphs
diverge not just a little, but wildly.
-10 -5 5 10
-4
-2
2
4
sinx versusPadéapproximant
Figure 3: Comparing sin(x) with an order 10 Padé approximant.
2.2 Tangent
You will recall that tan(x) = sin(x)/cos(x) but that it is undefined when cos(x) = 0, that is when x is a
multiple of pi/2. You could just do the division:
1 double tan(double x) { return sin(x) / cos(x); }
While doing the division is correct in the real numbers, remember that we are working with approxi-
mations and so you will be more accurate with a formula for directly computing your approximation for
tan(x). A Padé approximant for tan(x) is:
tan(x)≈ x(x
8−990x6+135135x4−4729725x2+34459425)
45(x8−308x6+21021x4−360360x2+765765)
.
You will want to rewrite it in Horner normal form before using it.
© 2020 Darrell Long 4
2.3 ex
The other important function is ex called the exponential function. While you might hope for a Padé
approximant, you are bound for disappointment. The domain of ex is unrestricted so we will not be able
to use any fixed formula. As we get further away from where we center our series, the more terms that
will need to get good results. Fortunately, the series for ex converges very quickly.
Let’s consider how well this works. In Figure 4, we will use our expansion to 10 terms and plot for
e0, . . . ,e10. After about 7, the approximation starts to diverge significantly. What this tells us is that 10
terms is not enough.
2 4 6 8 10
x
2000
4000
6000
8000
ex
TaylorApproximation
Figure 4: Comparing ex with its Taylor approximation centered at zero.
If we are naïve about computing the terms of the series we can quickly get into trouble — the values
of k! get large very quickly. We can do better if we observe that:
xn
n!
=
xn−1
(n−1)!
× x
n
.
At first, that looks like a recursive definition (and in fact, you could write it that way, but it would be
wasteful). As we progress through the computation, assume that we know the previous result. We then
just have to compute the next term and multiply it by the previous term. At each step we just need to
compute x/n, starting with n= 0! = 1 and multiply it by the previous value and add it into the total. It
turns into a simple for or while loop.
We can use an (epsilon) to halt the computation since |xn|Figure 5: briefly, xn dominates but is quickly overwhelmed by n! and so the ratio rapidly approaches
zero. For this assignment, = 10−9 will be sufficient.
© 2020 Darrell Long 5
2 4 6 8 10
5
10
15
20
25
Figure 5: Comparing xn/n! for x= 2,3,4,5.
Pre-lab Part 1
1. Write pseudo-code for approximating ex with either a for or while loop.
2. Write pseudo-code for printing the output for ex.
2.4 Command-line Arguments
Your program will use command-line arguments to select which one of your implemented functions to
run. You will use getopt() to parse command-line arguments, which is included under getopt.h. In
most C programs, you will see two parameters in the main() function: int argc and char **argv.
The parameter argc refers to the argument counter, or how many arguments were supplied on the
command-line. Arguments are delimited by white space, which includes spaces and tabs. The parameter
argv refers to the argument values, and is interpreted as an array of strings, where argv[0] corresponds
to the first argument on the command-line: the executable itself.
Try this code, and make sure that you understand it:
1 #include
2
3 int main(int argc , char **argv) {
4 for (int i = 0; i < argc; i += 1) {
5 printf("argv[%d] = %s\n", i, argv[i]);
6 }
7 return 0;
8 }
© 2020 Darrell Long 6
Command-line options must be defined in order for getopt() to parse them. These options are de-
fined in a string, where each character in the string corresponds to an option character that can be speci-
fied on the on the command-line. Upon running the executable, getopt() scans through the command-
line arguments, checking for option characters.
As an example, the following program supports two command-line options, ‘p’ and ‘i’. Note that the
option character ‘i’ in the defined option string OPTIONShas a colon following it. The colon signifies that,
when the ‘i’ option is enabled on the command-line using ‘-i’, getopt() is looking for an argument to
be supplied following it. An error is thrown by getopt() if an argument for a flag requiring one is not
supplied.
1 #include
2 #include
3 #include
4
5 #define OPTIONS "pi:"
6
7 int main(int argc , char **argv) {
8 int c = 0;
9 bool print = false;
10 char *infile = NULL;
11
12 while ((c = getopt(argc , argv , OPTIONS)) != -1) {
13 switch (c) {
14 case ’p’: // Print option.
15 print = true;
16 break;
17 case ’i’: // Input file option.
18 infile = optarg; // A pointer to the next element in argv.
19 break;
20 }
21 }
22
23 if (argc == 1) {
24 puts("Error: no arguments supplied!");
25 return -1;
26 }
27
28 if (! infile) {
29 puts("Error: input file required!");
30 return -1;
31 }
32
33 if (print) {
34 puts(infile);
35 }
36
© 2020 Darrell Long 7
37 return 0;
38 }
What this program does is check for the ‘p’ option to print and the ‘i’ option to specify some input file
name. The getopt() defined variable optarg points at the following argv-element, which in this case
should be the input file name. The condition argc == 1 checks if no other command-line arguments
were supplied other than the executable itself and errors out if so. The program also errors out in the
event that an input file name was not supplied. The supplied input file name is only printed if the print
option was specified on the command-line.
Example usage of program to specify input file and print its name (assume the executable name is
filename):
1 ./ filename -p -i
Your program must support the following command-line options:
1. -s : runs sin tests
2. -c : runs cos tests
3. -t : runs tan tests
4. -e : runs exp tests
5. -a : runs all tests
These options are mutually exclusive; only one may be chosen at a time. Is a bool the best choice, or
would an enum be better?
Pre-lab Part 2
1. What does getopt() return? Hint: check the man page.
2. Is a bool or an enum the best choice? Explain why.
3. Provide pseudo-code for your main function. Assume you have helper functions available
to you.
3 Deliverables
You will need to turn in:
1. math.c: This is your main program that determines whether exp, sin, cos, or tan will be run as
well as contain your implementations. getopt() must be used to handle command-line argu-
ments dictating which tests to run. The getopt() options you must support are ‘-s’ to run sin
tests, ‘-c’ to run cos tests, ‘-t’ to run tan tests, ‘e’ to run exp tests, and ‘-a’ to run all tests. Your
output for each function must match the form shown in Figures 1 and 2. Your compiled program
must be called math.
© 2020 Darrell Long 8
2. Makefile: This is a file that will allow the grader to type make to compile your program. Typing
make must build your program and running ./math alone as well as with option flags must run
your program.
• CFLAGS=-Wall -Wextra -Werror -Wpedantic must be included.
• CC=clang must be specified.
• make clean must remove all files that are compiler generated.
• make infermust run infer on your program. Complaints generated by infer should either
be fixed or explained in your README.
• make should build your program, as should make all.
3. README.md: This must be inmarkdown. This must describe how to use your program and Makefile.
This is also where you put any explanations for complaints generated by infer.
4. DESIGN.pdf: This must be a PDF. The design document should contain answers to the pre-lab
in the beginning and describe your design for your program with enough detail that a sufficiently
knowledgeable programmer would be able to replicate your implementation. This does not mean
copying your entire program in verbatim. You should instead describe how your program works
with supporting pseudo-code. For this assignment pay extra attention to how you describe your
representation of a Taylor series approximation in C. You must push the DESIGN.pdf before you
push any code.
5. WRITEUP.pdf: This must be a PDF. Your WRITEUP should be a discussion of the results for your
experiments. This means analyzing the differences in the output of your implementations versus
those in the math.h library. Include possible reasons for the differences between your implemen-
tation and the math.h version. Graphs can be especially useful in showing the differences and
backing up your arguments.
4 Submission
To submit your assignment, refer back to asgn0 for the steps on how to submit your assignment through
git. Remember: add, commit, and push!
Your assignment is turned in only after you have pushed. If you forget to push, you have not turned in
your assignment and you will get a zero. “I forgot to push” is not a valid excuse. It is highly recommended
to commit and push your changes often.
5 Supplemental Readings
• The C Programming Language by Kernighan & Ritchie
– Chapter 3 §3.4-3.7
– Chapter 4 §4.1 & 4.2
– Chapter 7 §7.2
– Appendix B §B4
© 2020 Darrell Long 9

欢迎咨询51作业君