Boost
C++ Libraries
...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
Boost
C++ Libraries
...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
By using MPL metafunctions and the template specializations for operations
on composite dimensions (defined in boost/units/dimension.hpp)
it is possible to perform compile time arithmetic according to the dimensional
analysis rules described above
to produce new composite dimensions :
typedef mpl::times<length_dimension,mass_dimension>::type LM_type; typedef mpl::divides<length_dimension,time_dimension>::type L_T_type; typedef static_root< mpl::divides<energy_dimension,mass_dimension>::type, static_rational<2> >::type V_type;
outputting (with symbol demangling, implemented in utility.hpp)
length_dimension = list<dim<length_base_dimension, static_rational<1l, 1l> >, dimensionless_type> mass_dimension = list<dim<mass_base_dimension, static_rational<1l, 1l> >, dimensionless_type> time_dimension = list<dim<time_base_dimension, static_rational<1l, 1l> >, dimensionless_type> energy_dimension = list<dim<length_base_dimension, static_rational<2l, 1l> >, list<dim<mass_base_dimension, static_rational<1l, 1l> >, list<dim<time_base_dimension, static_rational<-2l, 1l> >, dimensionless_type> > > LM_type = list<dim<length_base_dimension, static_rational<1l, 1l> >, list<dim<mass_base_dimension, static_rational<1l, 1l> >, dimensionless_type> > L_T_type = list<dim<length_base_dimension, static_rational<1l, 1l> >, list<dim<time_base_dimension, static_rational<-1l, 1l> >, dimensionless_type> > V_type = list<dim<length_base_dimension, static_rational<1l, 1l> >, list<dim<time_base_dimension, static_rational<-1l, 1l> >, dimensionless_type> >
(unit.cpp)
This example demonstrates the use of the simple but functional unit system implemented in test_system.hpp
const length L; const mass M; // needs to be namespace-qualified because of global time definition const boost::units::test::time T; const energy E;
We can perform various algebraic operations on these units, resulting in the following output:
L = m L+L = m L-L = m L/L = dimensionless meter*meter = m^2 M*(L/T)*(L/T) = m^2 kg s^-2 M*(L/T)^2 = m^2 kg s^-2 L^3 = m^3 L^(3/2) = m^(3/2) 2vM = kg^(1/2) (3/2)vM = kg^(2/3)
This example demonstrates how to use quantities of our toy unit system :
quantity<length> L = 2.0*meters; // quantity of length quantity<energy> E = kilograms*pow<2>(L/seconds); // quantity of energy
giving us the basic quantity functionality :
L = 2 m L+L = 4 m L-L = 0 m L*L = 4 m^2 L/L = 1 dimensionless L*meter = 2 m^2 kilograms*(L/seconds)*(L/seconds) = 4 m^2 kg s^-2 kilograms*(L/seconds)^2 = 4 m^2 kg s^-2 L^3 = 8 m^3 L^(3/2) = 2.82843 m^(3/2) 2vL = 1.41421 m^(1/2) (3/2)vL = 1.5874 m^(2/3)
As a further demonstration of the flexibility of the system, we replace the
double value type with a std::complex<double> value type (ignoring the question of
the meaningfulness of complex lengths and energies) :
quantity<length,std::complex<double> > L(std::complex<double>(3.0,4.0)*meters); quantity<energy,std::complex<double> > E(kilograms*pow<2>(L/seconds));
and find that the code functions exactly as expected with no additional work,
delegating operations to std::complex<double> and performing the appropriate dimensional
analysis :
L = (3,4) m L+L = (6,8) m L-L = (0,0) m L*L = (-7,24) m^2 L/L = (1,0) dimensionless L*meter = (3,4) m^2 kilograms*(L/seconds)*(L/seconds) = (-7,24) m^2 kg s^-2 kilograms*(L/seconds)^2 = (-7,24) m^2 kg s^-2 L^3 = (-117,44) m^3 L^(3/2) = (2,11) m^(3/2) 2vL = (2,1) m^(1/2) (3/2)vL = (2.38285,1.69466) m^(2/3)
This example provides a fairly extensive set of tests covering most of the
quantity
functionality. It uses the SI unit system defined in boost/units/systems/si.hpp.
If we define a few units and associated quantities,
/// scalar const double s1 = 2; const long x1 = 2; const static_rational<4,3> x2; /// define some units force u1 = newton; energy u2 = joule; /// define some quantities quantity<force> q1(1.0*u1); quantity<energy> q2(2.0*u2);
the various algebraic operations between scalars, units, and quantities give
S1 : 2 X1 : 2 X2 : (4/3) U1 : N U2 : J Q1 : 1 N Q2 : 2 J
Scalar/unit operations :
U1*S1 : 2 N S1*U1 : 2 N U1/S1 : 0.5 N S1/U1 : 2 m^-1 kg^-1 s^2
Unit/unit operations and integral/rational powers of units :
U1+U1 : N U1-U1 : N U1*U1 : m^2 kg^2 s^-4 U1/U1 : dimensionless U1*U2 : m^3 kg^2 s^-4 U1/U2 : m^-1 U1^X : m^2 kg^2 s^-4 X1vU1 : m^(1/2) kg^(1/2) s^-1 U1^X2 : m^(4/3) kg^(4/3) s^(-8/3) X2vU1 : m^(3/4) kg^(3/4) s^(-3/2)
Scalar/quantity operations :
Q1*S1 : 2 N S1*Q1 : 2 N Q1/S1 : 0.5 N S1/Q1 : 2 m^-1 kg^-1 s^2
Unit/quantity operations :
U1*Q1 : 1 m^2 kg^2 s^-4 Q1*U1 : 1 m^2 kg^2 s^-4 U1/Q1 : 1 dimensionless Q1/U1 : 1 dimensionless
Quantity/quantity operations and integral/rational powers of quantities :
+Q1 : 1 N -Q1 : -1 N Q1+Q1 : 2 N Q1-Q1 : 0 N Q1*Q1 : 1 m^2 kg^2 s^-4 Q1/Q1 : 1 dimensionless Q1*Q2 : 2 m^3 kg^2 s^-4 Q1/Q2 : 0.5 m^-1 Q1^X1 : 1 m^2 kg^2 s^-4 X1vQ1 : 1 m^(1/2) kg^(1/2) s^-1 Q1^X2 : 1 m^(4/3) kg^(4/3) s^(-8/3) X2vQ1 : 1 m^(3/4) kg^(3/4) s^(-3/2)
Logical comparison operators are also defined between quantities :
/// check comparison tests quantity<length> l1(1.0*meter), l2(2.0*meters);
giving
l1 == l2 false l1 != l2 true l1 <= l2 true l1 < l2 true l1 >= l2 false l1 > l2 false
Implicit conversion is allowed between dimensionless quantities and their corresponding value types :
/// check implicit unit conversion from dimensionless to value_type const double dimless = (q1/q1);
A generic function for computing mechanical work can be defined that takes force and distance arguments in an arbitrary unit system and returns energy in the same system:
/// the physical definition of work - computed for an arbitrary unit system template<class System,class Y> constexpr quantity<unit<energy_dimension,System>,Y> work(quantity<unit<force_dimension,System>,Y> F, quantity<unit<length_dimension,System>,Y> dx) { return F*dx; }
/// test calcuation of work quantity<force> F(1.0*newton); quantity<length> dx(1.0*meter); quantity<energy> E(work(F,dx));
which functions as expected for SI quantities :
F = 1 N dx = 1 m E = 1 J
The ideal gas law can also be implemented in SI units :
/// the ideal gas law in si units template<class Y> constexpr quantity<si::amount,Y> idealGasLaw(const quantity<si::pressure,Y>& P, const quantity<si::volume,Y>& V, const quantity<si::temperature,Y>& T) { using namespace boost::units::si; using namespace constants::codata; return (P*V/(R*T)); }
/// test ideal gas law quantity<temperature> T = (273.+37.)*kelvin; quantity<pressure> P = 1.01325e5*pascals; quantity<length> r = 0.5e-6*meters; quantity<volume> V = (4.0/3.0)*3.141592*pow<3>(r); quantity<amount> n(idealGasLaw(P,V,T));
with the resulting output :
r = 5e-07 m P = 101325 Pa V = 5.23599e-19 m^3 T = 310 K n = 2.05835e-17 mol R = 8.314472 m^2 kg s^-2 K^-1 mol^-1 (rel. unc. = 1.8e-06)
Trigonometric and inverse trigonometric functions can be implemented for
any unit system that provides an angular base dimension. For radians, these
functions are found in boost/units/cmath.hpp
These behave as one expects, with trigonometric functions taking an angular
quantity and returning a dimensionless quantity, while the inverse trigonometric
functions take a dimensionless quantity and return an angular quantity :
Defining a few angular quantities,
/// test trig stuff quantity<plane_angle> theta = 0.375*radians; quantity<dimensionless> sin_theta = sin(theta); quantity<plane_angle> thetap = asin(sin_theta);
yields
theta = 0.375 rd sin(theta) = 0.366273 dimensionless asin(sin(theta)) = 0.375 rd
Dealing with complex quantities is trivial. Here is the calculation of complex impedance :
quantity<electric_potential,complex_type> v = complex_type(12.5,0.0)*volts; quantity<current,complex_type> i = complex_type(3.0,4.0)*amperes; quantity<resistance,complex_type> z = complex_type(1.5,-2.0)*ohms;
giving
V = (12.5,0) V I = (3,4) A Z = (1.5,-2) Ohm I*Z = (12.5,0) V
User-defined value types that support the appropriate arithmetic operations
are automatically supported as quantity value types. The operators that
are supported by default for quantity value types are unary plus, unary
minus, addition, subtraction, multiplication, division, equal-to, not-equal-to,
less-than, less-or-equal-to, greater-than, and greater-or-equal-to. Support
for rational powers and roots can be added by overloading the power_typeof_helper
and root_typeof_helper
classes. Here we implement a user-defined measurement
class that models a numerical measurement with an associated measurement
error and the appropriate algebra and demonstrates its use as a quantity
value type; the full code is found in measurement.hpp.
Then, defining some measurement
quantity
variables
quantity<length,measurement<double> > u(measurement<double>(1.0,0.0)*meters), w(measurement<double>(4.52,0.02)*meters), x(measurement<double>(2.0,0.2)*meters), y(measurement<double>(3.0,0.6)*meters);
gives
x+y-w = 0.48(+/-0.632772) m w*x = 9.04(+/-0.904885) m^2 x/y = 0.666667(+/-0.149071) dimensionless
If we implement the overloaded helper classes for rational powers and roots then we can also compute rational powers of measurement quantities :
w*y^2/(u*x)^2 = 10.17(+/-3.52328) m^-1 w/(u*x)^(1/2) = 3.19612(+/-0.160431) dimensionless
This example demonstrates the various allowed conversions between SI and CGS units. Defining some quantities
quantity<si::length> L1 = quantity<si::length,int>(int(2.5)*si::meters); quantity<si::length,int> L2(quantity<si::length,double>(2.5*si::meters));
illustrates implicit conversion of quantities of different value types where implicit conversion of the value types themselves is allowed. N.B. The conversion from double to int is treated as an explicit conversion because there is no way to emulate the exact behavior of the built-in conversion. Explicit constructors allow conversions for two cases:
quantity
to a different value_type
:
quantity<si::length,int> L3 = static_cast<quantity<si::length,int> >(L1);
quantity
to a different unit :
quantity<cgs::length> L4 = static_cast<quantity<cgs::length> >(L1);
giving the following output :
L1 = 2 m L2 = 2 m L3 = 2 m L4 = 200 cm L5 = 5 m L6 = 4 m L7 = 200 cm
A few more explicit unit system conversions :
quantity<si::volume> vs(1.0*pow<3>(si::meter)); quantity<cgs::volume> vc(vs); quantity<si::volume> vs2(vc); quantity<si::energy> es(1.0*si::joule); quantity<cgs::energy> ec(es); quantity<si::energy> es2(ec); quantity<si::velocity> v1 = 2.0*si::meters/si::second, v2(2.0*cgs::centimeters/cgs::second);
which produces the following output:
volume (m^3) = 1 m^3 volume (cm^3) = 1e+06 cm^3 volume (m^3) = 1 m^3 energy (joules) = 1 J energy (ergs) = 1e+07 erg energy (joules) = 1 J velocity (2 m/s) = 2 m s^-1 velocity (2 cm/s) = 0.02 m s^-1
This example demonstrates the use of boost::math::quaternion
as a value type for quantity
and the converse. For the first case, we first define specializations of
power_typeof_helper
and root_typeof_helper
for powers and roots, respectively:
/// specialize power typeof helper template<class Y,long N,long D> struct power_typeof_helper<boost::math::quaternion<Y>,static_rational<N,D> > { // boost::math::quaternion only supports integer powers BOOST_STATIC_ASSERT(D==1); typedef boost::math::quaternion< typename power_typeof_helper<Y,static_rational<N,D> >::type > type; static type value(const boost::math::quaternion<Y>& x) { return boost::math::pow(x,static_cast<int>(N)); } };
/// specialize root typeof helper template<class Y,long N,long D> struct root_typeof_helper<boost::math::quaternion<Y>,static_rational<N,D> > { // boost::math::quaternion only supports integer powers BOOST_STATIC_ASSERT(N==1); typedef boost::math::quaternion< typename root_typeof_helper<Y,static_rational<N,D> >::type > type; static type value(const boost::math::quaternion<Y>& x) { return boost::math::pow(x,static_cast<int>(D)); } };
We can now declare a quantity
of a quaternion :
typedef quantity<length,quaternion<double> > length_dimension; length_dimension L(quaternion<double>(4.0,3.0,2.0,1.0)*meters);
so that all operations that are defined in the quaternion
class behave correctly. If rational powers were defined for this class, it
would be possible to compute rational powers and roots with no additional
changes.
+L = (4,3,2,1) m -L = (-4,-3,-2,-1) m L+L = (8,6,4,2) m L-L = (0,0,0,0) m L*L = (2,24,16,8) m^2 L/L = (1,0,0,0) dimensionless L^3 = (-104,102,68,34) m^3
Now, if for some reason we preferred the quantity to be the value
type of the quaternion class
we would have :
typedef quaternion<quantity<length> > length_dimension; length_dimension L(4.0*meters,3.0*meters,2.0*meters,1.0*meters);
Here, the unary plus and minus and addition and subtraction operators function
correctly. Unfortunately, the multiplication and division operations fail
because quaternion implements
them in terms of the *= and
/= operators, respectively,
which are incapable of representing the heterogeneous unit algebra needed
for quantities (an identical problem occurs with std::complex<T>,
for the same reason). In order to compute rational powers and roots, we need
to specialize power_typeof_helper
and root_typeof_helper
as follows:
/// specialize power typeof helper for quaternion<quantity<Unit,Y> > template<class Unit,long N,long D,class Y> struct power_typeof_helper< boost::math::quaternion<quantity<Unit,Y> >, static_rational<N,D> > { typedef typename power_typeof_helper< Y, static_rational<N,D> >::type value_type; typedef typename power_typeof_helper< Unit, static_rational<N,D> >::type unit_type; typedef quantity<unit_type,value_type> quantity_type; typedef boost::math::quaternion<quantity_type> type; static type value(const boost::math::quaternion<quantity<Unit,Y> >& x) { const boost::math::quaternion<value_type> tmp = pow<static_rational<N,D> >(boost::math::quaternion<Y>( x.R_component_1().value(), x.R_component_2().value(), x.R_component_3().value(), x.R_component_4().value())); return type(quantity_type::from_value(tmp.R_component_1()), quantity_type::from_value(tmp.R_component_2()), quantity_type::from_value(tmp.R_component_3()), quantity_type::from_value(tmp.R_component_4())); } };
/// specialize root typeof helper for quaternion<quantity<Unit,Y> > template<class Unit,long N,long D,class Y> struct root_typeof_helper< boost::math::quaternion<quantity<Unit,Y> >, static_rational<N,D> > { typedef typename root_typeof_helper< Y, static_rational<N,D> >::type value_type; typedef typename root_typeof_helper< Unit, static_rational<N,D> >::type unit_type; typedef quantity<unit_type,value_type> quantity_type; typedef boost::math::quaternion<quantity_type> type; static type value(const boost::math::quaternion<quantity<Unit,Y> >& x) { const boost::math::quaternion<value_type> tmp = root<static_rational<N,D> >(boost::math::quaternion<Y>( x.R_component_1().value(), x.R_component_2().value(), x.R_component_3().value(), x.R_component_4().value())); return type(quantity_type::from_value(tmp.R_component_1()), quantity_type::from_value(tmp.R_component_2()), quantity_type::from_value(tmp.R_component_3()), quantity_type::from_value(tmp.R_component_4())); } };
giving:
+L = (4 m,3 m,2 m,1 m) -L = (-4 m,-3 m,-2 m,-1 m) L+L = (8 m,6 m,4 m,2 m) L-L = (0 m,0 m,0 m,0 m) L^3 = (-104 m^3,102 m^3,68 m^3,34 m^3)
This example demonstrates how to implement a replacement complex
class that functions correctly both as a quantity value type and as a quantity
container class, including heterogeneous multiplication and division operations
and rational powers and roots. Naturally, heterogeneous operations are only
supported on compilers that implement typeof.
The primary differences are that binary operations are not implemented using
the op=
operators and use the utility classes add_typeof_helper,
subtract_typeof_helper,
multiply_typeof_helper,
and divide_typeof_helper.
In addition, power_typeof_helper
and root_typeof_helper
are defined for both cases :
namespace boost { namespace units { /// replacement complex class template<class T> class complex { public: typedef complex<T> this_type; constexpr complex(const T& r = 0,const T& i = 0) : r_(r),i_(i) { } constexpr complex(const this_type& source) : r_(source.r_),i_(source.i_) { } constexpr this_type& operator=(const this_type& source) { if (this == &source) return *this; r_ = source.r_; i_ = source.i_; return *this; } constexpr T& real() { return r_; } constexpr T& imag() { return i_; } constexpr const T& real() const { return r_; } constexpr const T& imag() const { return i_; } constexpr this_type& operator+=(const T& val) { r_ += val; return *this; } constexpr this_type& operator-=(const T& val) { r_ -= val; return *this; } constexpr this_type& operator*=(const T& val) { r_ *= val; i_ *= val; return *this; } constexpr this_type& operator/=(const T& val) { r_ /= val; i_ /= val; return *this; } constexpr this_type& operator+=(const this_type& source) { r_ += source.r_; i_ += source.i_; return *this; } constexpr this_type& operator-=(const this_type& source) { r_ -= source.r_; i_ -= source.i_; return *this; } constexpr this_type& operator*=(const this_type& source) { *this = *this * source; return *this; } constexpr this_type& operator/=(const this_type& source) { *this = *this / source; return *this; } private: T r_,i_; }; } } #if BOOST_UNITS_HAS_BOOST_TYPEOF #include BOOST_TYPEOF_INCREMENT_REGISTRATION_GROUP() BOOST_TYPEOF_REGISTER_TEMPLATE(boost::units::complex, 1) #endif namespace boost { namespace units { template<class X> constexpr complex<typename unary_plus_typeof_helper<X>::type> operator+(const complex<X>& x) { typedef typename unary_plus_typeof_helper<X>::type type; return complex<type>(x.real(),x.imag()); } template<class X> constexpr complex<typename unary_minus_typeof_helper<X>::type> operator-(const complex<X>& x) { typedef typename unary_minus_typeof_helper<X>::type type; return complex<type>(-x.real(),-x.imag()); } template<class X,class Y> constexpr complex<typename add_typeof_helper<X,Y>::type> operator+(const complex<X>& x,const complex<Y>& y) { typedef typename boost::units::add_typeof_helper<X,Y>::type type; return complex<type>(x.real()+y.real(),x.imag()+y.imag()); } template<class X,class Y> constexpr complex<typename boost::units::subtract_typeof_helper<X,Y>::type> operator-(const complex<X>& x,const complex<Y>& y) { typedef typename boost::units::subtract_typeof_helper<X,Y>::type type; return complex<type>(x.real()-y.real(),x.imag()-y.imag()); } template<class X,class Y> constexpr complex<typename boost::units::multiply_typeof_helper<X,Y>::type> operator*(const complex<X>& x,const complex<Y>& y) { typedef typename boost::units::multiply_typeof_helper<X,Y>::type type; return complex<type>(x.real()*y.real() - x.imag()*y.imag(), x.real()*y.imag() + x.imag()*y.real()); // fully correct implementation has more complex return type // // typedef typename boost::units::multiply_typeof_helper<X,Y>::type xy_type; // // typedef typename boost::units::add_typeof_helper< // xy_type,xy_type>::type xy_plus_xy_type; // typedef typename // boost::units::subtract_typeof_helper<xy_type,xy_type>::type // xy_minus_xy_type; // // BOOST_STATIC_ASSERT((boost::is_same<xy_plus_xy_type, // xy_minus_xy_type>::value == true)); // // return complex<xy_plus_xy_type>(x.real()*y.real()-x.imag()*y.imag(), // x.real()*y.imag()+x.imag()*y.real()); } template<class X,class Y> constexpr complex<typename boost::units::divide_typeof_helper<X,Y>::type> operator/(const complex<X>& x,const complex<Y>& y) { // naive implementation of complex division typedef typename boost::units::divide_typeof_helper<X,Y>::type type; return complex<type>((x.real()*y.real()+x.imag()*y.imag())/ (y.real()*y.real()+y.imag()*y.imag()), (x.imag()*y.real()-x.real()*y.imag())/ (y.real()*y.real()+y.imag()*y.imag())); // fully correct implementation has more complex return type // // typedef typename boost::units::multiply_typeof_helper<X,Y>::type xy_type; // typedef typename boost::units::multiply_typeof_helper<Y,Y>::type yy_type; // // typedef typename boost::units::add_typeof_helper<xy_type, xy_type>::type // xy_plus_xy_type; // typedef typename boost::units::subtract_typeof_helper< // xy_type,xy_type>::type xy_minus_xy_type; // // typedef typename boost::units::divide_typeof_helper< // xy_plus_xy_type,yy_type>::type xy_plus_xy_over_yy_type; // typedef typename boost::units::divide_typeof_helper< // xy_minus_xy_type,yy_type>::type xy_minus_xy_over_yy_type; // // BOOST_STATIC_ASSERT((boost::is_same<xy_plus_xy_over_yy_type, // xy_minus_xy_over_yy_type>::value == true)); // // return complex<xy_plus_xy_over_yy_type>( // (x.real()*y.real()+x.imag()*y.imag())/ // (y.real()*y.real()+y.imag()*y.imag()), // (x.imag()*y.real()-x.real()*y.imag())/ // (y.real()*y.real()+y.imag()*y.imag())); } template<class Y> complex<Y> pow(const complex<Y>& x,const Y& y) { std::complex<Y> tmp(x.real(),x.imag()); tmp = std::pow(tmp,y); return complex<Y>(tmp.real(),tmp.imag()); } template<class Y> std::ostream& operator<<(std::ostream& os,const complex<Y>& val) { os << val.real() << " + " << val.imag() << " i"; return os; } /// specialize power typeof helper for complex<Y> template<class Y,long N,long D> struct power_typeof_helper<complex<Y>,static_rational<N,D> > { typedef complex< typename power_typeof_helper<Y,static_rational<N,D> >::type > type; static type value(const complex<Y>& x) { const static_rational<N,D> rat; const Y m = Y(rat.numerator())/Y(rat.denominator()); return boost::units::pow(x,m); } }; /// specialize root typeof helper for complex<Y> template<class Y,long N,long D> struct root_typeof_helper<complex<Y>,static_rational<N,D> > { typedef complex< typename root_typeof_helper<Y,static_rational<N,D> >::type > type; static type value(const complex<Y>& x) { const static_rational<N,D> rat; const Y m = Y(rat.denominator())/Y(rat.numerator()); return boost::units::pow(x,m); } }; /// specialize power typeof helper for complex<quantity<Unit,Y> > template<class Y,class Unit,long N,long D> struct power_typeof_helper<complex<quantity<Unit,Y> >,static_rational<N,D> > { typedef typename power_typeof_helper<Y,static_rational<N,D> >::type value_type; typedef typename power_typeof_helper<Unit,static_rational<N,D> >::type unit_type; typedef quantity<unit_type,value_type> quantity_type; typedef complex<quantity_type> type; static type value(const complex<quantity<Unit,Y> >& x) { const complex<value_type> tmp = pow<static_rational<N,D> >(complex<Y>(x.real().value(), x.imag().value())); return type(quantity_type::from_value(tmp.real()), quantity_type::from_value(tmp.imag())); } }; /// specialize root typeof helper for complex<quantity<Unit,Y> > template<class Y,class Unit,long N,long D> struct root_typeof_helper<complex<quantity<Unit,Y> >,static_rational<N,D> > { typedef typename root_typeof_helper<Y,static_rational<N,D> >::type value_type; typedef typename root_typeof_helper<Unit,static_rational<N,D> >::type unit_type; typedef quantity<unit_type,value_type> quantity_type; typedef complex<quantity_type> type; static type value(const complex<quantity<Unit,Y> >& x) { const complex<value_type> tmp = root<static_rational<N,D> >(complex<Y>(x.real().value(), x.imag().value())); return type(quantity_type::from_value(tmp.real()), quantity_type::from_value(tmp.imag())); } }; } // namespace units } // namespace boost
With this replacement complex
class, we can declare a complex variable :
typedef quantity<length,complex<double> > length_dimension; const length_dimension L(complex<double>(2.0,1.0)*meters);
to get the correct behavior for all cases supported by quantity with a complex value type :
+L = 2 + 1 i m -L = -2 + -1 i m L+L = 4 + 2 i m L-L = 0 + 0 i m L*L = 3 + 4 i m^2 L/L = 1 + 0 i dimensionless L^3 = 2 + 11 i m^3 L^(3/2) = 2.56713 + 2.14247 i m^(3/2) 3vL = 1.29207 + 0.201294 i m^(1/3) (3/2)vL = 1.62894 + 0.520175 i m^(2/3)
and, similarly, complex with
a quantity
value type
typedef complex<quantity<length> > length_dimension; const length_dimension L(2.0*meters,1.0*meters);
gives
+L = 2 m + 1 m i -L = -2 m + -1 m i L+L = 4 m + 2 m i L-L = 0 m + 0 m i L*L = 3 m^2 + 4 m^2 i L/L = 1 dimensionless + 0 dimensionless i L^3 = 2 m^3 + 11 m^3 i L^(3/2) = 2.56713 m^(3/2) + 2.14247 m^(3/2) i 3vL = 1.29207 m^(1/3) + 0.201294 m^(1/3) i (3/2)vL = 1.62894 m^(2/3) + 0.520175 m^(2/3) i
This example provides an ad hoc performance test to verify that zero runtime
overhead is incurred when using quantity
in place of double. Note that
performance optimization and testing is not trivial, so some care must be
taken in profiling. It is also critical to have a compiler capable of optimizing
the many template instantiations and inline calls effectively to achieve
maximal performance. Zero overhead for this test has been verified using
gcc 4.0.1, and icc 9.0, 10.0, and 10.1 on Mac OS 10.4 and 10.5, and using
msvc 8.0 on Windows XP.
This example demonstrates the implementation of two non-SI units of length, the nautical mile :
namespace nautical { struct length_base_unit : boost::units::base_unit<length_base_unit, length_dimension, 1> { static std::string name() { return "nautical mile"; } static std::string symbol() { return "nmi"; } }; typedef boost::units::make_system<length_base_unit>::type system; /// unit typedefs typedef unit<length_dimension,system> length; static const length mile,miles; } // namespace nautical // helper for conversions between nautical length and si length BOOST_UNITS_DEFINE_CONVERSION_FACTOR(nautical::length_base_unit, boost::units::si::meter_base_unit, double, 1.852e3);
and the imperial foot :
namespace imperial { struct length_base_unit : boost::units::base_unit<length_base_unit, length_dimension, 2> { static std::string name() { return "foot"; } static std::string symbol() { return "ft"; } }; typedef boost::units::make_system<length_base_unit>::type system; /// unit typedefs typedef unit<length_dimension,system> length; static const length foot,feet; } // imperial // helper for conversions between imperial length and si length BOOST_UNITS_DEFINE_CONVERSION_FACTOR(imperial::length_base_unit, boost::units::si::meter_base_unit, double, 1.0/3.28083989501312);
These units include conversions between themselves and the meter. Three functions for computing radar beam height from radar range and the local earth radius are defined. The first takes arguments in one system and returns a value in the same system :
template<class System,typename T> constexpr quantity<unit<boost::units::length_dimension,System>,T> radar_beam_height(const quantity<unit<length_dimension,System>,T>& radar_range, const quantity<unit<length_dimension,System>,T>& earth_radius, T k = 4.0/3.0) { return quantity<unit<length_dimension,System>,T> (pow<2>(radar_range)/(2.0*k*earth_radius)); }
The second is similar, but is templated on return type, so that the arguments are converted to the return unit system internally :
template<class return_type,class System1,class System2,typename T> constexpr return_type radar_beam_height(const quantity<unit<length_dimension,System1>,T>& radar_range, const quantity<unit<length_dimension,System2>,T>& earth_radius, T k = 4.0/3.0) { // need to decide which system to use for calculation return pow<2>(static_cast<return_type>(radar_range)) / (2.0*k*static_cast<return_type>(earth_radius)); }
Finally, the third function is an empirical approximation that is only valid for radar ranges specified in nautical miles, returning beam height in feet. This function uses the heterogeneous unit of nautical miles per square root of feet to ensure dimensional correctness :
constexpr quantity<imperial::length> radar_beam_height(const quantity<nautical::length>& range) { return quantity<imperial::length> (pow<2>(range/(1.23*nautical::miles/root<2>(imperial::feet)))); }
With these, we can compute radar beam height in various unit systems :
const quantity<nautical::length> radar_range(300.0*miles); const quantity<si::length> earth_radius(6371.0087714*kilo*meters); const quantity<si::length> beam_height_1(radar_beam_height(quantity<si::length>(radar_range),earth_radius)); const quantity<nautical::length> beam_height_2(radar_beam_height(radar_range,quantity<nautical::length>(earth_radius))); const quantity<si::length> beam_height_3(radar_beam_height< quantity<si::length> >(radar_range,earth_radius)); const quantity<nautical::length> beam_height_4(radar_beam_height< quantity<nautical::length> >(radar_range,earth_radius));
giving
radar range : 300 nmi earth radius : 6.37101e+06 m beam height 1 : 18169.7 m beam height 2 : 9.81085 nmi beam height 3 : 18169.7 m beam height 4 : 9.81085 nmi beam height approx : 59488.4 ft beam height approx : 18132.1 m
Mixed units and mixed unit conversions.
This code:
quantity<si::length> L(1.5*si::meter); quantity<cgs::mass> M(1.0*cgs::gram); std::cout << L << std::endl << M << std::endl << L*M << std::endl << L/M << std::endl << std::endl; std::cout << 1.0*si::meter*si::kilogram/pow<2>(si::second) << std::endl << 1.0*si::meter*si::kilogram/pow<2>(si::second)/si::meter << std::endl << std::endl; std::cout << 1.0*cgs::centimeter*si::kilogram/pow<2>(si::second) << std::endl << 1.0*cgs::centimeter*si::kilogram/pow<2>(si::second)/si::meter << std::endl << std::endl;
gives
1.5 m 1 g 1.5 m g 1.5 m g^-1 1 N 1 kg s^-2 1 cm kg s^-2 1 cm m^-1 kg s^-2
Arbitrary conversions also work:
quantity<si::area> A(1.5*si::meter*cgs::centimeter); std::cout << 1.5*si::meter*cgs::centimeter << std::endl << A << std::endl << std::endl;
yielding
1.5 cm m 0.015 m^2
This example demonstrates using of absolute temperatures and relative temperature
differences in Fahrenheit and converting between these and the Kelvin temperature
scale. This issue touches on some surprisingly deep mathematical concepts
(see Wikipedia
for a basic review), but for our purposes here, we will simply observe that
it is important to be able to differentiate between an absolute temperature
measurement and a measurement of temperature difference. This is accomplished
by using the absolute
wrapper class.
First we define a system using the predefined fahrenheit base unit:
typedef temperature::fahrenheit_base_unit::unit_type temperature; typedef get_system<temperature>::type system; BOOST_UNITS_STATIC_CONSTANT(degree,temperature); BOOST_UNITS_STATIC_CONSTANT(degrees,temperature);
Now we can create some quantities:
quantity<absolute<fahrenheit::temperature> > T1p( 32.0*absolute<fahrenheit::temperature>()); quantity<fahrenheit::temperature> T1v( 32.0*fahrenheit::degrees); quantity<absolute<si::temperature> > T2p(T1p); quantity<si::temperature> T2v(T1v);
Note the use of absolute
to wrap a unit. The resulting output is:
{ 32 } F { 273.15 } K { 273.15 } K [ 32 ] F [ 17.7778 ] K [ 17.7778 ] K
(runtime_conversion_factor.cpp)
The Boost.Units library does not require that the conversion factors be compile time constants, as is demonstrated in this example:
using boost::units::base_dimension; using boost::units::base_unit; static const long currency_base = 1; struct currency_base_dimension : base_dimension<currency_base_dimension, 1> {}; typedef currency_base_dimension::dimension_type currency_type; template<long N> struct currency_base_unit : base_unit<currency_base_unit<N>, currency_type, currency_base + N> {}; typedef currency_base_unit<0> us_dollar_base_unit; typedef currency_base_unit<1> euro_base_unit; typedef us_dollar_base_unit::unit_type us_dollar; typedef euro_base_unit::unit_type euro; // an array of all possible conversions double conversion_factors[2][2] = { {1.0, 1.0}, {1.0, 1.0} }; double get_conversion_factor(long from, long to) { return(conversion_factors[from][to]); } void set_conversion_factor(long from, long to, double value) { conversion_factors[from][to] = value; conversion_factors[to][from] = 1.0 / value; } BOOST_UNITS_DEFINE_CONVERSION_FACTOR_TEMPLATE((long N1)(long N2), currency_base_unit<N1>, currency_base_unit<N2>, double, get_conversion_factor(N1, N2));
It is also possible to define base units that have derived rather than base dimensions:
struct imperial_gallon_tag : base_unit<imperial_gallon_tag, volume_dimension, 1> { }; typedef make_system<imperial_gallon_tag>::type imperial; typedef unit<volume_dimension,imperial> imperial_gallon; struct us_gallon_tag : base_unit<us_gallon_tag, volume_dimension, 2> { }; typedef make_system<us_gallon_tag>::type us; typedef unit<volume_dimension,us> us_gallon;
If a unit has a special name and/or symbol, the free functions name_string and symbol_string
can be overloaded directly.
std::string name_string(const cgs::force&) { return "dyne"; } std::string symbol_string(const cgs::force&) { return "dyn"; }
In this case, any unit that reduces to the overloaded unit will be output with the replacement symbol.
Special names and symbols for the SI and CGS unit systems are found in boost/units/systems/si/io.hpp
and boost/units/systems/cgs/io.hpp,
respectively. If these headers are not included, the output will simply follow
default rules using the appropriate fundamental dimensions. Note that neither
of these functions is defined for quantities because doing so would require
making assumptions on how the corresponding value type should be formatted.
Three ostream formatters,
symbol_format, name_format, and typename_format
are provided for convenience. These select the textual representation of
units provided by symbol_string
or name_string in the first
two cases, while the latter returns a demangled typename for debugging purposes.
Formatting of scaled unit is also done correctly.
It is often desirable to scale a unit
automatically, depending on its value, to keep the integral part in a limited
range, usually between 1 and 999.
For example, using engineering notation prefixes,
"1234.5 m" is more helpfully displayed as "1.234 km" "0.000000001234 m" is more clearly displayed as "1.2345 nanometer".
The iostream manipulators engineering_prefixes
or binary_prefixes make this
easy.
using boost::units::binary_prefix; using boost::units::engineering_prefix; using boost::units::no_prefix; quantity<length> l = 2.345 * meters; // A quantity of length, in units of meters. cout << engineering_prefix << l << endl; // Outputs "2.345 m". l = 1000.0 * l; // Increase it by 1000, so expect a k prefix. // Note that a double 1000.0 is required - an integer will fail to compile. cout << engineering_prefix << l << endl; // Output autoprefixed with k to "2.345 km". quantity<energy> e = kilograms * pow<2>(l / seconds); // A quantity of energy. cout << engineering_prefix << e << endl; // 5.49902 MJ cout << name_format << engineering_prefix << e << endl; // 5.49902 megaJoule
(The complete set of engineering and scientific multiples is not used (not centi or deci for example), but only powers of ten that are multiples of three, 10^3).
Similarly, the equivalent binary prefixes used for displaying computing kilobytes, megabytes, gigabytes...
These are the 2^10 = 1024, 2^20 = 1 048 576, 2^30 ... multiples.
(See also Prefixes for binary multiples
This scale is specified in IEC 60027-2, Second edition, 2000-11, Letter symbols to be used in electrical technology - Part 2: Telecommunications and electronics).
// Don't forget that the units name or symbol format specification is persistent. cout << symbol_format << endl; // Resets the format to the default symbol format. quantity<byte_base_unit::unit_type> b = 2048. * byte_base_unit::unit_type(); cout << engineering_prefix << b << endl; // 2.048 kb cout << symbol_format << binary_prefix << b << endl; // "2 Kib"
But note that scalar dimensionless values, like int, float and double, are not prefixed automatically by the engineering_prefix or binary_prefix iostream manipulators.
const double s1 = 2345.6; const long x1 = 23456; cout << engineering_prefix << s1 << endl; // 2345.6 cout << engineering_prefix << x1 << endl; // 23456 cout << binary_prefix << s1 << endl; // 2345.6 cout << binary_prefix << x1 << endl; // 23456
You can output the name or symbol of a unit (rather than the most common quantity of a unit).
const length L; // A unit of length (but not a quantity of length). cout << L << endl; // Default length unit is meter, // but default is symbol format so output is just "m". cout << name_format << L << endl; // default length name is "meter".
Note too that all the formatting flags are persistent, so that if you set engineering_prefix, then it applies to all future outputs, until you select binary_prefix, or explicitly switch autoprefix off. You can specify no prefix (the default of course) in two ways:
no_prefix(cout); // Clear any prefix flag. cout << no_prefix << endl; // Clear any prefix flag using `no_prefix` manipulator.
And you can get the format flags for diagnosing problems.
cout << boost::units::get_autoprefix(cout) << endl; // 8 is `autoprefix_binary` from `enum autoprefix_mode`. cout << boost::units::get_format(cout) << endl; // 1 is `name_fmt` from `enum format_mode`.
This code demonstrates the use of the conversion_factor
free function to determine the scale factor between two units.
double dyne_to_newton = conversion_factor(cgs::dyne,si::newton); std::cout << dyne_to_newton << std::endl; double force_over_mass_conversion = conversion_factor(si::newton/si::kilogram,cgs::dyne/cgs::gram); std::cout << force_over_mass_conversion << std::endl; double momentum_conversion = conversion_factor(cgs::momentum(),si::momentum()); std::cout << momentum_conversion << std::endl; double momentum_over_mass_conversion = conversion_factor(si::momentum()/si::mass(),cgs::momentum()/cgs::gram); std::cout << momentum_over_mass_conversion << std::endl; double acceleration_conversion = conversion_factor(cgs::gal,si::meter_per_second_squared); std::cout << acceleration_conversion << std::endl;
Produces
1e-005 100 1e-005 100 0.01
This example shows how to implement an interface that allow different units at runtime while still maintaining type safety for internal calculations.
namespace { using namespace boost::units; using imperial::foot_base_unit; std::map<std::string, quantity<si::length> > known_units; } quantity<si::length> calculate(const quantity<si::length>& t) { return(boost::units::hypot(t, 2.0 * si::meters)); } int main() { known_units["meter"] = 1.0 * si::meters; known_units["centimeter"] = .01 * si::meters; known_units["foot"] = conversion_factor(foot_base_unit::unit_type(), si::meter) * si::meter; std::string output_type("meter"); std::string input; while((std::cout << "> ") && (std::cin >> input)) { if(!input.empty() && input[0] == '#') { std::getline(std::cin, input); } else if(input == "exit") { break; } else if(input == "help") { std::cout << "type \"exit\" to exit\n" "type \"return 'unit'\" to set the return units\n" "type \"'number' 'unit'\" to do a simple calculation" << std::endl; } else if(input == "return") { if(std::cin >> input) { if(known_units.find(input) != known_units.end()) { output_type = input; std::cout << "Done." << std::endl; } else { std::cout << "Unknown unit \"" << input << "\"" << std::endl; } } else { break; } } else { try { double value = boost::lexical_cast<double>(input); if(std::cin >> input) { if(known_units.find(input) != known_units.end()) { std::cout << static_cast<double>( calculate(value * known_units[input]) / known_units[output_type]) << ' ' << output_type << std::endl; } else { std::cout << "Unknown unit \"" << input << "\"" << std::endl; } } else { break; } } catch(...) { std::cout << "Input error" << std::endl; } } } }
The header boost/units/lambda.hpp
provides overloads and specializations needed to make Boost.Units usable
with the Boost.Lambda library.
int main(int argc, char **argv) { using namespace std; namespace bl = boost::lambda; namespace bu = boost::units; namespace si = boost::units::si; //////////////////////////////////////////////////////////////////////// // Mechanical example: linear accelerated movement //////////////////////////////////////////////////////////////////////// // Initial condition variables for acceleration, speed, and displacement bu::quantity<si::acceleration> a = 2.0 * si::meters_per_second_squared; bu::quantity<si::velocity> v = 1.0 * si::meters_per_second; bu::quantity<si::length> s0 = 0.5 * si::meter; // Displacement over time boost::function<bu::quantity<si::length> (bu::quantity<si::time>) > s = 0.5 * bl::var(a) * bl::_1 * bl::_1 + bl::var(v) * bl::_1 + bl::var(s0); cout << "Linear accelerated movement:" << endl << "a = " << a << ", v = " << v << ", s0 = " << s0 << endl << "s(1.0 * si::second) = " << s(1.0 * si::second) << endl << endl; // Change initial conditions a = 1.0 * si::meters_per_second_squared; v = 2.0 * si::meters_per_second; s0 = -1.5 * si::meter; cout << "a = " << a << ", v = " << v << ", s0 = " << s0 << endl << "s(1.0 * si::second) = " << s(1.0 * si::second) << endl << endl; //////////////////////////////////////////////////////////////////////// // Electrical example: oscillating current //////////////////////////////////////////////////////////////////////// // Constants for the current amplitude, frequency, and offset current const bu::quantity<si::current> iamp = 1.5 * si::ampere; const bu::quantity<si::frequency> f = 1.0e3 * si::hertz; const bu::quantity<si::current> i0 = 0.5 * si::ampere; // The invocation of the sin function needs to be postponed using // bind to specify the oscillation function. A lengthy static_cast // to the function pointer referencing boost::units::sin() is needed // to avoid an "unresolved overloaded function type" error. boost::function<bu::quantity<si::current> (bu::quantity<si::time>) > i = iamp * bl::bind(static_cast<bu::dimensionless_quantity<si::system, double>::type (*)(const bu::quantity<si::plane_angle>&)>(bu::sin), 2.0 * pi * si::radian * f * bl::_1) + i0; cout << "Oscillating current:" << endl << "iamp = " << iamp << ", f = " << f << ", i0 = " << i0 << endl << "i(1.25e-3 * si::second) = " << i(1.25e-3 * si::second) << endl << endl; //////////////////////////////////////////////////////////////////////// // Geometric example: area calculation for a square //////////////////////////////////////////////////////////////////////// // Length constant const bu::quantity<si::length> l = 1.5 * si::meter; // Again an ugly static_cast is needed to bind pow<2> to the first // function argument. boost::function<bu::quantity<si::area> (bu::quantity<si::length>) > A = bl::bind(static_cast<bu::quantity<si::area> (*)(const bu::quantity<si::length>&)>(bu::pow<2>), bl::_1); cout << "Area of a square:" << endl << "A(" << l <<") = " << A(l) << endl << endl; //////////////////////////////////////////////////////////////////////// // Thermal example: temperature difference of two absolute temperatures //////////////////////////////////////////////////////////////////////// // Absolute temperature constants const bu::quantity<bu::absolute<si::temperature> > Tref = 273.15 * bu::absolute<si::temperature>(); const bu::quantity<bu::absolute<si::temperature> > Tamb = 300.00 * bu::absolute<si::temperature>(); boost::function<bu::quantity<si::temperature> (bu::quantity<bu::absolute<si::temperature> >, bu::quantity<bu::absolute<si::temperature> >)> dT = bl::_2 - bl::_1; cout << "Temperature difference of two absolute temperatures:" << endl << "dT(" << Tref << ", " << Tamb << ") = " << dT(Tref, Tamb) << endl << endl; return 0; }