C++ For Mathematicians (2006) [eng]

46	C++ for Mathematicians

else (such as extended_gcd)? Second, we can pass values to a procedure in its list of arguments, but the procedure cannot change the value of the arguments. This fact is known as call by value and was discussed earlier in this chapter (see Section 3.2). The good news is that neither of these is a significant hurdle.

First, two different procedures may have the same name, but there is a caveat: The procedures must have different types of arguments. For example, suppose we want to create procedures for (a) finding the slope of a line segment from the origin to a given point and (b) finding the slope of a line segment joining two given points. In some programming languages, you would be required to give these different names; in C++, we may name them both slope. The declarations for these procedures (which we put in a header file named, say, slope.h) are these.

double slope(double x, double y);

double slope(double x1, double y1, double x2, double y2);

The first is for the origin-to-slope version of the procedure and the second for the point-to-point version. Because the first takes two double arguments and the second takes four double arguments, the C++ compiler recognizes these as different procedures.

The definitions of the procedures (which we place in a file named slope.cc) look like this:

double slope(double x, double y) { return y/x;

}

double slope(double x1, double y1, double x2, double y2) { return (y1-y2)/(x1-x2);

}

When we use the same name for different procedures, these procedures ought to be closely related and serve essentially the same purpose.

We mathematicians often use the same symbol for two closely related (but different) objects. For example, if we have a function f : R → R, then f (x) only makes sense when x is a real number. However, we often extend this notation. If A R, then f (A) (same f ) means { f (x) : x A}.

In C++, we refer to this ability to name different procedures with the same name as overloading.

Second, it is true that C++ passes data to procedures using call by value; this means that the arguments to the procedure are copies of the original, and the invoked procedure cannot change the original values. For example, consider this procedure.

void swap(long a, long b) { long tmp;

tmp = a; a = b; b = tmp;

}

The procedure’s return type is void. This means that the procedure does not return any value. The code in swap takes the values held in a and b and exchanges them.

Greatest Common Divisor

47

If, when the procedure is called a holds 5 and b holds 17, then at the end, a holds 17 and b holds 5. However, this procedure has no effect whatsoever on the calling procedure. Suppose the calling procedure contained this code:

long a = 5; long b = 17; swap(a,b);

cout << "a = " << a << endl; cout << "b = " << b << endl;

The output of this would be as follows.

a = 5b = 17

There is a way to modify the procedure swap so that it is capable of modifying its arguments. Instead of passing the value of a to the procedure, we pass a reference to a. To do this, we add the character & to the end of the type name, like so:

void swap(long& a, long& b) { long tmp;

tmp = a; a = b; b = tmp;

}

The & does this: The argument a is a reference to a variable of type long. Not only is a of type long, it is more than a copy of the long. This syntax means: instead of sending a copy of a to this procedure, send the variable itself (and not a clone). When a procedure receives a reference to a variable, it is working on the original and not a copy.

With the rewritten swap procedure, the code

long a = 5; long b = 17; swap(a,b);

cout << "a = " << a << endl; cout << "b = " << b << endl;

		produces this output:
		a = 17
		b = 5

Call by reference is the mechanism we need so that the new gcd procedure can deliver three answers: the gcd returned via a return statement and the values x and y. The declaration for the procedure (which we add to gcd.h) looks like this:

long gcd(long a, long b, long& x, long& y);

A recursive approach to solving the extended gcd problem works well. Suppose we want to solve the extended gcd problem for positive integers a and b. We start by

48	C++ for Mathematicians

solving the extended gcd problem for b and c where c = a mod b. Let’s say we have that solution so

d = gcd(b,c) = bx0 + cy0

for some integers x0,y0. We know that c is the remainder when we divide a by b; that is,

a = qb + c

where 0 ≤ c < b. Writing c = a −qb we have

d= bx0 + cy0

=bx0 + (a −qb)y0

=ay0 + b(x0 −qy0)

=ax + by

where x = y0 and y = x0 −qy0.

These ideas form the heart of the recursion. To construct the procedure, we need to check special cases (a or b might be zero or negative). The code follows.

Program 3.11: Code for the extended gcd procedure.

1 long gcd(long a, long b, long &x, long &y) {

2long d; // place to hold final gcd

3
4	// in	case	b = 0, we have d=|a|, x=1 or -1, y arbitrary (say, 0)
5	if (b==0) {
6	if (a<0)		{
7	d	= -a;
8	x	= -1;
9	y	= 0;

10}

11else {

12d = a;

13x = 1;

14y = 0;

15}

16return d;

17}

18

19// if b is negative, here is a workaround

20if (b<0) {

21d = gcd(a,-b,x,y);

22y = -y;

23return d;

24}

25

26// if a is negative, here is a workaround

27if (a<0) {

28d = gcd(-a,b,x,y);

29x = -x;

30return d;

31}

Greatest Common Divisor

49

32

33// set up recursion

34long aa = b;

35long bb = a%b;

36long qq = a/b;

37long xx,yy;

38

39 d = gcd(aa,bb,xx,yy);

40

41x = yy;

42y = xx - qq*yy;

43

44return d;

45}

3.7Exercises

3.1Use the Euclidean algorithm to find d = gcd(51,289) and to find integers x and y so that d = 51x + 289y.

3.2On page 39 we presented the following mostly correct procedure for computing factorials.

long factorial(long n) { if (n==0) return 1; return n*factorial(n-1);

}

This procedure, however, still has a bug. What is it?

3.3 Write a procedure to calculate the sign of its argument (the signum function). That is, your procedure should calculate

+1 if x > 0,

sgnx = 0 if x = 0, and

−1 if x < 0.

3.4Write two procedures for producing the nth Fibonacci number when the number n is given as input. In one procedure, use a for loop and produce the answer using iteration. In the second version, use recursion. The input argument and return value should be type long.

Use the following definition of Fibonacci numbers: F0 = F1 = 1, and Fn = Fn−1 + Fn−2 for all n ≥ 2.

Have your procedure return −1 if the input argument is negative. Don’t worry about overflow.

50	C++ for Mathematicians

3.5Write a main to test the Fibonacci procedures you created in Exercise 3.4. Use it to evaluate F20, F30, and F40.

You should find that one version is much faster than the other. Explain why. Include in your explanation an analysis of how many times the recursive version of your procedure gets called in order to compute Fn.

3.6Write two procedures to calculate

	N	1
	f (N) = ∑		.
		2
	k=1	k

Note that for N large, this approaches ζ (2) = π2/6.

The first procedure should calculate the sum in the usual order

	1		1	1
1 +		+		+ ···+
	4		9		N2

and the second should calculate the sum in the reverse order

1		+	1	+ ···+	1	+ 1.
	N2		(N −1)2		4

In both cases, use float variables for all real values.

Evaluate these sums for N = 106 and report which gives the better result. (π2/6 ≈ 1.6449340668482264365.)

Explain.

3.7Write procedures for converting between rectangular and polar coordinates named xy2polar and polar2xy. Invoking xy2polar(x,y,r,t); should take the rectangular coordinates (x,y) and save the resulting polar coordinates; (r,θ) are r and t, respectively. The procedure polar2xy(r,t,x,y); should have the reverse effect. Be sure to handle the origin in a sensible manner.

Of course, you need trigonometry functions to accomplish this task; consult Appendix C (especially Appendix C.6.1) where you can find useful functions such as atan2.

3.8Every positive integer has a multiple that, when expressed in base ten, is comprised entirely of zeros and ones. Write a procedure to find the least multiple of n, the input parameter, of the required form.

Warnings: If the long long type1 is available on your computer, use it. The least multiple of the required form may be much larger than n. For example, what is the least multiple of 9 that contains only zeros and ones? What is the least such multiple of 99?

Design your procedure to detect overflow and return −1 if it is unable to find the required multiple.

1Or int64.

Greatest Common Divisor

51

3.9 Write a procedure to find the gcd of three integers of the form

long gcd(long a, long b, long c);

and an extended version of the form

long gcd(long a, long b, long c, long& x, long& y, long& z);

that returns d = gcd(a,b,c) and populates the last three arguments with values x,y,z such that ax + by + cz = d.

You may use the gcd procedures developed in this chapter as part of your solution.

Chapter 4

Random Numbers

In this chapter we discuss the generation random numbers in C++. The motivation is the problem introduced in Section 3.1: let pn be the probability that two numbers, sampled independently and uniformly from {1,2,...,n} are relatively prime. What can we say about pn as n → ∞?

In Chapter 3 we used direct enumeration to calculate pn for various values of n. However, as n approaches 100,000, the time for the computer to complete the calculation becomes excessive. This motivates us to find another attack. The approach we take is to sample pairs of integers in {1,2,...,n} at random and keep track of how often we find that they are relatively prime. To do this, however, we require a mechanism to generate (pseudo) random numbers.

4.1Pseudo random number generation

We need a procedure that produces uniform random values in the set {1,2,...,n}. In this section, we develop such a procedure as well as a procedure to produce random real values in an interval.

Before we begin, however, we need to understand that the computer does not produce random values. Instead, it provides a pseudo random number generator. This is a procedure whose output looks random, but is actually deterministic. The simplest type of random number generator is a linear congruential generator or LCG for short. An LCG produces a series of integers x0,x1,x2,... by the following calculation,

xn+1 = (axn + b) mod c

where a,b,c are fixed positive integers. The method is specified fully once we choose a value for x0. This first value is known as the seed.

Pseudo random numbers produced by an LCG behave in some ways as uniform random values from {0,1,...,c − 1}, but there are problems. For example, by the pigeonhole principle, the sequence x0,x1,x2,... must repeat itself over and over. The hope is that by taking c sufficiently large, this is not an issue.

A variety of more sophisticated pseudo random number generators have been proposed and implemented each with various advantages and disadvantages. The important point we need to keep in mind is that pseudo random numbers are not random,

53

54	C++ for Mathematicians

so we need to be mildly skeptical of the results they suggest and careful in how we use them.

Standard C++ includes the rand procedure for producing pseudo random numbers. To use rand it is necessary to include the cstdlib library. A call to rand() returns an integer value between 0 and a constant named RAND_MAX (inclusive). The value of RAND_MAX is defined in the cstdlib header file. A typical value is 231 −1.

4.2Uniform random values

It is possible to write programs so they call rand directly whenever a random value is required. However, it is a good idea to create our own procedures that serve as intermediaries between the program that requires a random value and rand. Why do we need a middle man? Here are three good reasons.

•First, rand produces random values in a large discrete set. We might want a random integer in a smaller set (such as {0,1} if we want to simulate a coin flip) or we might want a uniform continuous value in the interval [0,1]. It is convenient to have procedures that do each of these.

•Second, there are preferred methods to extract random values from rand. For example, if we want to simulate a coin flip, it is tempting to write code such as this: flip = rand()%2;. The problem is that the lowest-order bit of random values from a pseudo random number generator might not behave as you would expect. In the worst case, this last bit might simply oscillate between 0 and 1, and such coin flips do not look random at all. We create a better alternative one time and use that procedure henceforth.

•Some day you might decide that you prefer a different random number generator. Rather than editing all your programs to replace rand with the new procedure, you only need to rewrite the intermediaries.

What do we want these procedures to do for us? We want to produce a random real value uniformly in the interval [0,1], or more generally in an interval [a,b] where a,b R. And we want to produce a random integer chosen uniformly from a finite set of the form {1,2,...,n}.

Because all of these return a random value that is uniformly distributed over its domain, we name all three of these unif. This procedure name overloading is permissible because the three versions have different type arguments.

Program 4.1: Header file, uniform.h for procedures to generate uniform random values.

1 #ifndef UNIFORM_H

2 #define UNIFORM_H

Random Numbers

55

3

4/**

5* Generate a random number between 0 and 1.

6* @return a uniform random number in [0,1].

7*/

8double unif();

9

10/**

11* Generate a random number in a real interval.

12* @param a one end point of the interval.

13* @param b the other end point of the interval.

14* @return a uniform random number in [0,1].

15*/

16double unif(double a, double b);

17

18/**

19* Generate a random integer between 1 and a given value.

20* @param n the largest value this procedure can produce.

21* @return a uniform random value in {1,2,...,n}.

22*/

23long unif(long a);

24

25/**

26* Reset the random number generator based on the system clock.

27*/

28void seed();

29

30 #endif

The three flavors of unif are declared on lines 8, 16, and 23. The first produces a uniform real value in the interval [0,1], the second generalizes this and produces a real value in an arbitrary interval [a,b], and the third produces an integer value uniformly in {1,2,...,n}.

In addition, we have declared a procedure named seed on line 28. The purpose of seed is to initialize the random number generator from the system clock (more on this later).

Let’s turn to implementing each of these starting with double unif(). A call to rand() returns an integer between 0 and RAND_MAX. To convert this to a real value in [0,1] we simply compute rand() / double(RAND_MAX);. (See lines 7–9 of the file uniform.cc in Program 4.2 below.)

Once we have the double unif() version written, we use it to write the second version. We simply multiply unif() by (b-a) and add a. See lines 11–13 in uniform.cc.

For the integer version, we multiply the continuous unif() by a and convert to an integer. This gives a value in {0,...,a − 1}, so we add 1 to place the value into the desired set. What if the user gives a negative or zero value for a? We have a choice as to handling these undefined situations. See lines 15–19 of uniform.cc to see our decision.

Now we consider the seed procedure. Because rand (and the procedures that we wrote based on rand) returns an unpredictable stream of values, we might ex-