C++ For Mathematicians (2006) [eng]

186	C++ for Mathematicians

10.5 Class and file organization for PPoint and PLine

Because of the duality between points and lines in the projective plane, most of the data and code we use to represent these concepts in C++ are the same. If at all possible, we should avoid writing the same code twice for two reasons. First, the initial work in creating the programs is doubled. More important, maintaining the code is also made more difficult; if a change is required to the code, we need to remember to make that change twice. When we are fussing with our programs and making a number of minor modifications, it is easy to forget to update both versions. To illustrate this, we deliberately include a subtle flaw in the first version of the program; we then repair the problem once (not twice).

To this end, we define three classes: a parent class named PObject that contains as much of the code as possible and two derived classes, PPoint and PLine, that include code particular to each. These classes are defined in two files each: a header

.h file and a code .cc file. Figure 10.3 illustrates the organization.

			PObject
			PObject.h
			PObject.cc
PPoint						PLine
PPoint.h						PLine.h
PPoint.cc						PLine.cc

Figure 10.3: Hierarchy of the PObject classes.

The header file PObject.h is used to declare the PObject class. Both PPoint.h and PLine.h require the directive #include "PObject.h". Programs that use projective geometry need all three. Rather than expecting the user to type multiple #include directives, we create a convenient header file that includes everything we need. We call this header file Projective.h; here it is.

Program 10.3: Header file for all projective geometry classes, Projective.h.

1#ifndef PROJECTIVE_H

2#define PROJECTIVE_H

3#include "PPoint.h"

4#include "PLine.h"

5#endif

The Projective Plane

187

Notice that we do not require #include "PObject.h" because that file is already #included by PPoint.h and PLine.h. All these files have the usual structure to prevent multiple inclusion.

All told, we have seven files that implement points and lines in the projective plane.

PObject.h The header file for the PObject class.

PObject.cc The C++ code for the PObject class.

PPoint.h The header file for the PPoint class.

PPoint.cc The C++ code for the PPoint class.

PLine.h The header file for the PLine class.

PLine.cc The C++ code for the PLine class.

Projective.h The only header file subsequent programs need to include to use projective points and lines.

The decision to write the program in seven different files is based on the principle of breaking a problem down to manageable sizes and working on each piece separately. We could have packed all this work into two larger files (one .h and one

.cc).

10.6 The parent class PObject

Data and functionality common to PPoint and PLine are implemented in the class PObject. Using the ideas presented in Section 10.2, we map out the class

PObject.

Data A point or line in RP2 is represented by homogeneous coordinates: (x,y,z) or [x,y,z]. To hold these coordinates, we declare three private double variables, x, y, and z. Let us write x,y,z for the homogeneous coordinates of a generic projective object (either a point of a line). (See Program 10.4, line 9.)

Because 2,1,−5 and −4,−2,10 name the same object, it is useful to choose a canonical triple. One idea is to make sure that x,y,z is a unit vector (but then we have a sign ambiguity). The representation we use is to make the last nonzero coordinate of the triple equal to 1. The motivation for this choice is that a point in the Euclidean plane at coordinates (x,y) corresponds to the point (x,y,1) in the projective plane.

If later we are unhappy with this decision, we can choose another manner to store the homogeneous coordinates. Because the coordinates are private members of PObject, we would only need to update the code for PObject.

188	C++ for Mathematicians

We provide public methods getX(), getY(), and getZ() to inspect (but not modify) the values held in x, y, and z, respectively. (See Program 10.4 lines 33–35.)

We reserve 0,0,0 to stand for an invalid projective object. We provide a public method is_invalid() to check if an object is invalid. (See Program 10.4 lines 37–39.)

Constructors It is natural to define a constructor for a PObject with three parameters that set the homogeneous coordinates of the object. The user might invoke the constructor like this:

PObject P(2., 3., 5.);

Rather than holding this point as 2,3,5 , we use the canonical representation0.4,0.6,1 . (See Program 10.4 lines 25–30.)

To facilitate the conversion of user-supplied coordinates to canonical coordinates, we create a private helper procedure scale. (See Program 10.4 line 10 and Program 10.5 lines 4–19.)

All C++ classes ought to define a zero-argument constructor that creates a default object of that class. In this case, a sensible choice is the object 0,0,1 . This corresponds to the origin (as a point) or the line at infinity (as a line). (See Program 10.4 lines 21–24.)

Relations We want to be able to compare projective points (and lines) for equality, inequality, and < (for sorting). To this end, we might be tempted to define operator== in the public portion of PObject like this:

public:

bool operator==(const PObject& that) const {

return ( (x==that.x) && (y==that.y) && (z==that.z) );

}

Although this would give the desired result when comparing points to points or lines to lines, it would also provide the unfortunate ability to compare points to lines. Were we to use the above method for testing equality, then we might run into the following situation.

PPoint P(2., 3., 5.);

PLine L(2., 3., 5.); if (P == L) {

cout << "They are equal." << endl;

}

When this code runs, the message They are equal would be written on the computer’s screen.

There are at least two problems with this approach. First, the point (2,3,5) and the line [2,3,5] are not equal even though they have the same homogeneous

The Projective Plane

189

coordinates. Second, lines and points should not even be comparable by ==. Code that compares a point to a line is almost certainly a bug; the best thing in this situation is for the compiler to flag such an expression as an error.

We therefore take a different approach to equality testing. We want the fundamental code that checks for equality to reside inside the PObject class because that code is common to both PPoint and PLine. We do this by declaring a protected method called equals inside PObject (see Program 10.4 line 14):

protected:

bool equals(const PObject& that) const;

In the file PObject.cc we give the code for this method (see Program 10.5 lines 34–36):

bool PObject::equals (const PObject& that) const { return ( (x==that.x) && (y==that.y) && (z==that.z) );

}

Then, in the class definitions for PPoint and PLine we give the necessary operator definitions. For example, in PLine we have this:

public:

bool operator==(const PLine& that) const { return equals(that);

}

bool operator!=(const PLine& that) const { return !equals(that);

}

In this way, if (indeed, when) we decide to change the manner in which we test for equality, we only need to update PObject’s protected equals method.

Similarly, rather than defining a < operator for PObject, we define a protected less method. The children can access less to define their individual, public operator< methods. (See Program 10.4 line 15 and Program 10.5 lines 38– 45.)

Meet/Join Operation Given two PPoints P and Q, we want P+Q to return the line through those points. Dually, given PLines L and M, we want L*M to return the point of intersection of these lines.

In both cases, the calculations are the same: given the triples x1,y1,z1 andx2,y2,z2 we need to find a new triple x3,y3,z3 that is orthogonal to the first two. To do this, we calculate the cross product:

x3,y3,z3 = x1,y1,z1 × x2,y2,z2 .

Therefore, in PObject we declare a protected method called op that is invoked by operator+ in PPoint and operator* in PLine. (See Program 10.4 line 18 and Program 10.5 lines 138–148.)

190	C++ for Mathematicians

Incidence Given a point and a line, we want to be able to determine if the point lies on the line. If the coordinates for these are (x,y,z) and [a,b,c], respectively, then we simply need to test if ax + by + cz = 0.

It does not make sense to ask if one line is incident with another, so we do not make an “is incident with” method publicly available in PObject. Rather, we make a protected incident method (that calls, in turn, a private dot method for calculating dot product). (See Program 10.4 lines 11,16 and Program 10.5 lines 21–23,47–49.)

Then, PPoint and PLine can declare their own public methods that access incident. Details on this later.

Collinearity/Concurrence Are three given points collinear? Are three given lines concurrent? If the three objects have coordinates x1,y1,z1 , x2,y2,z2 , andx3,y3,z3 , then the answer is yes if and only if the vectors are linearly dependent. We check this by calculating

	x1	y1	z1
det	x2	y2	z2
	x3	y3	z3

and seeing if we get zero. In PObject.h we declare a procedure dependent that takes three PObject arguments and returns true or false. We then use dependent to implement procedures collinear and concurrent for the classes PPoint and PLine (respectively). (See Program 10.4 line 44 and Program 10.5 lines 60–77.)

Random points/lines There is no way to generate a point uniformly at random in the Euclidean plane, but there is a sensible way in which we can do this for the projective plane. Recall that points in RP2 correspond to lines through the origin in R3. Thus, to select a point at random in RP2 we generate a point uniformly at random on the unit ball centered at the origin. An efficient way to perform this latter task is to select a vector v uniformly in [−1,1]3. If v > 1, then we reject v and try again.

The method for generating a random line is precisely the same.

We therefore include a public method randomize that resets the coordinates of the PObject by the algorithm we just described. (See Program 10.4 line 31 and Program 10.5 lines 25–32.)

To choose a random line through a point is similar. Suppose the point is (x,y,z). The line [a,b,c] should be chosen so that [a,b,c] is orthogonal to (x,y,z). To do this, we find an orthonormal basis for (x,y,z) that we denote {(a1,b1,c1),(a2,b2,c2)}. We then choose t uniformly at random in [0,2π]. The random line [a,b,c] is given by

a = a1 cost + a2 sint b = b1 cost + b2 sint c = c1 cost + c2 sint.

The Projective Plane

191

Rather than generate t and then compute two trigonometric functions, we can obtain the pair (cost,sint) by choosing a point uniformly at random in the unit disk (in R2) and then scaling.

The technique for selecting a random point on a given line is exactly the same.

Thus, we define a protected rand_perp method in PObject (Program 10.4 line 17 and Program 10.5 lines 79–136).

The rand_perp method is used by rand_point in PLine and rand_line in PPoint.

Output Finally, PObject.h declares a procedure for writing to an output stream. The format is <x,y,z>. The ostream operator<< defined in PObject is overridden by like-named procedures in PPoint and PLine. (See Program 10.4 line 42 and Program 10.5 lines 51–58.)

With these explanations in place, we now give the header and code files for the class PObject.

Program 10.4: Header file for the PObject class (version 1).

1#ifndef POBJECT_H

2#define POBJECT_H

3

4#include <iostream>

5using namespace std;

6

7class PObject {

8private:

9double x,y,z;

10void scale();

11double dot(const PObject& that) const;

12

13protected:

14bool equals(const PObject& that) const;

15bool less(const PObject& that) const;

16bool incident(const PObject& that) const;

17PObject rand_perp() const;

18PObject op(const PObject& that) const;

19

20public:

21PObject() {

22x = y = 0.;

23z = 1.;

24}

25PObject(double a, double b, double c) {

26x = a;

27y = b;

28z = c;

29scale();

30}

31void randomize();

32

33 double getX() const { return x; }

192	C++ for Mathematicians

34double getY() const { return y; }

35double getZ() const { return z; }

36

37bool is_invalid() const {

38return (x==0.) && (y==0.) && (z==0.);

39}

40};

41

42 ostream& operator<<(ostream& os, const PObject& A);

43

44 bool dependent(const PObject& A, const PObject& B, const PObject& C);

45

46 #endif

Program 10.5: Program file for the PObject class (version 1).

1#include "PObject.h"

2#include "uniform.h"

3

4void PObject::scale() {

5if (z != 0.) {

6x /= z;

7y /= z;

8z = 1.;

9return;

10}

11if (y != 0.) {

12x /= y;

13y = 1.;

14return;

15}

16if (x != 0) {

17x = 1.;

18}

19}

20

21double PObject::dot(const PObject& that) const {

22return x*that.x + y*that.y + z*that.z;

23}

24

25void PObject::randomize() {

26do {

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

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

29z = unif(-1.,1.);

30} while (x*x + y*y + z*z > 1.);

31scale();

32}

33

34bool PObject::equals(const PObject& that) const {

35return ( (x==that.x) && (y==that.y) && (z==that.z) );

36}

37

38bool PObject::less(const PObject& that) const {

39if (x < that.x) return true;

The Projective Plane

193

40if (x > that.x) return false;

41if (y < that.y) return true;

42if (y > that.y) return false;

43if (z < that.z) return true;

44return false;

45}

46

47bool PObject::incident(const PObject& that) const {

48return dot(that)==0.;

49}

50

51ostream& operator<<(ostream& os, const PObject& A) {

52os << "<"

53<< A.getX() << ","

54<< A.getY() << ","

55<< A.getZ()

56<< ">";

57return os;

58}

59

60bool dependent(const PObject& A, const PObject& B, const PObject& C){

61double a1 = A.getX();

62double a2 = A.getY();

63double a3 = A.getZ();

64

65double b1 = B.getX();

66double b2 = B.getY();

67double b3 = B.getZ();

68

69double c1 = C.getX();

70double c2 = C.getY();

71double c3 = C.getZ();

72
73	double det =	a1*b2*c3	+	a2*b3*c1	+	a3*b1*c2
74	-	a3*b2*c1	-	a1*b3*c2	-	a2*b1*c3;
75

76return det == 0.;

77}

78

79PObject PObject::rand_perp() const {

80if (is_invalid()) return PObject(0,0,0);

81
82	double x1,y1,z1;	//	One vector orthogonal to (x,y,z)
83	double x2,y2,z2;	//	Another orthogonal to (x,y,z) and (x1,y1,z1)
84

85if (z == 0.) { // If z==0, take (0,0,1) for (x1,y1,y2)

86x1 = 0;

87y1 = 0;

88z1 = 1;

89}

90else {

91if (y == 0.) { // z != 0 and y == 0, use (0,1,0)

92x1 = 0;

93y1 = 1;

94z1 = 1;

95}

194	C++ for Mathematicians

96else { // y and z both nonzero, use (0,-z,y)

97x1 = 0;

98y1 = -z;

99z1 = y;

100}

101}

102

103// normalize (x1,y1,z1)

104double r1 = sqrt(x1*x1 + y1*y1 + z1*z1);

105x1 /= r1;

106y1 /= r1;

107z1 /= r1;

108

109// (get x2,y2,z2) by cross product with (x,y,z) and (x1,y1,z1)

110x2 = -(y1*z) + y*z1;

111y2 = x1*z - x*z1;

112z2 = -(x1*y) + x*y1;

113

114// normalize (x2,y2,z2)

115double r2 = sqrt(x2*x2 + y2*y2 + z2*z2);

116x2 /= r2;

117y2 /= r2;

118z2 /= r2;

119

120// get a point uniformly on the unit circle

121double a,b,r;

122do {

123a = unif(-1.,1.);

124b = unif(-1.,1.);

125r = a*a + b*b;

126} while (r > 1.);

127r = sqrt(r);

128a /= r;

129b /= r;

130

131double xx = x1 * a + x2 * b;

132double yy = y1 * a + y2 * b;

133double zz = z1 * a + z2 * b;

134

135return PObject(xx,yy,zz);

136}

137

138

139 PObject PObject::op(const PObject& that) const {

140

141 if (equals(that)) return PObject(0,0,0);

142

143double c1 = y*that.z - z*that.y;

144double c2 = z*that.x - x*that.z;

145double c3 = x*that.y - y*that.x;

146

147return PObject(c1,c2,c3);

148}

The Projective Plane

195

10.7 The classes PPoint and PLine

With the code for PObject in place, we are ready to finish our work by writing the files for the classes PPoint and PLine. We do this work in four files: PPoint.h, PPoint.cc, PLine.h, and PLine.cc (Programs 10.6 through 10.9).

As we start to write these files, we meet a chicken-and-egg problem. Which do we define first: the class PPoint or the class PLine? For us, this is more than a philosophical conundrum. Some PLine methods need to refer to PPoints. For example, operator* acts on lines to produce points and PLine’s rand_point method returns a PPoint. Dually, some of the PPoint methods require the PLine class.

Here is how we solve this dilemma. In the PPoint.h file, before we give the definition of PPoint we have the following statement (see line 6 of Program 10.6),

class PLine;

This lets the C++ compiler know that a class named PLine is defined somewhere else. Therefore, when the compiler sees PLine, it knows that this identifier refers to some class. However, it does not know anything else about PLine other than the fact it is a class. This means that we may not give an inline definition for any method that refers to objects of type PLine. Later, in the file PPoint.cc, we include all the headers (conveniently with the single directive #include "Projective.h"). Therefore, the code in PPoint.cc can make full use of the PLine class.

Dually, we include the statement class PPoint; in the file PLine.h.

We focus our attention on the class PPoint. The analysis of PLine is similar. For the class PPoint we have five constructors. At first, this may seem to be too

many, but we show how to do this easily.

Of course, we want a three-argument constructor PPoint(a,b,c) that creates the point (a,b,c). We also want a zero-argument constructor PPoint() that creates the point (0,0,1) corresponding to the origin.

It makes sense to define a two-argument constructor PPoint(a,b) to create the point (a,b,1). In a sense, this maps the Euclidean point (a,b) to its natural correspondent in the RP2.

What should be the action of a single-argument constructor? The real number a can be identified with the point (a,0) on the x-axis, and this in turn corresponds to (a,0,1) in RP2. In summary, we have the following constructors and their effects.

Constructor form	Point created
PPoint P();	(0,0,1)
PPoint P(a);	(a,0,1)
PPoint P(a,b);	(a,b,1)
PPoint P(a,b,c);	(a,b,c)

The great news is that we can implement these four constructors with a single definition using default parameters.