C++ For Mathematicians (2006) [eng]

86	C++ for Mathematicians

19long* primes;

20primes = new long[TABLE_SIZE];

21long np;

22np = sieve(2*N,primes);

23

24	long long count=0;	// sum of phi(d) from 1 to n
25

26cout << setprecision(20);

27for (long k=1; k<=N; k++) {

28count += totient(k, primes);

29if (k%100000==0) {

30cout << k/1000 << " thousand \t";

31cout << double(2*count-1) / (double(k)*double(k)) << endl;

32}

33}

34return 0;

35}

Notice there is a new header included at line 4. The iomanip header provides devices to change the output style. In our case, we want to print out more digits of pn than usual. This occurs on line 26. The cout << setprecision(20); statement modifies cout so that it prints up to 20 decimal digits for double real numbers (trailing zeros, if any, are not printed). The iomanip header is needed to define setprecision. (See Section 14.6 for more information on adjusting output format.)

Lines 19–22 are used to generate the table of primes.

The variable count, declared on line 24, is used to accumulate the sum of ϕ(k). As discussed, this needs to be type long long because the final sum exceeds the maximum value a long can hold.

The core of the program is on lines 27–33. This is a for loop that increments count by ϕ(k) as k goes from one to one million. On line 29 we check if k is divisible by 100 thousand; if so, we report pk at that time.

Here is the output of the program (which took about 25 minutes to run on my

		computer).
		100 thousand		0.6079301507
		200	thousand	0.60792994587500004
		300	thousand	0.60792774407777783
		400	thousand	0.60792759136874996
		500	thousand	0.607928317404
		600	thousand	0.6079276484527778
		700	thousand	0.60792730424285712
		800	thousand	0.60792796007343752
		900	thousand	0.6079273649074074
		1000 thousand		0.60792710478300005

Arrays

87

5.8The answer

It certainly appears that pn is converging and that the limit, to six decimal places, is 0.607927. The next step, necessarily, takes us beyond C++; we need to recognize this number to formulate (and prove!) a conjecture.

Fortunately, there are good tools for this step. Neil Sloane’s On-Line Encyclopedia of Integer Sequences is a remarkable resource that takes us directly to the answer. Visit http://www.research.att.com/˜njas/sequences/ and enter the sequence of digits into the sequence search engine: 6 0 7 9 2 7 and press the SEARCH button. After a brief delay, the answer emerges:

16

ζ(2) = π2 = 0.6079271018540...

is the number we seek. (The site also gives several references to this well-known problem.)

We close this chapter with a sketch of the proof.

Sweeping all worries about convergence under the rug, consider two large integers. What is the probability they are not both even (a necessary condition for the numbers to be relatively prime)? Each has a 12 chance of being even, so the probability neither has a factor of 2 is 1 − 14 . More generally, the probability neither has a prime p as a common factor is 1 −1/p2 . So the limit of pn is

1

∏p 1 − p2

where the product is over all primes. Recall that ζ (2) is given by

∞ 1

ζ (2) = ∑ .

n=1 n2

This can be expressed as a product. The idea is to factor n2 into even powers of its prime divisors. The product representation is

ζ (2) = ∏p	1 + p2	+ p4	+ p6 + ···
	1	1	1

where the product is over primes p. To see why this works, consider the term (in the sum) 1/1202. We factor 120 as 23 × 3 × 5. Expanding the product representation, the term 1/602 appears by taking 1/26 from the first factor, 1/32 from the second, 1/52 from the third, and 1 from all the other factors.

88	C++ for Mathematicians

Notice that the factors in the product representation are geometric series. Therefore

ζ (2) = ∏p				1			!
		1			1
					2
				− p
and so		1 − p2
ζ (2) = ∏p
1				1

as desired.

5.9Exercises

5.1Calculate ϕ(100), ϕ(29), and ϕ(5!).

5.2What is wrong with this program and how can the mistake be repaired?

#include <iostream> using namespace std; int main() {

int n;

cout << "Enter n: "; cin >> n;

int vals[n];

// do stuff with the vals array return 0;

}

5.3Solve the pair of congruences x ≡ 3 (mod 20) and x ≡ 5 (mod 9).

5.4Write a procedure to solve problems such as Exercise 5.3. It may be declared like this:

long crt(long a1, long n1, long a2, long n2);

(The name crt stands for Chinese Remainder Theorem.)

The procedure should be designed to solve the pair of recurrences

x ≡ a1 (mod n1)

x ≡ a2 (mod n2)

where n1 and n2 are relatively prime positive integers. The return value is the solution mod n1n2.

How should the procedure handle a situation in which n1 and n2 are not of this form?

Arrays

89

5.5In Program 5.7 we defined two values named blank and marked and used them to populate an array named theSieve. For the marks we used char type values (and the array was declared to be char*). However, it would be more logical to use bool values because each cell in theSieve takes only one of two possible values and the Boolean true or false makes perfect sense in this context.

Why did we use char type values instead of bool?

5.6Write a program that fills an array with Fibonacci numbers F0 through F20 and then prints them out in a chart.

5.7Write a program that fills two arrays a and b with integers according to this recurrence:

a0 = b0 = 1 an = bn−1 bn = an−1 + 2bn−1.

The program should then print out a table in which each row is of the form

k	ak	bk	ak	.
			bk

Conjecture a value for limn→∞ an/bn. And, of course, prove your conjecture.

5.8Write a procedure to find the maximum value in an array of long integers. The inputs to the procedure should be the array and the number of elements in the array. The output should be the maximum value in the array.

5.9Write a procedure that generates an array of Fibonacci numbers as its return value. The input to the procedure should be an integer n ≥ 2 that specifies the desired size of the array. Here is how such a procedure would appear in a main():

int main() { long* fibs;

fibs = make_fibs(10);

for (int k=0; k<10; k++) cout << fibs[k] << " "; cout << endl;

return 0;

}

This code should produce the output:

1 1 2 3 5 8 13 21 34 55

In addition, there is a subtle bug in the main(). What is it?

5.10Create a procedure long fibs(int n) to return the nth Fibonacci number. The procedure should work as follows. The first time the procedure is called, it creates a table of Fibonacci numbers holding Fn through F40. Then it returns the value held in the table it created. On all subsequent calls, it does not need

90	C++ for Mathematicians

to recompute any Fibonacci numbers, but simply returns the value in the table it built during its first invocation.

If the input parameter is out of range (either less than 0 or greater than 40) the procedure should return −1.

5.11What happens when new asks for more memory than your computer can provide? Write a program that repeatedly requests new for large blocks of memory without ever releasing those blocks with delete[]. That is, your program should have a severe, deliberate memory leak.

5.12A computational experiment yields the following result: 5.8598744820. Propose a conjecture.

Part II

Objects

Chapter 6

Points in the Plane

The first part of this book introduces the fundamental data types (long, double, bool, etc.), arrays of these types, and procedures.

C++ provides the ability to define new types of data that can be used just as the basic types are. New data types that we create are called classes. In this chapter we create a class called Point that represents a point in the Euclidean plane. When we want a variable to represent an integer quantity, we declare it to be type long (or one of the other integer types). Likewise, once we have the Point class set up, we can declare a variable to represent a point in the plane like this:

Point X;

We call X an object of type Point.

6.1Data and methods

A class definition specifies the data that describe an object of the class and the methods to inspect and manipulate objects.

For the class Point we need to hold data that specify a point’s location in the plane. There are two natural ways we might do this: rectangular coordinates (x,y) or polar coordinates (r,θ). Later in this chapter, when we write the C++ files to create the Point class, we choose the rectangular representation.

A class is more than a way to bundle data together. A class also specifies operations that may be performed on its objects. Here are some things we might want to know about points and actions we might want to perform on points.

•Learn a point’s rectangular coordinates (x and y).

•Learn a point’s polar coordinates (r and θ).

•Change one (or both) of a point’s rectangular coordinates.

•Change one (or both) of a point’s polar coordinates.

•Rotate a point about the origin through a given angle.

•Check if two points are equal.

93

94	C++ for Mathematicians

•Find the distance between two points.

•Find the midpoint between two points.

•Print a point’s coordinates using a statement of the form cout<<P<<endl;.

We perform these various tasks by means of procedures. Many of the procedures are part of the definition of the class itself and are invoked by a different syntax. For example, we create a procedure to change a point’s x coordinate called setX. To set a point P’s x coordinate to −4.5, we use the following special syntax,

P.setX(-4.5);

The setX procedure is part of the definition of the class Point. Computer scientists call such procedures member functions but as we reserve the word function for its mathematical meaning, in this book we call such procedures methods. The term method is also used by many computer scientists.

A C++ class is a bundle that combines data that describe its objects and methods for inspecting and manipulating the objects.

Once a class is defined, we can use it in procedures just as with any other C++ data type. For example, here is a procedure named dist that computes the distance between two points.

double dist(Point P, Point Q) { double dx = P.getX() - Q.getX(); double dy = P.getY() - Q.getY(); return sqrt(dx*dx + dy*dy);

}

Notice that dist is a procedure that works on two arguments (both of type Point) and returns a real number answer of type double. Inside this procedure we use Point’s getX and getY methods to access the x and y coordinates of the two points. (The sqrt procedure is defined in a C++ header file so this code requires a #include <cmath> directive.)

The dist procedure is not a method (i.e., not a member of the Point class); it is simply a procedure that uses the Point data type. It is invoked using the usual syntax; to calculate the distance between two points we call dist(P,Q). However, getX is a method of the Point class; it is used to reveal the point’s x coordinate. Because it is a class method, it is invoked using the special syntax P.getX().

Let’s see how this all works.

6.2Declaring the Point class

When we create a new procedure we break the definition into two files: a header file (whose name ends with .h) that declares the procedure and a code file (whose name ends with .cc) that specifies the algorithm.

Points in the Plane

95

Likewise, a class is defined in two files: a .h file that declares the class and a .cc file that gives the algorithms for its various methods.

The declaration for a class looks like this:

class ClassName {

declarations for data and methods;

};

Notice the required semicolon after the class declaration’s close brace—it is easy to forget.

We now present the declaration for the Point class. There are several new ideas present in this declaration and we examine each. As an aid to readability, we omit all the documentation from the file Point.h. It is extremely important to include comments in the header file that explain what the class represents and what each of its various methods does. The time you take to write these comments will be repaid tenfold when you subsequently write programs that use your classes.

Here is the header file.

Program 6.1: Header file Point.h for the Point class (condensed version).

1 #ifndef POINT_H

2#define POINT_H

3 #include <iostream>

4using namespace std;

5

6class Point {

7

8private:

9 double x;

10 double y;

11

12public:

13Point();

14Point(double xx, double yy);

15double getX() const;

16double getY() const;

17void setX(double xx);

18void setY(double yy);

19double getR() const;

20void setR(double r);

21double getA() const;

22void setA(double theta);

23void rotate(double theta);

24bool operator==(const Point& Q) const;

25bool operator!=(const Point& Q) const;

26

27 };

28

29double dist(Point P, Point Q);

30Point midpoint(Point P, Point Q);

31ostream& operator<<(ostream& os, const Point& P);

32

33 #endif