C++ For Mathematicians (2006) [eng]
.pdf
36 |
C++ for Mathematicians |
37return d;
38}
Lines 23–25 handle the special case in which a is zero. After that, we declare a variable d to hold the answer.
Lines 31–35 do the bulk of the work of this program. We begin with a new C++ keyword: for. The general form of the for statement is this:
for( starting statement ; test condition ; advancing statement) { statements to be done on each iteration;
}
The starting statement is executed the first time the for statement is reached. In our case, the starting statement is long t=1;. This declares a new long integer variable named t and assigns it the value 1.
Next the test condition is evaluated; if this condition evaluates to TRUE, the statements enclosed in the curly braces are executed. In our case, the test condition is t<=a;. As long as t is less than or equal to a, we do the statements enclosed in the curly braces.
After the enclosed statements are executed, the advancing statement is executed; in our case, that statement is t++;. This means that t is increased by 1.
Now the entire process repeats. For our example, t now holds the value 2. As long as this is less than or equal to a, the enclosed statements are executed. Then t is advanced to 3, then to 4, and so forth, until the test condition is FALSE. At that point, the for loop is exhausted and we proceed to the next statement after the close brace; in our example that’s at line 37 and the statement is return d;.
In other words, the code
for (long t=1; t<=a; t++) { statements;
}
executes the statements between the curly braces a times with t equal to 1, then 2, then 3, and so on, until t equals a.
This style of for statement is common in programs. Of course, we can use the for statement to step down through values as in this code:
for (long s=n; s > 0; s--) { statements;
}
Alternatively, we could step only through odd values of the index:
for (long j=1; j <= n; j += 2) { statements;
}
For our gcd procedure, letting t take the values, 1, 2, 3, and so on, until a is precisely what we want. For each value of t, we check if t is a divisor of both a and b. We do this with the conditional if( (a%t==0) && (b%t==0) ). If this condition is satisfied, we update d to the current value of t, and this is what happens on line 33.
After the loop is finished, the value held in d is returned.
Greatest Common Divisor |
37 |
Now that our gcd program is finished, it’s time to try it out. In a separate file, that we name gcd-tester.cc, we write a simple main() procedure to try out our code.
Program 3.6: A program to test the gcd procedure.
1#include "gcd.h"
2 #include <iostream>
3using namespace std;
4
5/**
6* A program to test the gcd procedure.
7*/
8 int main() {
9long a,b;
10
11cout << "Enter the first number --> ";
12cin >> a;
13cout << "Enter the second number --> ";
14cin >> b;
15
16cout << "The gcd of " << a << " and " << b << " is "
17<< gcd(a,b) << endl;
18return 0;
19}
3.3Euclid’s method
The program we developed in Section 3.2 works, but it is a slow, inefficient algorithm. Suppose we want to calculate the greatest common divisor of numbers that are around one billion? This is not too large for a long integer, but in order to find the answer, the trial-division algorithm runs for billions of iterations. There is a much better way that was developed by Euclid.
The key idea is the following result.
Proposition 3.1. Let a,b be positive integers and let c = a mod b. Then gcd(a,b) = gcd(b,c).
Proof. Let a,b be positive integers and let c = a mod b; that is, a = qb + c where q,c Z and 0 ≤ c < b.
Note that if d is a common divisor of a and b, then d is also a divisor of c because c = a−qb. Conversely, if d is a common divisor of b and c, then (because a = qb+c) d is also divisor of a. Hence the set of common divisors of a and b is the same as the set of common divisors of b and c. Therefore gcd(a,b) = gcd(b,c). 
We use this to develop an algorithm. Suppose we want to find gcd(80,25). By Proposition 3.1, we calculate 80 mod 25 = 5 and so gcd(80,25) = gcd(25,5). To
38 |
C++ for Mathematicians |
find gcd(25,5), we again apply Proposition 3.1. Because 25 mod 5 = 0, we have gcd(25,5) = gcd(5,0). At this point Proposition 3.1 does not apply because 0 is not positive. However, we know that gcd(5,0) = 5 and so we have
gcd(80,25) = gcd(25,5) = gcd(5,0) = 5
and we only needed to do two divisions (not 25 or 80).
In this section we present two programs for computing the greatest common divisor by this method. Here is the first.
Program 3.7: A recursive procedure for gcd.
1#include "gcd.h"
2 #include <iostream>
3using namespace std;
4
5long gcd(long a, long b) {
6 // Make sure a and b are both nonnegative 7 if (a<0) a = -a;
8if (b<0) b = -b;
9
10// if a and b are both zero, print an error and return 0
11if ( (a==0) && (b==0) ) {
12cerr << "WARNING: gcd called with both arguments equal to zero."
13<< endl;
14return 0;
15}
16
17// If b is zero, the answer is a
18if (b==0) return a;
19// If a is zero, the answer is b
20if (a==0) return b;
21
22 long c = a%b;
23
24return gcd(b,c);
25}
Lines 7 and 8 ensure that a and b are nonnegative. We use a slightly different syntax for these if statements. When an if is followed by exactly one statement, the curly braces may be omitted. The single-line statement if (a<0) a = -a; is equivalent to this:
if (a<0) { a = -a;
}
Line 11 checks to see if the arguments are both zero; if they are, we issue a warning and return zero.
Lines 18 and 20 check if one of the arguments is zero; if so, the other argument is the answer we desire.
The real work of this procedure is on lines 22 and 24. If the program reaches line 22, then we know that both a and b are positive integers and Proposition 3.1
Greatest Common Divisor |
39 |
applies. We calculate c to be a mod b, so the answer to this problem is gcd(b,c) and we return that.
The question is: How is gcd(b,c) computed? The answer is that the gcd procedure calls itself. This is known as recursion. If we call gcd(80,25), then when we reach line 24, c holds the value 5. At this point we issue a call to gcd(25,5). The previous call to gcd(80,25) goes “on hold” pending the result of gcd(25,5). When this second invocation of gcd reaches line 24 we have a equal to 25, b equal to 5, and c equal to 0. A third call to gcd is generated at this point requesting gcd(5,0) and the second call is also placed on hold. During the call to gcd(5,0), we come to line 18 (because b==0 evaluates to TRUE) and so 5 is returned. This is passed back to the second and then to the first call to gcd and the final answer, 5, is returned.
Recursion is a powerful idea. When it applies, such as in this case, it can make for particularly simple code. The primary danger in using this technique is infinite recursion: if the program does not check adequate boundary conditions, it may run forever without returning an answer. This is akin to neglecting the basis case in a proof by induction.
A classic example of recursion is a program to calculate factorials. Here is the general idea presented in an incorrect program.
long factorial(long n) { return n*factorial(n-1);
}
To calculate factorial(5) the computer makes a call to factorial(4) which (if all goes well) returns the value 24. We then multiply this by 5 to get our answer.
The mistake, however, is that factorial(4) calls factorial(3) which calls factorial(2), and so on forever. We need to catch this process someplace, and the right place is when n is zero.
Here’s a better version.
long factorial(long n) { if (n==0) return 1; return n*factorial(n-1);
}
This gives the correct result for factorial(5), but is not entirely free of the infinite recursion trap. Figure out for yourself the problem and what you can do to address it. (See Exercise 3.2.)
Program 3.7 is a perfectly good, efficient gcd procedure. At this point, it would be proper to move on to the problem at hand (described at the beginning of this chapter). However, it is possible to make the program slightly more efficient and, in so doing, we can study another C++ feature: while loops.
There are some minor inefficiencies in Program 3.7. The tests to ensure that a and b are nonnegative and not both zero are run at every iteration. One can prove, however, that we do not need to worry about this in the embedded calls to gcd. The extra tests are performed unnecessarily. Also, the computer needs to do a modest amount of work every time a new procedure is called. To calculate the greatest common
40 |
C++ for Mathematicians |
divisor of two large numbers may result in dozens of calls to gcd. This overhead is not a serious problem, but if we plan to call gcd repeatedly, small improvements are worthwhile.
We now present another version of the gcd procedure that does not use recursion.
Program 3.8: An iterative procedure for gcd.
1#include "gcd.h"
2 #include <iostream>
3using namespace std;
4
5long gcd(long a, long b) {
6 // Make sure a and b are both nonnegative 7 if (a<0) a = -a;
8if (b<0) b = -b;
9
10// if a and b are both zero, print an error and return 0
11if ( (a==0) && (b==0) ) {
12cerr << "WARNING: gcd called with both arguments equal to zero."
13<< endl;
14return 0;
15}
16
17 long new_a, new_b; // place to hold new versions of a and b
18
19/*
20* We use the fact that gcd(a,b) = gcd(b,c) where c = a%b.
21* Note that if b is zero, gcd(a,b) = gcd(a,0) = a. If a is zero,
22* and b is not, we get a%b equal to zero, so new_b will be zero,
23* hence b will be zero and the loop will exit with a == 0, which
24* is what we want.
25*/
26
27while (b != 0) {
28new_a = b;
29new_b = a%b;
30
31a = new_a;
32b = new_b;
33}
34
35return a;
36}
The beginning of this program is the same as Program 3.7; we make sure the arguments are nonnegative and not both zero.
At line 17 we declare two new variables: new_a and new_b. The idea is to let
a0 ← b and b0 ← a mod b
and to continue with a0 and b0 in lieu of a and b. We keep doing this until b reaches zero, and then a will hold the answer. For example, starting with a = 80 and b = 25, we have this:
Greatest Common Divisor |
41 |
|||||
|
|
|
|
|
|
|
|
Iteration |
1 |
2 |
|
3 |
|
|
value of a |
80 |
25 |
|
5 |
|
|
value of b |
25 |
5 |
|
0 |
|
These steps take place on lines 27–33. The control statement is while. The general form of a while statement is this:
while (condition) { statements to perform;
}
When a program encounters a while statement, it checks to see if condition is TRUE or FALSE. If condition is FALSE, it skips all the statements enclosed in the curly braces. However, if condition is TRUE, then the statements inside the braces are executed and then we start the loop over again, checking to see if condition is TRUE or FALSE.
In this manner, the loop is executed over and over again until such time as the condition specified in the while statement is FALSE.
This is precisely what we need here. As long as b is not zero, we replace a and b with b and a%b, respectively. We enlist the help of the temporary variables new_a and new_b to do this.
When the while statement terminates, a holds the greatest common divisor of the original a and b, and so we return that at line 35.
It is interesting to note that we do not need two temporary variables for the while loop. The following works as well.
while (b != 0) { new_b = a%b; a = b;
b = new_b;
}
Although this is correct, it is more difficult to understand than the while loop in Program 3.8. The slightly more verbose version shows clearly how a and b are updated from the previous values of a and b. Clarity is often preferred to cleverness because clarity is more likely to be correct.
3.4Looping with for, while, and do
We have introduced two of C++’s looping statements: for and while. Here we introduce you to a third: do.
Recall that a while loop is structured like this:
while (condition) { statements to execute;
}
42 |
C++ for Mathematicians |
The condition is tested before the enclosed statements are ever executed. If condition is FALSE, the statements will never be executed. Incidentally, the while structure can be replaced by a for loop as follows.
for(; condition; ) { statements to execute;
}
The starting and advancing statements are missing, so they have no effect. We mention this mostly as a curiosity; the while version is preferable because it is clearer.
Sometimes you may find that the condition you want to check doesn’t make sense until some series of instructions has been performed at least once. For example, a program might read in data from a file and should stop reading when a negative value is encountered. We need to read at least one data value before the condition makes sense.
For such situations, the following structure can be used.
do {
statements to execute; } while (condition);
The statements to execute are performed at least once. At the end of the loop, the condition is checked. If it is TRUE, then we return to the start of the loop for another round; if the condition is FALSE, the loop is finished and control passes to the next statement after the loop. The do-while loop structure is not used as frequently as while and for loops.
There are two special statements that can be used to modify the execution of a loop: break and continue.
The break command causes the loop to exit immediately with control passing to the next statement after the loop. Consider this code.
for (long a=0; a<100; a++) { statement1;
statement2;
if (x > 0) break; statement3; statement4;
}
statement5;
If after statement2 the variable x holds a positive quantity, we skip statement3 and statement4, and go directly to statement5.
The break command can be used with while and do loops, too.
The continue command directs the computer to go to the end of the loop and then do what is natural at that point. Consider this code.
for (long a=0; a<100; a++) { statement1;
statement2;
if (x > 0) continue;
Greatest Common Divisor |
43 |
statement3;
statement4;
}
statement5;
In this case, all 100 iterations of the loop take place. However, if during some iteration, after we execute statement2, we find x to hold a positive quantity, we skip both statement3 and statement4. At this point we increase a (because the advancing statement is a++). If a is less than 100, another pass through the loop is performed.
3.5An exhaustive approach to the GCD problem
It is now time to tackle the problem from the beginning of this chapter. Let pn be the probability that two integers, chosen independently and uniformly from {1,2,...,n}, are relatively prime.
The following program counts the number of pairs (a,b) with 1 ≤ a,b ≤ n that satisfy gcd(a,b) = 1 and divides by n2.
Program 3.9: A program to calculate pn.
1 #include <iostream>
2#include "gcd.h"
3using namespace std;
4
5/**
6* Find the probability that two integers in {1,...,n} are relatively
7* prime.
8 |
*/ |
9 |
|
10int main() {
11long n;
12cout << "Enter n --> ";
13cin >> n;
14
15 long count = 0;
16
17for (long a=1; a<=n; a++) {
18for (long b=1; b<=n; b++) {
19if (gcd(a,b) == 1) {
20count++;
21}
22}
23}
24
25cout << double(count) / (double(n) * double(n)) << endl;
26return 0;
27}
44 |
C++ for Mathematicians |
The code is straightforward. We ask the user to enter n and set a counter (named count) to zero. The main action takes place in lines 17–23. Two nested for loops run through all possible values of a and b with 1 ≤ a,b ≤ n. On each iteration, if a and b are found to be relatively prime, the counter is increased by one.
At the end, we divide the counter by n2 to get the probability. Because we want a floating point answer, we cast the appropriate terms into type double (see line 25).
To run this program, we need to use three files: gcd.h, gcd.cc, and this file (let’s call it exhaust.cc). We can either load all three into an IDE or (on a UNIX system) save all three in a directory and compile with a command such as
g++ gcd.cc exhaust.cc -o exhaust
and the program is saved in a file named exhaust. (See Appendix A for more details.)
Here is a typical run of the program.
Enter n --> 1000.6087
Thus, p100 = 0.6087.
If we run the program for various values of n, we generate the following results.
n |
pn |
100 |
0.6087 |
500 |
0.608924 |
1000 |
0.608383 |
5000 |
0.608037 |
10000 |
0.60795 |
50000 |
0.607939 |
It appears that limn→∞ pn exists and is converging to a value around 0.6079. It would be useful to extend this table further. However, the last entry in this table took a long time for my computer to calculate. Is there a more efficient method? There is one modest modification we can make to the program. Notice that we are computing the same values twice in most instances. That is, we compute both gcd(10,15) and gcd(15,10). We can make our program twice as fast by calculating only one of these. Also, we do not have to bother calculating gcd(a,a); the only instance in which gcd(a,a) equals one is when a = 1.
Here is the modified version of the program.
Program 3.10: A slightly better program to calculate pn.
1 #include <iostream>
2#include "gcd.h"
3using namespace std;
4
5/**
6* Find the probability that two integers in {1,...,n} are relatively
7 * prime.
8*/
Greatest Common Divisor |
45 |
9
10int main() {
11long n;
12cout << "Enter n --> ";
13cin >> n;
14
15 long count = 0;
16
17for (long a=1; a<=n; a++) {
18for (long b=a+1; b<=n; b++) {
19if (gcd(a,b) == 1) {
20count++;
21}
22}
23}
24count = 2*count + 1;
25
26cout << double(count) / (double(n) * double(n)) << endl;
27return 0;
28}
Notice that the second for loop (line 18) begins with b=a+1, so the inner loop runs from a+1 up to n. Thus the program calculates gcd(5,10), but not gcd(10,5). We need to double count at the end to correct. Also, we do not calculate gcd(t,t) for any t, so we add 1 to (the doubled) count at the end to correct for that. These corrections are on line 24.
In the following chapters, we approach the problem using two other methods: a Monte Carlo randomized algorithm (giving us an opportunity to learn about random numbers in C++) and a more intelligent exhaust using Euler’s totient ϕ (giving us an opportunity to learn about arrays). However, before we move on to those other approaches, we make our gcd procedure better still and use that opportunity to learn additional C++ features.
3.6Extended gcd, call by reference, and overloading
The Euclidean algorithm can be used not only to find the greatest common divisor of two integers a and b, but also to express the gcd as an integer linear combination of a and b. That is, if a and b are integers (not both zero) then there exist x,y Z such that ax + by = gcd(a,b).
Our goal is a procedure that can be called like this: d = gcd(a,b,x,y);
where a,b are the numbers whose gcd we desire, d is gcd(a,b), and x,y have the property that d = ax + by. At first glance this appears impossible for two reasons. First, we already have a procedure named gcd. Shouldn’t we name this something