C++ For Mathematicians (2006) [eng]

56	C++ for Mathematicians

pect that the programs we write would behave differently every time we run them. Interestingly (and with good reason) this is not the case. Every time a program containing rand is run, the rand procedure gives the same sequence of pseudo random values. The reason is that the first time rand is invoked, it contains a fixed starting value called the seed of the random number generator. The reason this behavior is desirable is reproducibility. If you perform a computational experiment and wish to report it in a journal, it is important that others can run the same program as you and see the same results.

However, there are times when this reproducible behavior is undesirable. For example, the motivating purpose of this chapter is to write code to generate many pairs of integers in a set {1,...,n} to see how frequently they are relatively prime. We might want to run our program a few times with different streams of random values so we can compare results.

The standard library (cstdlib) provides the procedure srand that is used to set the seed used by rand. Calling srand(s), where s is a long integer, resets to the seed to s. We could ask the user to provide a seed value like this:

long s;

cout << "Enter a seed for the random number generator --> "; cin >> s;

srand(s);

Another solution, that is easier for the user, is to use the computer’s clock to provide the seed. As long as we don’t run the program twice in the same second, a different value is used for the seed.

The header ctime defines the procedure time. Calling time(0) returns the number of seconds that have elapsed since a specific date and time (on many computers, that date and time is January 1, 1970 at 12:00 A.M. UTC). Our seed procedure simply takes the value1 returned by time(0) as input to srand. See lines 21–23 of uniform.cc.

Program 4.2: Definitions of the unif procedures in uniform.cc.

1 #include "uniform.h"

2 #include <cstdlib>

3 #include <ctime>

4#include <cmath>

5using namespace std;

6

7double unif() {

8 return rand() / double(RAND_MAX);

9}

10

11double unif(double a, double b) {

12return (b-a)*unif() + a;

13}

1To be precise, the procedure time returns a value of type time t. Your compiler might require you to convert this to an unsigned integer before using it as an argument to srand. Replace line 22 with srand(unsigned(time(0))); in this case.

Random Numbers

57

14

15long unif(long a) {

16if (a < 0) a = -a;

17if (a==0) return 0;

18return long(unif()*a) + 1;

19}

20

21void seed() {

22srand(time(0));

23}

4.3 More on pseudo random number generation

The integer version of unif used the expression long(unif()*a)+1 to produce a random value in {1,...,a}. This is mildly convoluted because we have the additional call to the continuous version of unif that calls rand and divides by RAND_MAX. It is both simpler and more efficient simply to calculate 1 + rand()%a, but this latter approach is less reliable.

The problem is that some older versions of the rand pseudo random number generator are purported to be unreliable, and the lower-order bits of the values produced by rand do not behave well. One way to rectify this situation is to replace rand by a better procedure.

To understand why the low-order bits of a random number generator might not be good approximations of randomness, we develop our own linear congruential generator (LCG) here. (And this gives us an opportunity to introduce additional C++ concepts.)

Recall that an LCG produces a stream of values x0,x1,x2,... where

xn+1 = (axn + b) mod c.

The LCG is fully specified once we select values for a, b, c, and x0. For our example, we take

a = 17, b = 3, c = 64, and x0 = 0.

Suppose we call our procedure lcg. Each time we call lcg() it should return the next xj value in the sequence. However, for this the procedure lcg needs to remember the previous x value. The following code does not work.

int lcg() {

int state = 0;

return (17*state + 3) % 64;

}

Every time this procedure is called, it returns the value 3. We need to indicate that the variable state should only be initialized to zero the first time lcg is called.

58	C++ for Mathematicians

To do this, we declare the variable state to be static. The declaration looks like this instead: static int state = 0;. With this declaration, the variable state retains its value after the procedure lcg ends. (Without the static modifier, the usual behavior is for variables to cease to exist once the procedure ends.) The variable state is initialized to zero only the first time lcg is called. Henceforth, it retains the value it held when lcg last terminated.

The return statement is fine as written (return (17*state + 3) % 64) but there is a better way. The constants 17, 3, and 64 should be given names and declared at the beginning of the procedure. For a short simple procedure such as lcg, this is not an important issue. However, for more complicated procedures, giving constants specific names makes the program easier to read and easier to modify. Imagine that we write a program in which we are often reducing numbers modulo 64. Rather than explicitly typing 64 repeatedly in the code, we can define a variable named theMod set equal to 64 instead. Although typing theMod is longer than typing 64, if we ever want to change the program so that theMod is now, say, 128, we only have to change one line. The declaration of theMod looks like this:

const int theMod = 64;

The const modifier means that the procedure does not modify the value of theMod. This enables the C++ compiler to generate more efficient object code and prevents you from mistakenly putting theCode on the left-hand side of an assignment statement.

Returning to the lcg procedure, we declare three variables (named a, b, and c) to be type const int (see lines 12–14 of Program 4.3). The return statement is then return (a*state+b) % c; which is easier to read. Later, if we wish to modify the lcg program, we simply need to change the values we ascribe to a, b, or c at the top of the procedure.

With the lcg procedure written, we create a short main() to test it out. The main() calls lcg 20 times, reducing each return value modulo two. It then calls lcg another 20 times, reducing each of those results mod four. Here is the full program.

Program 4.3: A program to illustrate the problem with lower-order bits in an LCG.

1 #include <iostream>

2using namespace std;

3

4/**

5* A sample linear congruential pseudo random number generator that

6* returns values in {0,1,...,63}.

7*/

8int lcg() {

9 static int state = 0;

10const long a = 17;

11const long b = 3;

12const long c = 64;

13

Random Numbers

59

14state = (a*state+b) % c;

15return state;

16}

17

18/**

19* This main calls lcg twenty times and prints out the value modulo

20* two, and then prints twenty more values taken modulo four.

21*/

22int main() {

23cout << "Values mod 2: ";

24for (int k=0; k<20; k++) {

25cout << lcg()%2 << " ";

26}

27cout << endl;

28

29cout << "Values mod 4: ";

30for (int k=0; k<20; k++) {

31cout << lcg()%4 << " ";

32}

33cout << endl;

34

35return 0;

36}

		When this program is run, the following output is printed on the screen.
		Values mod 2:	1 0	1 0	1 0 1 0 1	0 1 0 1 0	1 0 1 0 1	0
		Values mod 4:	3 2	1 0	3 2 1 0 3	2 1 0 3 2	1 0 3 2 1	0

Clearly, the sequence of values we produce in this manner is far from random! However, we can change the way we extract zeros and ones from lcg. Consider the following alternative main.

int main() {

cout << "Values mod 2: "; for (int k=0; k<20; k++) { double x = lcg() / 64.; cout << int(2*x) << " ";

}

cout << endl; return 0;

}

Here we first produce a double value x in the range [0,1) by dividing the output of lcg by 64. We then multiply x by 2 (so we are now somewhere in [0,2)) and thencast to type int (so the value is now either 0 or 1). Here’s the result.

Values mod 2: 0 1 0 1 1 1 0 0 0 0 1 0 0 1 0 1 1 1 0 0

As you can see, the output appears, at least superficially, more random.

60	C++ for Mathematicians

4.4A Monte Carlo program for the GCD problem

We now present a simple program to estimate pn. The user selects n (thereby specifying the set {1,2,...,n} from which pairs of integers are to be drawn) and the number of repetitions. The program generates the pairs, counts the number that are relatively prime, and reports the frequency.

Program 4.4: A Monte Carlo approach to calculating pn.

1 #include "uniform.h"

2#include "gcd.h"

3 #include <iostream>

4using namespace std;

5

6/**

7* This main generates many pairs of values from the set {1,2,...,n}

8* and reports how often the pairs are relatively prime. The value n

9* and the number of pairs are specified by the user.

10	*/
11
12	int main() {
13	long n;	// max el’t in the set {1,2,...,n}
14	long reps;	// number of times we perform the experiment
15	long a,b;	// values chosen from {1,2,...,n}
16	long count;	// number of pairs that are relatively prime
17
18	count = 0;
19

20cout << "Enter n (maximum el’t of the set) --> ";

21cin >> n;

22

23cout << "Enter the number of pairs to sample --> ";

24cin >> reps;

25

26for (long k=1; k<=reps; k++) {

27a = unif(n);

28b = unif(n);

29if (gcd(a,b) == 1) ++count;

30}

31

32 cout << count / (double(reps));

33

34return 0;

35}

To begin, let’s run the program with n = 1000 for 10,000 repetitions. The session

		looks like this:
		Enter n (maximum el’t of the set) --> 1000
		Enter the number of pairs to sample --> 10000
		0.6117

Random Numbers

61

The program estimates p1000 ≈ 0.6117 which is reasonably close to the actual value 0.608383. If, however, we run the code for a million repetitions we get p1000 ≈ 0.608932 which is correct to three decimal places.

We can use the program to estimate pn for n = 106. Running the program for 108 repetitions, we find pn ≈ 0.607979. How good is this estimate? When we perform

r independent repetitions of an experiment whose probability of success is p, the

√

expected number of successes is rp, but the standard deviation is on the order of r.

√

In this case, r = 108, so r = 104. So it is reasonable to suppose that the value of pn we found is correct to four decimal places. To obtain greater accuracy, we need either to increase the number of repetitions or (as we do in the next chapter) find a better method for calculating pn.

One hundred million repetitions take a modest amount of time. It would be feasible to run for a billion repetitions, but that would only give us another “half” digit of accuracy. We might consider running for a trillion repetitions, but that would take too long.

Here’s a reasonable rule of thumb for today’s computers. Programs that do millions of operations using millions of bytes of memory are quick. Programs that do billions of operations take some time, but are feasible to run. Holding billions of bytes in memory, however, is at the limit of today’s personal computers. However, trillions of operations takes too long and most computers cannot hold trillions of bytes of data in memory.

4.5Normal random values

Before we leave the subject of random numbers, we expand our repertoire of the types of random variables we can simulate. In Section 4.2 we developed code to produce discrete and continuous uniform random variables. Here, we produce normal (Gaussian) random values.

The distribution of a normal random variable is based on the classic bell curve, f (t) = exp(−t2). We need to modify the bell curve slightly so that the total area under the curve is 1, the mean of the Gaussian is 0, and the standard deviation is 1. To that end, we use the following for the density.

	1			exp −t2/2 .
	f (t) =	√
			2π
From this it follows (with a bit of work) that					f (t) dt = 1,	t f (t) dt = 0, and
R	t2 f (t) dt = 1 (integration over all of R).				R	R
	We therefore define the Gaussian random variable X using					f as its density, or
equivalently
	Φ(x) = P[X ≤ x] = Z x				f (t) dt.

−∞

62	C++ for Mathematicians

The integrals at the end of the previous paragraph imply that X has mean zero and standard deviation one; for this reason, X is also known as N(0,1)—a normal random variable with mean zero and standard deviation one.

There is an efficient algorithm for producing Gaussian random values known as the polar or Box–Muller method. One begins by generating a point (x,y) uniformly at random in the unit disc. We then let

r = x2 + y2, µ = r		−2r	, Z1 = µx, and Z2 = µy.
		logr

Then Z1 and Z2 are independent N(0,1) Gaussian random variables. Here is a C++ program based on this method:

#include <cmath> #include "uniform.h" using namespace std;

double randn() { double x,y,r; do {

x = unif(-1.,1.); y = unif(-1.,1.); r = x*x + y*y;

} while (r >= 1.);

double mu = sqrt(-2.0 * log(r) / r);

return mu*x;

}

The do/while loop generates points (x,y) uniformly in the square [−1,1]2 until one that is interior to the unit disk is found. Each pass through the loop has a π/4 chance of succeeding, so after just a few iterations we are assured of finding a point chosen uniformly from the unit disk. Once the point (x,y) has been found, the rest of the algorithm follows the Box–Muller method.

There is an inefficiency in this implementation. The algorithm is capable of producing two independent normal random values. In our implementation we find one of these values and simply ignore the other. We could make this procedure twice as fast if there were a way to preserve that second normal value for the next call to this procedure. To do this, we enlist the help of static variables.

The header file randn.h looks like this:

#ifndef RANDN_H #define RANDN_H

#include "uniform.h"

double randn();

#endif

Here is the program file (randn.cc).

Random Numbers

63

Program 4.5: A program to generate Gaussian random values.

1 #include "randn.h"

2#include <cmath>

3using namespace std;

4

5double randn() {

6 static bool has_saved = false;

7static double saved;

8

9if (has_saved) {

10has_saved = false;

11return saved;

12}

13

14double x,y,r;

15do {

16x = unif(-1.,1.);

17y = unif(-1.,1.);

18r = x*x + y*y;

19} while (r >= 1.);

20

21 double mu = sqrt(-2.0 * log(r) / r);

22

23saved = mu*y;

24has_saved = true;

25

26return mu*x;

27}

The new version includes two static variables: a Boolean value has_saved and a real value saved. The has_saved variable is set to false to show that the procedure does not currently hold a saved Gaussian value (in saved). See lines 6–7.

At line 9 we check if there is a saved Gaussian value that we have not used yet. If so, then we change has_saved to false and return the value in saved.

Otherwise (starting at line 14), we generate the two Gaussian values Z1 and Z2 by the polar method. We save Z2 in saved and set the flag has_saved to true. Finally, we return Z1.

The statements on lines 14–21 are the time-consuming part of this procedure. By saving the second Gaussian value generated by the polar method, the slow part is only executed on every other invocation of the procedure. This makes the procedure nearly twice as fast.

4.6Exercises

4.1Suppose two points are chosen uniformly at random within the unit square [0,1]2. Write a program to estimate the expected (average) length of the seg-

64	C++ for Mathematicians

ment joining such points.

4.2Buffon’s Needle. Imagine a needle is dropped at random onto a floor painted with equally spaced parallel lines. The length of the needle is the same as the distance between the lines.

Write a program that simulates dropping the needle and that counts the number of times the needle crosses one of the lines drawn on the floor.

If you wish to use standard trigonometry functions, such as cos, or the floor function, insert a #include <cmath> directive at the start of your program. See Appendix C.6.1.

4.3Sylvester’s Four-Point Problem. Let K be a compact convex subset of the plane with nonempty interior. Let P(K) denote the probability that when four points are chosen independently and uniformly in K, then they lie at the vertices of a convex quadrilateral (as in the left portion of the illustration but not the right).

Write procedures to (a) generate a point uniformly at random inside a circle,

(b) generate a point uniformly at random inside a triangle, and (c) test whether four points determine the vertices of a convex quadrilateral.

Use your procedures to estimate P(K) for the cases where K is a circle or a triangle.

4.4Random point on a circle. Suppose you wish to generate a point uniformly at random on the unit circle, {(x,y) : x2 + y2 = 1}. Suppose the procedure is

declared like this:

void point_on_circle(double& x, double& y);

Here are two ways you might implement this procedure.

First, you can generate a uniform [0,1] value that you multiply by 2π; call the result θ. Then set x and y to be cosθ and sinθ.

Alternatively, you can generate a point in the interior of the unit ball by the rejection method. [That is, pick points (x,y) [0,1]2 until you find one that

p satisfies x2 + y2 ≤ 1.] Then rescale by dividing by x2 + y2.

Which do you think is faster? Create two versions of the procedure and time how long it takes each to generate 100 million points.

Random Numbers

65

4.5Random point on a sphere. Continuing Exercise 4.4, we consider the issue of generating a point at random on the surface of a sphere in R3 and in higher

dimensions.

One way to generate a point on the unit sphere in R3 is to generate x,y,z

p

uniformly in [−1,1]. If x2 + y2 + z2 ≤ 1, then we scale by 1/ x2 + y2 + z2 to generate the point. Otherwise, we generate another triple and try again.

What is the probability of successfully generating a point during a given iteration?

Next, explain why this procedure is woefully inefficient in high-dimensional space. Try to create an alternative method.

4.6Create a procedure int random_walk(); that simulates a random walk on the integers. That is, random_walk() returns the position of a particle whose initial position is 0. Each time random_walk() is called, the particle moves one step left or right, each with probability 50%. The return value is the new position of the particle.

		For example, calling random_walk 20 times might produce the following
		values.
	1 0 1 2 3 4 5 4 5 4 3 4 3 4 3 4 3 2 3 2

4.7Write a pair of procedures int up() and int down() that behave as follows. When up() is called, a certain value is increased by 1 and the new value is returned. Similarly, when down() is invoked, the value is decreased by 1 and that new value is returned. At the start of the program, the value is zero.

For example, suppose we run the following code.

cout << up() << " "; cout << up() << " "; cout << up() << " "; cout << down() << " "; cout << up() << " "; cout << down() << endl;

		Then the following output is produced.
			2
	1 2 3 2 3