C++中编写物理公式

物理计算中，数值都不是独立的，大多数物理量都附带有量纲。随着计算变得越来越复杂，维持物理量的正确量纲能够避免诸如“将质量赋给长度”之类不经意间犯下的错误。我们接下来就要为数值建立一个类型系统。

国际标准量纲制规定了物理量的基本量纲为：质量，长度，时间，电荷，温度，密度以及物质的量。其他量纲都是在此基础上复合而成的。我们将要实现的类型系统必须能够表示这七种以及其组合的量纲。

一个复合量纲可以看成是几个基本量纲的乘积。我们打算用一个数组表示他们，其值表示对应量纲的幂次： 

typedef int dimension[7]; // m l t ...
dimension const mass    = {1, 0, 0, 0, 0, 0, 0};
dimension const length  = {0, 1, 0, 0, 0, 0, 0};
dimension const time    = {0, 0, 1, 0, 0, 0, 0};
...

根据这种表示法，力的量纲表示如下:

dimension const force = {1, 1, -2, 0, 0, 0, 0};

即\(mlt^{-2}\)。然而，这些数组的使用并不能达到我们的目的：它们的类型都是相同的，而质量和长度的类型不同。

使用MPL中的vector可以创建一组不同类型的数组。 

#include <boost/mpl/vector.hpp>

typedef boost::mpl::vector<signed char, short, int, long> signed_types;

使用整形外覆器可以用类型表示数量。 MPL提供了int_<N>类模板，它以一个內嵌的::value来表示他的整形参数N： 

#include <boost/mpl/int.hpp>

namespace mpl = boost::mpl;
static int const five = mpl::int_<5>::value;

MPL还提供了不同的vector让我们写出简短的代码: 

#include <boost/mpl/vector_c.hpp>

typedef mpl::vector_c<int, 1, 0, 0, 0, 0, 0, 0> mass;
typedef mpl::vector_c<int, 0, 1, 0, 0, 0, 0, 0> length;
typedef mpl::vector_c<int, 0, 0, 1, 0, 0, 0, 0> time;
typedef mpl::vector_c<int, 0, 0, 0, 1, 0, 0, 0> charge;
typedef mpl::vector_c<int, 0, 0, 0, 0, 1, 0, 0> temperature;
typedef mpl::vector_c<int, 0, 0, 0, 0, 0, 1, 0> intensity;
typedef mpl::vector_c<int, 0, 0, 0, 0, 0, 0, 1> amount_of_substance;

如果我们愿意，还可以定义一些复合量纲： 

typedef mpl::vector_c<int, 0, 1, -1, 0, 0, 0, 0> velocity;       // l/t
typedef mpl::vector_c<int, 0, 1, -2, 0, 0, 0, 0> acceleration;   // l/(t^2)
typedef mpl::vector_c<int, 1, 1, -1, 0, 0, 0, 0> momentum;       // ml/t
typedef mpl::vector_c<int, 1, 1, -2, 0, 0, 0, 0> force;          // ml/(t^-2)
typedef mpl::vector_c<int, 0, 0, 0, 0, 0, 0, 0> scalar;

我们要进行真实的计算，还要进行类型检查，这需要一个简单的外覆类： 

template<class T, class Dimension>
struct quantity {
    explicit quantity(T x): m_value(x) {}
    T value() const {return m_value;}
    private:
    T m_value;
};

现在，我们把数值和量纲联系起来了： 

quantity<float, length> l(1.0f);
quantity<float, mass> m(2.0f);

m = l; // compile error

仅仅这样我们还不能对数据进行计算。下面我们先来实现加法和减法： 

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());
}

这样就可以正确地实现两个物理量的加法了: 

quantity<float, length> len1(1.0f);
quantity<float, length> len2(2.0f);

len1 = len1 + len2; // ok
len1 = len2 + quantity<float, mass>(3.7);// compile error

接下来，为了实现乘法和除法，我们需要对量纲进行计算，获得计算结果的量纲。plus和minus两个boost提供的模板类可以满足要求： 

#include <boost/mpl/plus.hpp>
#include <boost/mpl/int.hpp>
#include <iostream>

typedef mpl::plus< mpl::int_<2>, mpl::int_<3> >::type pls2_3;
int main() {
    std::cout << pls2_3::value << std::endl;
    return 0;
}

到这里应该就可以了。我想读者可以很容易猜到transform的用途。 

#include <boost/mpl/int.hpp>
#include <boost/mpl/vector_c.hpp>
#include <boost/mpl/plus.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/equal.hpp>
#include <boost/mpl/transform.hpp>
#include <boost/static_assert.hpp>
#include <iostream>

namespace mpl = boost::mpl;

typedef mpl::vector_c<int, 1, 0, 0, 0, 0, 0, 0> mass;
typedef mpl::vector_c<int, 0, 1, 0, 0, 0, 0, 0> length;
typedef mpl::vector_c<int, 0, 0, 1, 0, 0, 0, 0> time_;
typedef mpl::vector_c<int, 0, 0, 0, 1, 0, 0, 0> charge;
typedef mpl::vector_c<int, 0, 0, 0, 0, 1, 0, 0> temperature;
typedef mpl::vector_c<int, 0, 0, 0, 0, 0, 1, 0> intensity;
typedef mpl::vector_c<int, 0, 0, 0, 0, 0, 0, 1> amount_of_substance;

typedef mpl::vector_c<int, 0, 1, -1, 0, 0, 0, 0> velocity;       // l/t
typedef mpl::vector_c<int, 0, 1, -2, 0, 0, 0, 0> acceleration;   // l/(t^2)
typedef mpl::vector_c<int, 1, 1, -1, 0, 0, 0, 0> momentum;       // ml/t
typedef mpl::vector_c<int, 1, 1, -2, 0, 0, 0, 0> force;          // ml/(t^2)
typedef mpl::vector_c<int, 0, 0, 0, 0, 0, 0, 0> scalar;

template<class T, class Dimensions>
struct quantity {
    explicit quantity(T x): m_value(x) {}
    template<class OtherDimensions>
    quantity(quantity<T, OtherDimensions>const& rhs): m_value(rhs.value()) {
        BOOST_STATIC_ASSERT((mpl::equal<Dimensions, OtherDimensions>::type::value));
    }
    T value() const {return m_value;}
    private:
    T m_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 D>
quantity<T, D>
operator-(quantity<T,D> x, quantity<T, D> y) {
    return quantity<T, D>(x.value() - y.value());
}

struct plus_f {
    template<class T1, class T2>
    struct apply: mpl::plus<T1, T2> {};
};
struct minus_f {
    template<class T1, class T2>
    struct apply: mpl::minus<T1, T2> {};
};

template<class T, class D1, class D2>
quantity<T, typename mpl::transform<D1, D2, plus_f>::type >
operator*(quantity<T, D1> x, quantity<T, D2> y) {
    typedef typename mpl::transform<D1, D2, plus_f>::type dim;
    return quantity<T, dim>(x.value() * y.value);
}

template<class T, class D1, class D2>
quantity<T, typename mpl::transform<D1, D2, plus_f>::type >
operator/(quantity<T, D1> x, quantity<T, D2> y) {
    typedef typename mpl::transform<D1, D2, minus_f>::type dim;
    return quantity<T, dim>(x.value() / y.value);
}

By .