Tuesday, 7 October 2008

Contraction of Four-Vectors: A C++ Exercise

NOTE: The examples as typed have not been compiled or checked, however they are directly based on examples I talked through while teaching C++, so there should not be any particularly glaring errors or omissions. Note that I've avoided talking too much about object-oriented design, advanced concepts or why certain things are done the way they are. This article is intended mainly as an introduction to the abstraction mechanisms available with C++ class-based design, and is not intended to introduce any new concepts to the world of C++. In particular, I do not claim to be the first person to have thought of such concepts (as a number of comments have suggested in the past); I'm merely providing an example for those interested!

Without going into the physics behind it, if we have two objects A and B, each with four components, we can define their contraction as:
A_{\mu}B^{\mu} = A_0 B_0 - (A_1 B_1 + A_2 B_2 + A_3 B_3)


I intend now to use this concept to develop a C++ program from a `first draft' to a finished object-oriented program. This is indented primarily as an introductory-level tutorial for those new to C++ but familiar with many aspects of C or similar languages, so I will not go into a lot of detail about the basics, instead focussing on the C++ specific stuff.

The first approach could look something like this:

#include <iostream>

int main()
{
// Create some variables
double a0, a1, a2, a3;
double b0, b1, b2, b3;

// Prompt and read eight components
std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
std::cin >> a0 >> a1 >> a2 >> a3;
std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
std::cin >> b0 >> b1 >> b2 >> b3;

// Calculate contraction and display result
double result = a0*b0 - a1*b1 - a2*b2 - a3*b3;
std::cout << "A.B = " << result << std::endl;

return 0;
}


This is a fairly simple program, achieving exactly what we set out to do. However, from the perspective of simplicity, there's a lot of overhead in the main program which detracts from the understanding of what is actually going on. We can make things slightly clearer by writing the contraction as a function; then we'll be able to call it repeatedly from our code, instead of defining what a contraction of four-vectors means at every point it is used:

double contract(double a0, double a1, double a2, double a3,
double b0, double b1, double b2, double b3)
{
return a0*b0 - a1*b1 - a2*b2 - a3*b3;
}


Now, our main() looks like this:

int main()
{
// Create some variables
double a0, a1, a2, a3;
double b0, b1, b2, b3;

// Prompt and read eight components
std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
std::cin >> a0 >> a1 >> a2 >> a3;
std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
std::cin >> b0 >> b1 >> b2 >> b3;

// Calculate contraction and display result
double result = contract(a0, a1, a2, a3, b0, b1, b2, b3);
std::cout << "A.B = " << result << std::endl;

return 0;
}

No major improvement there, but we've moved the contraction code out into a function so that we can reuse it. Still, it's pretty ugly to pass eight doubles into the function. Much nicer would be to pass two four-vectors; after all, that is what the contraction works with, at a higher level:

#include <iostream>

class FourVector {
public:
double x0;
double x1;
double x2;
double x3;
};

double contract(FourVector a, FourVector b)
{
return a.x0*b.x0 - a.x1*b.x1 - a.x2*b.x2 - a.x3*b.x3;
}

int main()
{
// Create some variables
Fourvector a, b;

// Prompt and read eight components
std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
std::cin >> a.x0 >> a.x1 >> a.x2 >> a.x3;
std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
std::cin >> b.x0 >> b.x1 >> b.x2 >> b.x3;

// Calculate contraction and display result
double result = contract(a, b);
std::cout << "A.B = " << result << std::endl;

return 0;
}


Things are a little neater now, but still not ideal. For a start, we have public member variables in our class. This is generally considered bad design, so we should make them private and instead provide alternative mechanisms for accessing them. In this case, the class lends itself perfectly towards the use of operator overloading. The idea here is to define what it means to "multiply" two four-vectors together. We can define operator*() as a member-function of the FourVector class:

#include <iostream>

class FourVector {
public:
double x0;
double x1;
double x2;
double x3;
double operator*(const FourVector& b);
};

double FourVector::operator*(const FourVector& b)
{
return x0*b.x0 - x1*b.x1 - x2*b.x2 - x3*b.x3;
}

int main()
{
// Create some variables
Fourvector a, b;

// Prompt and read eight components
std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
std::cin >> a.x0 >> a.x1 >> a.x2 >> a.x3;
std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
std::cin >> b.x0 >> b.x1 >> b.x2 >> b.x3;

// Calculate contraction and display result
double result = a * b;
std::cout << "A.B = " << result << std::endl;

return 0;
}


As you can see, we have now defined enough code that the statement result = a * b; works as expected, producing the contraction of the two. The member function, when called, is bound to object a (in this case), so x0, x1, x2, x3 all refer to a's variables. We have to refer to the variables in b (passed as an argument to the operator) explicitly as b.x0 etc.

We still haven't hidden the member variables. To do this, we need to add a couple more things: a constructor, and a stream operator.

#include <iostream>

class FourVector {
public:
FourVector() : x0(0), x1(0), x2(0), x3(0) {};
FourVector(double i0, double i1, double i2, double i3)
: x0(i0), x1(i1), x2(i2), x3(i3) {};
double operator*(const FourVector& b);
friend std::istream& operator>>(std::istream& in, FourVector& v);
private:
double x0;
double x1;
double x2;
double x3;
};

double FourVector::operator*(const FourVector& b)
{
return x0*b.x0 - x1*b.x1 - x2*b.x2 - x3*b.x3;
}

std::istream& operator>>(std::istream& in, FourVector& v)
{
in >> v.x0 >> v.x1 >> v.x2 >> v.x3;
return in;
}

int main()
{
// Create some variables
Fourvector a, b;

// Prompt and read eight components
std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
std::cin >> a;
std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
std::cin >> b;

// Calculate contraction and display result
double result = a * b;
std::cout << "A.B = " << result << std::endl;

return 0;
}

There are a few things to note here. First, the default constructor:
FourVector() : x0(0), x1(0), x2(0), x3(0) {};


Constructors are named identically to the class they construct, and have no return type. This particular one takes no arguments, has an initialiser list which sets all the member variables to zero, and has an empty function body -- {};

The other constructor allows one to construct a FourVector object by passing four doubles into the constructor:

double E, px, py, pz;
// Some code which sets E and p goes here
FourVector EnergyMomentum(E, px, py, pz);


Finally, operator>> is a little different. Because the order is stream >> FourVector, we can't make the operator>> a member of FourVector (since it would expect the FourVector object to be passed first!). Furthermore, we need access to the private member variables, or we have nothing to read into. This is a bit of a dilemma, until we use the `friend' statement. This tells the compiler to let a function which is defined outside the class have access to internal class members; precisely what we want here!

You'll notice that the mechanism by which we got here is quite convoluted; we've added constructors and defined operators, and the class is far from complete -- for it to be generally useful we'd need an output operator, other mathematical operations and more. The amount of code has increased, as has the complexity of the overall solution. But look at the code inside main()... this is now considerably simpler. If we ignore comments and don't prompt the user, we can now write exactly the code which corresponds to the mathematical concept:

int main()
{
FourVector a, b;
std::cin >> a >> b;
double result = a * b;
std::cout >> result >> std::endl;
return 0;
}


The best part is that the FourVector class is reusable, and sufficiently powerful to allow this abstraction to be carried forward to develop a whole system for tensor calculus, expressed in terms of the mathematical objects involved!