Wednesday, January 24, 2007

units and dimensions with C++ templates

I personally consider templates to be the most powerful and distinguishing feature of C++ across all languages (note that I surely don't know all programming languages).

In this post, I shall show you the control they provide in handling units and dimensions in programming in a very elegant way. I got the motivation to read and think about it from Scott Meyers Most Important Aha moments - Scott Meyers Most Important Aha Moments (Understanding the use of non-type template parameters in Barton’s and Nackman’s approach to dimensional analysis).

I wrote a piece of sample code to illustrate it on codeguru discussion forum. The approach gives compile time validation in various computations so that they are dimensionally correct and various units are handled in a less error prone manner.

I shall not blabber much about it as you may find a lot to read about it and understand it from the web and will drop off the sample code for you to take a look at. It is a very initial attempt and I seriously hope that you can better it with adding multiple dimensions apart from just 3 that I used (mass, length and time).

Here is a working piece of code illustrating the approach.

[CODE]
#include<iostream>

//internally for everything use SI units.
//mass in kilograms,
//time in seconds
//length in meters

template <int M, int L, int T>
class Unit;

template <int M, int L, int T>
const Unit<M,L,T> operator+(const Unit<M,L,T>&, const Unit<M,L,T>&);

template <int M, int L, int T>
const Unit<M,L,T> operator-(const Unit<M,L,T>&, const Unit<M,L,T>&);

template<int M, int L, int T>
class Unit{
               public:
                               Unit(double val=0.0) : value(val){}
                               Unit& operator+=(const Unit& rhs){
                                               value+=rhs.value;
                                               return *this;
                               }
                               Unit& operator-=(const Unit& rhs){
                                               value-=rhs.value;
                                               return *this;
                               }
                               double getValue() const{
                                               return value;
                               }
                               private:
                                               double value;
                               friend const Unit<M,L,T> operator+<M, L, T>(const Unit<M,L,T>&, const Unit<M,L,T>&);
                               friend const Unit<M,L,T> operator-<M, L, T>(const Unit<M,L,T>&, const Unit<M,L,T>&);

};

template <int M, int L, int T>
const Unit<M,L,T> operator+(const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return Unit<M,L,T>(lhs)+=rhs;
}
template <int M, int L, int T>
const Unit<M,L,T> operator-(const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return Unit<M,L,T>(lhs)-=rhs;
}
template <int M1, int L1, int T1, int M2, int L2, int T2>
const Unit<M1+M2,L1+L2,T1+T2> operator*(const Unit<M1,L1,T1>& lhs, const Unit<M2,L2,T2>& rhs)
{
               return Unit<M1+M2,L1+L2,T1+T2>(lhs.getValue()*rhs.getValue());
}
template <int M, int L, int T>
const Unit<M,L,T> operator*(const double& lhs, const Unit<M,L,T>& rhs)
{
               return Unit<M,L,T>(lhs*rhs.getValue());
}

template <int M1, int L1, int T1, int M2, int L2, int T2>
const Unit<M1-M2,L1-L2,T1-T2> operator/(const Unit<M1,L1,T1>& lhs, const Unit<M2,L2,T2>& rhs)
{
               return Unit<M1-M2,L1-L2,T1-T2>(lhs.getValue()/rhs.getValue());
}

template <int M, int L, int T>
bool operator==(const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return (lhs.getValue()==rhs.getValue());
}
template <int M, int L, int T>
bool operator<=(const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return lhs.getValue()<=rhs.getValue();
}
template <int M, int L, int T>
bool operator>=(const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return lhs.getValue()>=rhs.getValue();
}
template <int M, int L, int T>
bool operator< (const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return lhs.getValue()<rhs.getValue();
}
template <int M, int L, int T>
bool operator> (const Unit<M,L,T>& lhs, const Unit<M,L,T>& rhs)
{
               return lhs.getValue()>rhs.getValue();
}

typedef Unit<1,0,0> Mass;
typedef Unit<0,1,0> Length;
typedef Unit<0,0,1> Time;

//length units - verify the constants
const Length meter = 1.0;
const Length centimeter = 0.01;
const Length feet = 0.3048*meter;
const Length mile = 5280*feet;

//time units
const Time second = 1.0;

//mass units
const Mass kilogram = 1.0;

//Can have more constants for all other units in terms of meter, second and kilogram
//feet,yard,inches,micrometer,hour,minute,grams,pounds, etc etc etc....

int main()
{
               Length distanceinCMs = 100000;
               Length distanceinFeets = distanceinCMs*centimeter/feet;
               std::cout << "distanceinCMs = " << distanceinCMs.getValue() << std::endl;
               std::cout << "distanceinFeets = " << distanceinFeets.getValue() << std::endl;
               //get in terms of miles and feets
               double total_in_miles = (distanceinCMs*centimeter/mile).getValue();
               double total_in_feets = (distanceinCMs*centimeter/feet).getValue();

               std::cout << "above distance = " << static_cast<int>(total_in_miles)
               << " miles and " << static_cast<int>(((total_in_miles -
                                               static_cast<int>(total_in_miles))*mile/feet).getValue())
                                               << " feet (approx.)" << std::endl;
               return 0;
}


If you run this code, with the sample main that I provided, you will see an output as follows:

[OUTPUT]
distanceinCMs = 100000
distanceinFeets = 3280.84
above distance = 0 miles and 3280 feet (approx.)


Isn't that very elegant and easy to use? The best part is, it takes care of everything in SI units for the calculations but you can provide constants for any unit (as defined above, for example 'feet') and you would just need to specify it with the values/variables (as shown in the sample above).

4 comments:

Anonymous said...

Your units example is simple and elegant. I'm writing a general-purpose chunk of code to make calculations in C++ type safe. If you'd like to review/comment/contribute, please contact me via donp at gdssw.com.

Anonymous said...

The idea seems to have been quite popular some years ago. I had a solution similar to yours. but you sholud read this paper
http://www.oonumerics.org/tmpw01/brown.pdf
it gives a solution with rational exponents.

Rushikesh Deshmukh said...

It would be nice to restrict operations of +, -, *, /, = etc. for dimensions having same units. It might result in detecting errors at compile time itself. Is this possible with above design?

Unknown said...

Rushikesh,
+ and - are restricted to objects with same 'dimensions' (template arguments). The point of the library is that your pick an internal representation that works for everything (for example meters internally for all lengths).

The real problem comes with temperature... C to F is not linear. Even Kelvin to C is not linear.