Dear OCamlers,
In 1994, Barton and Nackman in their book
'Scientific Engineering in
C++' [1] demonstrated how one could encode the
rules of Dimensional
Analysis [2] into the C++ type system enabling
compile-time checking
(no runtime-cost) of the plausibility (at
least up to the dimensional
correctness) of computations.
In 2004, Abrahams & Gurtovy in 'C++
Template Metaprogramming' [3]
showed the Barton Nackman technique to be
elegantly implementable
using compile time type sequences encoding
integer constants. At the
end of this post, I provide a complete listing
of their example
program [4].
The key properties of the system (as I see it)
are:
- Encoding of integers as types;
- Compile time manipulation of sequences of
these integer encodings
to deduce/produce new derived types.
Now, it is not immediately obvious to me how
to approach this problem
in OCaml. It irks me some that I can't
immediately produce a yet more
elegant OCaml program for this problem and
leaves me feeling like C++
has "got something over on us" here ;)
My question therefore is: Does anyone have
suggestions/pointers
on how to approach automatic dimensional
analysis via the OCaml type
system?
Best,
[1] John J. Barton and Lee R. Nackman.
Scientific and Engineering C++:
an Introduction with Advanced Techniques
and Examples. Reading,
MA: Addison Wesley. ISBN 0-201-53393-6.
1994.
[3] David Abrahams and Aleksey Gurtovy C++
Template Metaprogramming:
Concepts, Tools, and Techniques from Boost
and Beyond (C++ in
Depth Series), Addison-Wesley
Professional. ISBN:0321227255. 2004.
[4] Code listing:
//"c:/program files (x86)/Microsoft Visual
Studio 10.0/vc/vcvarsall.bat" x64
//cl /Fedimension.exe /EHsc /I
d:/boost_1_55_0 dimension.cpp
#include <boost/mpl/vector_c.hpp>
#include <boost/mpl/transform.hpp>
#include
<boost/mpl/placeholders.hpp>
#include <boost/mpl/equal.hpp>
#include <boost/mpl/plus.hpp>
#include <boost/mpl/minus.hpp>
#include <boost/static_assert.hpp>
typedef
boost::mpl::vector_c<int,1,0,0,0,0,0,0> mass;
typedef
boost::mpl::vector_c<int,0,1,0,0,0,0,0> length;
typedef
boost::mpl::vector_c<int,0,0,1,0,0,0,0> time;
typedef
boost::mpl::vector_c<int,0,0,0,1,0,0,0> charge;
typedef
boost::mpl::vector_c<int,0,0,0,0,1,0,0> temperature;
typedef
boost::mpl::vector_c<int,0,0,0,0,0,1,0> intensity;
typedef
boost::mpl::vector_c<int,0,0,0,0,0,0,1> angle;
typedef
boost::mpl::vector_c<int,0,1,-1,0,0,0,0> velocity;
// l/t
typedef
boost::mpl::vector_c<int,0,1,-2,0,0,0,0>
acceleration; // l/(t2)
typedef
boost::mpl::vector_c<int,1,1,-1,0,0,0,0> momentum;
// ml/t
typedef
boost::mpl::vector_c<int,1,1,-2,0,0,0,0> force;
// ml/(t2)
typedef
boost::mpl::vector_c<int,0,0,0,0,0,0,0> scalar;
template <class T, class Dimensions>
class quantity
{
public:
explicit quantity (T val)
: val (val)
{}
template <class OtherDimensions>
quantity (quantity<T,
OtherDimensions> const& other)
: val (other.value ()) {
BOOST_MPL_ASSERT(
(boost::mpl::equal<Dimensions, OtherDimensions>));
}
T value () const { return val; }
private:
T val;
};
template <class T, class D>
quantity<T, D>
operator + (quantity<T, D> x,
quantity<T, D> y )
{
return quantity<T, D>(x.value () +
y.value ());
}
template <class T, class D>
quantity<T, D>
operator - (quantity<T, D> x,
quantity<T, D> y )
{
return quantity<T, D>(x.value () -
y.value ());
}
template <class T, class D1, class
D2>
quantity <
T
, typename boost::mpl::transform<
D1, D2, boost::mpl::plus<
boost::mpl::placeholders::_1
,
boost::mpl::placeholders::_2> >::type
>
operator* (quantity<T, D1> x,
quantity <T, D2> y)
{
typedef typename
boost::mpl::transform<
D1, D2, boost::mpl::plus<
boost::mpl::placeholders::_1
,
boost::mpl::placeholders::_2> >::type D;
return quantity<T, D> (x.value ()
* y.value ());
}
template <class T, class D1, class
D2>
quantity <
T
, typename boost::mpl::transform<
D1, D2, boost::mpl::minus<
boost::mpl::placeholders::_1
,
boost::mpl::placeholders::_2> >::type
>
operator/ (quantity<T, D1> x,
quantity <T, D2> y)
{
typedef typename
boost::mpl::transform<
D1, D2, boost::mpl::minus<
boost::mpl::placeholders::_1
,
boost::mpl::placeholders::_2> >::type D;
return quantity<T, D> (x.value ()
/ y.value ());
}
// -- test
#include <iostream>
#include <limits>
#include <cassert>
int main ()
{
quantity<float, mass> m (5.0f);
quantity<float, acceleration>
a(9.8f);
quantity<float, force> f = m * a;
quantity<float, mass> m2 = f / a;
assert ((std::abs ((m2 - m).value ()))
<= std::numeric_limits<double>::epsilon ());
return 0;
}