# Dimensional-aware scalar and vectorial operations in C++

Eduard George Stan Michele Renda Călin Alexa<EMAIL_ADDRESS>IFIN-HH,
Particles Physics Department, Măgurele, Romania  University of Bucharest,
Faculty of Physics, Măgurele, Romania

###### Abstract

We review the problem of a dimensional check as a language feature for most of
the high-level programming languages and the possible solutions of this issue.
We introduce our C++ library designed to make 2D and 3D vectorial operations
using units of measurement. Eventually, we present our plans for future
improvements of the current implementation.

††journal: Journal of Computational Physics

##  1 Introduction

Dimensional analysis in one of the most powerful tools available in physics to
verify the correctness of mathematical formulas describing physical processes.
However, this tool alone is not able to spot all possible mistakes, it allows
only to perform a fast-check and to easily highlight the faulty operations and
provide a consistence check for the whole process.

Nowadays, complex simulations and derivations are performed using computing
algorithms taking advantage from high-speed calculation and numerical
correctness due to the advances of computing technology. However, when using
complex algorithms, spanning for thousands of source code lines divided in
tens of files, we tend not to automatically trust the results or, much worse,
to blindly accept them. Because there is no bug-free software implementation
we have to mitigate the impact of possible errors using count-less checks in
code and using proper testing procedures, which often are as time-consuming as
the development itself. An appropriate approach is to use development
techniques that could help reduce the overall cost of software development
[1]. Due to historical reasons, most high-levels programming languages lacks
of a dimensional check as a language feature, they work on floats and
integers, leaving to the programmer mercy the responsibility to keep trace of
their semantic meaning, leading to many potential mistakes.

Several libraries try to provide the dimensional analysis check at compile
time or at run-time. However, as shown by Mc.Keever [2], there is a clear
reticence to use such solutions due to:
\- Difficulty to use unfamiliar and complex libraries
\- Fear to add external dependencies
\- Concerns about the performance impact
\- Possible limitations on using external libraries or complex data structures
such as matrices and vectors.

In section 2 we review the existing solutions available in the most common
programming languages, we will describe our proposed solution and concrete
examples in section 3, and we present a performance analysis in section 4.
Limitations of the current development and future plans for improvements are
discussed in section 5 and eventually, conclusions are drawn.

##  2 Existing solutions

For a comprehensive review of dimensional analysis in programming languages we
refer to [3]. Apple’s Swift [4] and F# [5] are two widespread programming
languages that provide native support for UoM-aware (units of measure)
quantities. Promoting UoM analysis to a language feature has the benefit of
immediate access to meaningful error messages. Unfortunately, there are only a
few programming languages that offer native support on this matter, most
languages have an external library which relies on generic programming
techniques.

The most popular solution for C++ is Boost.Units [6], which makes possible the
UoM validation at compile time. Its main advantage is that it is included in
the Boost framework, a comprehensive set of libraries which is used as
launchpad for features that are later included in the standard library of the
language. Unfortunately, Boost is quite often seen as a heavy dependency.

For the FORTAN language there is the PHYSUNITS module [7] which provides
support for dimensional-aware routine creation. While this solution is
versatile, allowing working with both F77 and F90 code, because dimension
information is carried alongside numerical values PHYSUNITS introduces a run-
time execution overhead in term of speed and memory usage.

Dimensional analysis is not limited to compiled languages, Python supports
run-time dimensional analysis check thanks to libraries like Pint [8].

Almost all the existing libraries do not provide native support for vectorial
computation. However, there are two solutions in high energy physics (HEP) for
vectorial operations: CLHEP [9] and Eigen [10]. They provide representations
for both 2D and 3D vectorial quantities and n-dimensional matrices. However,
none of them has built-in support for dimensional-aware units, relying on
primitive numeric types available in C++.

##  3 Our UoM implementation

Integrating Boost.Units, Univec is our C++ library that allows 2D and 3D
vectorial operations. We provide a tool for dimension-aware vector
calculations which has a compile-time check of dimensional correctness.

The source code is available under open-source license
https://gitlab.com/micrenda/univec

Our implementation makes heavy use of generic template programming, shifting
the cost of the dimensional check at compiling time, removing, in principle,
any run-time overhead. The drawback is an increased usage complexity, which
may discourage non-experienced C++ developers. We tried to mitigate this issue
providing a set of default quantity alias, which covers the most common
quantities used in physics. This allows easy-to-remember type names to define
a variable, such as:

⬇

dfpe::QtySiVelocity my_var;

using a shorter name then its canonical name:

⬇

boost::units::quantity<boost::units::si::velocity> my_var;

The full list of aliases defined can be found at the address:
https://gitlab.com/micrenda/qtydef

###  3.1 Usage examples

Using the library we developed, together with our custom type definitions, we
can write compact and semantically clear code with fully benefit of
dimensinal-check at any instant. For example, to calculate the energy of an
object, we can use the well-known formulas:

⬇

VectorC3D<QtySiVelocity> vec1;

\parvec1 = VectorC3D<QtySiVelocity> (

23.5 * si::meter_per_second,

45.0 * si::meter_per_second,

63.0 * si::meter_per_second);

\parQtySiEnergy kinetic_energy =

1./2. * si::kilogram * vec1.normSquared();

\parcout << ”Kinetic energy: ” << kinetic_energy << endl;

\parQtySiLength height(3. * si::metre);

QtySiAcceleration g(9.8 * si::metre_per_second_squared);

QtySiEnergy gravitational_enegry = mass * height * g;

\parcout << ”Gravitational energy: ”

<< gravitational_enegry << endl;

\parQtySiVelocity c(299792458. * si::metre_per_second);

QtySiEnergy einstein_energy = mass * c * c;

\parcout << ”Einstein’s energy: ” << einstein_energy << endl;

however, if we perform any definition mistake, such as:

⬇

// The initialisation of the variable needs to coincide

// with it’s unit of measurement

QtySiLength mass(1. * si::kilogram);

or a plain calculus error, like this one:

⬇

QtySiMass mass(1. * si::kilogram);

\parVectorC3D<QtySiVelocity> vec1 (

23.5 * si::meter_per_second,

45.0 * si::meter_per_second,

63.0 * si::meter_per_second);

\par// This operation needs to include the module

// of the velocity vector squared

QtySiEnergy kinetic_energy = 1./2. * mass * (vec1 + vec1);

\parcout << ”Kinetic energy: ” << kinetic_energy << endl;

we will get a compilation error that will not allow to compile and use this
code (as shown in table 2).

Operation |  Method
---|---
$|\textbf{v}|=0$ |  v.isNull()
$|\textbf{v}|$ |  v.norm()
$|\textbf{v}|^{2}$ |  v.normSquared()
$\textbf{v}\cdot\textbf{u}$ |  v.dot()
$\textbf{v}\times\textbf{u}$ |  v.cross(u)
$\varangle(\textbf{v},\textbf{u})$ |  v.angle(u)
$\textbf{v}/|\textbf{v}|$ |  v.versor()
Table 1: List of functions implemented for each type of vector C++ Code | $2D$ | $3D$ | Quaternion
---|---|---|---
u + v | $(u_{x}+v_{x},u_{y}+v_{y})$ | $(u_{x}+v_{x},u_{y}+v_{y},u_{z}+v_{z})$ | $(u_{a}+v_{a},u_{b}+v_{b},u_{c}+v_{c},u_{d}+v_{d})$
u - v | $(u_{x}-v_{x},u_{y}-v_{y})$ | $(u_{x}-v_{x},u_{y}-v_{y},u_{z}-v_{z})$ | $(u_{a}-v_{a},u_{b}-v_{b},u_{c}-v_{c},u_{d}-v_{d})$
u * v | $(u_{x}\>v_{x},u_{y}\>v_{y})$ | $(u_{x}\>v_{x},u_{y}\>v_{y},u_{z}\>v_{z})$ | ($u_{a}\>v_{a},u_{b}\>v_{b},u_{c}\>v_{c},u_{d}\>v_{d})$
u / v | $(u_{x}/v_{x},u_{y}/v_{y})$ | $(u_{x}/v_{x},u_{y}/v_{y},u_{z}/v_{z})$ | ($u_{a}/v_{a},u_{b}/v_{b},u_{c}/v_{c},u_{d}/v_{d})$
k * v | $(k\>v_{x},k\>v_{y})$ | $(k\>v_{x},k\>v_{y},k\>v_{z})$ | $(k\>v_{a},k\>v_{b},k\>v_{c},k\>v_{d})$
k / v | $(k/v_{x},k/v_{y})$ | $(k/v_{x},k/v_{y},k/v_{z})$ | $(k/v_{a},k/v_{b},k/v_{c},k/v_{d})$
v * k | $(v_{x}\>k,v_{y}\>k)$ | $(v_{x}\>k,v_{y}\>k,v_{z}\>k)$ | $(v_{a}\>k,v_{b}\>k,v_{c}\>k,v_{d}\>k)$
v / k | $(v_{x}/k,v_{y}/k)$ | $(v_{x}/k,v_{y}/k,v_{z}/k)$ | $(v_{a}/k,v_{b}/k,v_{c}/k,v_{d}/k)$
u += v | $(u_{x}+v_{x},u_{y}+v_{y})$ | $(u_{x}+v_{x},u_{y}+v_{y},u_{z}+v_{z})$ | $(u_{a}+v_{a},u_{b}+v_{b},u_{c}+v_{c},u_{d}+v_{d})$
u -= v | $(u_{x}-v_{x},u_{y}-v_{y})$ | $(u_{x}-v_{x},u_{y}-v_{y},u_{z}-v_{z})$ | $(u_{a}-v_{a},u_{b}-v_{b},u_{c}-v_{c},u_{d}-v_{d})$
u *= v | $(u_{x}\>v_{x},u_{y}\>v_{y})$ | $(u_{x}\>v_{x},u_{y}\>v_{y},u_{z}\>v_{z})$ | $(u_{a}\>v_{a},u_{b}\>v_{b},u_{c}\>v_{c},u_{d}\>v_{d})$
u /= v | $(u_{x}/v_{x},u_{y}/v_{y})$ | $(u_{x}/v_{x},u_{y}/v_{y},u_{z}/v_{z})$ | $(u_{a}/v_{a},u_{b}/v_{b},u_{c}/v_{c},u_{d}/v_{d})$
u > v | $u_{x}^{2}+u_{y}^{2}>v_{x}^{2}+v_{y}^{2}$ | $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}>v_{x}^{2}+v_{y}^{2}+v_{z}^{2}$ | $u_{a}^{2}+u_{b}^{2}+u_{c}^{2}+u_{d}^{2}>v_{a}^{2}+v_{b}^{2}+v_{c}^{2}+v_{d}^{2}$
u >= v | $u_{x}^{2}+u_{y}^{2}\geq v_{x}^{2}+v_{y}^{2}$ | $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\geq v_{x}^{2}+v_{y}^{2}+v_{z}^{2}$ | $u_{a}^{2}+u_{b}^{2}+u_{c}^{2}+u_{d}^{2}\geq v_{a}^{2}+v_{b}^{2}+v_{c}^{2}+v_{d}^{2}$
u < v | $u_{x}^{2}+u_{y}^{2}<v_{x}^{2}+v_{y}^{2}$ | $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}<v_{x}^{2}+v_{y}^{2}+v_{z}^{2}$ | $u_{a}^{2}+u_{b}^{2}+u_{c}^{2}+u_{d}^{2}<v_{a}^{2}+v_{b}^{2}+v_{c}^{2}+v_{d}^{2}$
u <= v | $u_{x}^{2}+u_{y}^{2}\leq v_{x}^{2}+v_{y}^{2}$ | $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\leq v_{x}^{2}+v_{y}^{2}+v_{z}^{2}$ | $u_{a}^{2}+u_{b}^{2}+u_{c}^{2}+u_{d}^{2}\leq v_{a}^{2}+v_{b}^{2}+v_{c}^{2}+v_{d}^{2}$
u == v | $u_{x}=v_{x}\textrm{ and }u_{y}=v_{y}$ | $u_{x}=v_{x}\textrm{ and }u_{y}=v_{y}\textrm{ and }u_{z}=v_{z}$ | $u_{a}=v_{a}\textrm{ and }u_{b}=v_{b}\textrm{ and }u_{c}=v_{c},u_{d}=v_{d}$
u != v | $u_{x}\neq v_{x}\textrm{ or }u_{y}\neq v_{y}$ | $u_{x}\neq v_{x}\textrm{ or }u_{y}\neq v_{y}\textrm{ or }u_{z}\neq v_{z}$ | $u_{a}\neq v_{a}\textrm{ or }u_{b}\neq v_{b}\textrm{ or }u_{c}\neq v_{c}\textrm{ or }u_{d}\neq v_{d}$
Table 2: List of operators implemented for each type of $2D$, $3D$-vector and
operations. $u$ and $v$ symbols represents vector quantities while $k$
represent a pure-scalar vector. The operations +=, -=, *= and /= represents
the operator-assignment operators in C++. It is important to notice that u * v
and u / v does represent, respectively, element-wise multiplication and
division, which is not equivalent to scalar or vector product (these
operations can be performed using the dot and cross methods, instead).

⬇

main.cpp:21:13: error: no viable conversion from ’VectorC3D<typename
boost::units::multiply_typeof_helper<…>’ (aka ’VectorC3D<quantity<unit<…>>>’)
to ’dfpe::QtySiEnergy’ (aka ’quantity<unit<…>>’)

QtySiEnergy kinetic_energy = 1./2. * mass * (vec1 + vec1);

^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

/usr/include/boost/units/quantity.hpp:106:25: note: candidate constructor not
viable: no known conversion from ’VectorC3D<typename
boost::units::multiply_typeof_helper<…>’ (aka ’VectorC3D<quantity<unit<…>>>’)
to
’boost::units::quantity<boost::units::unit<…>>>::unspecified_null_pointer_constant_type’
(aka ’void (boost::units::quantity<boost::units::unit<…>>>::*)(int *******)’)
for 1st argument

BOOST_CONSTEXPR quantity(unspecified_null_pointer_constant_type) : val_()

^

/usr/include/boost/units/quantity.hpp:111:25: note: candidate constructor not
viable: no known conversion from ’VectorC3D<typename
boost::units::multiply_typeof_helper<…>’ (aka ’VectorC3D<quantity<unit<…>>>’)
to ’const boost::units::quantity<boost::units::unit<…>>>::this_type &’ (aka
’const quantity<boost::units::unit<…>>, double> &’) for 1st argument

BOOST_CONSTEXPR quantity(const this_type& source) : val_(source.val_)

^

1 error generated.

Example of units mismatch between VectorC3D and scalar value. Parts of this
error were omitted for readability.

###  3.2 Reference frame transformations

While using Univec, we found that formulas are simplified for coordinate
transformation. One clear example comes from the simulation of microscopic
electron-molecule non-relativistic interactions in the center-of-momentum
reference frame. Any linear transformation can be represented by a
transformation-matrix: complex transformations can then be easily represented
and computed by chaining multiple simpler transformations such as
translations, rotations and scaling operations.

Class name |  Description
---|---
TranslateFrameC3D |  Linear translation by a compatible vector
RotateFrameC3D |  Rotation according a Quaternion3D or an EulerRotation3D
ScaleFrameC3D |  Scaling by a scalar or by a different value for each axis
CompositeFrameC3D |  Wrap two transformation in a single one
Table 3: List of ready to use reference frame transformations, that can be
used to describe complex frame transformations. The class CompositeFrameC3D is
rarely used directly but can be created combining the other transformations
using the >> operator.

One can make a frame transformation using a concrete implementation of the
abstract class BaseFrameC3D (in table 3 we present ready-to-use
implementations). One can transform any vector in and out the new reference
frame:

⬇

VectorC3D<QtySiLength> translation (

QtySiLenth(),

QtySiLenth(1. * si::meter),

QtySiLenth());

\parTranslateFrameC3D<QtySiLength> frame(translation);

\parVectorC3D<QtySiLength> globalPosition; // ( 0, 0, 0) m

\parVectorC3D<QtySiLength> localPosition

= frame.forward(globalPosition); // ( 0,-1, 0) m

\parlocalPosition += Vector( QtySiLenth(1. * si::meter),

QtySiLenth(1. * si::meter),

QtySiLenth(1. * si::meter)); // ( 1, 0, 1) m

\parVectorC3D<QtySiLength> globalPosition

= frame.backward(localPosition); // ( 1, 1, 1) m

Multiple transformations can be combined to describe complex scenarios, using
the syntax:

⬇

VectorC3D<QtySiLengh> translation;

translation.y = quantity<si::length>(1. * si::meter);

\parEulerRotation3D<QtySiLength> rotation(

QtySiPlaneAngle(45. * degree::degrees),

QtySiPlaneAngle(),

QtySiPlaneAngle(),

RotationMode::INTRINSIC,

RotationAxis::X,

RotationAxis::Y,

RotationAxis::Z);

\parQtySiDimensionless scale(0.5);

\parauto composite = TranslateFrameC3D<QtySiLength>(translation)

>> RotateFrameC3D<QtySiLength>(rotation)

>> ScaleFrameC3D<QtySiDimensionless>(scale);

##  4 Performance study

Our performance testing consists in comparing three types of methods for
calculating the same operation. The first one is Raw, where the computation
has no dimensional check, using only C++ primitive types such as double. For
the second method Semi, we employ the Boost.Units method without the support
of vectorial operation. The third method is Univec, where our library is
performing UoM-aware vectorial operation. We aim to investigate the
performance loss due to UoM framework.

In fig. 1 we see that while Semi and Univec are quite similar, there is a
significant performance difference between Raw and both Semi and Univec, Raw
has a mean run time of about a third of the other two.

We conclude that the performance loss is coming from the underling UoM library
Boost.Units and not from Univec.

The values shown in fig. 1 are averages obtained for 1000 runs with 10000
random vectors each for the following operations: Norm Squared, Norm, isNull,
Dot Product and Cross Product for a 3D Vector.

Fig. 1: Time needed to perform the specified operation for 10 000 random
vectors. The value displayed is the mean value for 1000 runs, while error bars
represent the standard deviations. The operations are described in table 2, in
three different implementations. The raw implementation uses only primitive
C++ types such as double, and, for each vector, the $x$, $y$ and $z$
components are saved in a separate variable. The semi implementation uses
Boost.Units for UoM validation, but it still stores each vector component in a
separate variable. Finally, the univec implementation uses our UoM-aware
vector library, which employs a more compact syntax and provides results
similar to the semi approach. In this plot lower values means better
performance.

##  5 Limitations and future plans

The use of Boost.Units as underlying library for Univec has several non-
negligible weak points that, unfortunately, seem not to have a solution in the
near future. The first limitation is that this library is available only as
part of a more comprehensive framework, Boost, which is often seen by many
software developers as a heavy dependency. The second limitation is that the
current implementation is based on heavy use of Boost Metaprogramming Library
(MPL) and of template meta-programming techniques, which makes the analysis of
error messages (expecialy on $GCC$ compiler) very hard to understand for non-
experienced developers. In addition, despite the fact that this library is
advertised as zero-overhead, but, as we can observe in fig. 1 that there is
still a non-negligible performance cost for its usage, due to missing compiler
optimizations.

There is a proposal to introduce the support for UoM-aware units directly in
the standard language (see Mateusz Pusz’s talk at CppCon 2019 [11]). This
approach will solve the problems presented above and will provide a more
reliable solution with a less steep learning curve. We plan to switch our
Univec implementation from Boost.Units to any standard library implementation,
when it will become available.

##  6 Conclusions

In this paper we presented Univec, our solution for the issue of UoM-
validation in software development, aiming to improve the heavy vectorial
calculations.

After a brief overview of the existing UoM solutions, we discussed the
limitations which make them unsuitable. We presented an user-case scenario of
our solution and we have estimated the performance.

We plan to move to a native solution when an implementation of a UoM-aware
library will be published in the C++ standard library.

The developing of this library started as an internal project during the
designing phase of our Betaboltz [12] project. While working on this project,
we found the benefits of integrating the UoM analysis directly in our source
code, virtually removing the most usual causes of common mistakes. Even if
this solution does not remove all the sources of errors, Univec allowed us to
focus on the main workflow of the development, increasing the confidence on
our implementation.

We hope that in the future similar solutions will be widely accepted in the
development of software solutions for the scientific community.

## Acknowledgments

This study was supported by PN19060104 and ATLAS CERN-RO grants.

## References

  * Mayerhofer et al. [2016] T. Mayerhofer, M. Wimmer, A. Vallecillo,  Adding uncertainty and units to quantity types in software models,  in: Proceedings of the 2016 ACM SIGPLAN International Conference on Software Language Engineering, SLE 2016, Association for Computing Machinery, New York, NY, USA, 2016, pp. 118–131. URL: https://doi.org/10.1145/2997364.2997376. doi:10.1145/2997364.2997376.
  * McKeever et al. [2021] S. McKeever, O. Bennich-Björkman, O.-A. Salah,  Unit of measurement libraries, their popularity and suitability,  Software: Practice and Experience 51 (2021) 711–734.
  * Preussner [2018] G. M. Preussner, Dimensional Analysis in Programming Languages, 2018. URL: https://gmpreussner.com/research/dimensional-analysis-in-programming-languages.
  * Apple Inc [2022] Apple Inc, Swift Language, 2022\. URL: https://swift.org.
  * Kennedy [2010] A. Kennedy,  Types for Units-of-Measure: Theory and Practice,  in: Z. Horváth, R. Plasmeijer, V. Zsók (Eds.), Central European Functional Programming School: Third Summer School, CEFP 2009, Budapest, Hungary, May 21-23, 2009 and Komárno, Slovakia, May 25-30, 2009, Revised Selected Lectures, Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2010, pp. 268–305. URL: https://doi.org/10.1007/978-3-642-17685-2_8. doi:10.1007/978-3-642-17685-2_8.
  * Schabel and Watanabe [2020] M. C. Schabel, S. Watanabe, Chapter 44. Boost.Units 1.1.0 - 1.77.0, 2020\. URL: https://www.boost.org/doc/libs/1_77_0/doc/html/boost_units.html.
  * Petty [2001] G. W. Petty,  Automated computation and consistency checking of physical dimensions and units in scientific programs,  Software: Practice and Experience 31 (2001) 1067–1076. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/spe.401.
  * Grecco [2022] H. Grecco, Pint, 2022. URL: https://github.com/hgrecco/pint, original-date: 2012-07-13T16:27:11Z.
  * Lönnblad [1994] L. Lönnblad,  CLHEP—a project for designing a C++ class library for high energy physics,  Computer Physics Communications 84 (1994) 307–316.
  * Guennebaud et al. [2010] G. Guennebaud, B. Jacob, et al., Eigen v3, 2010. URL: http://eigen.tuxfamily.org.
  * Pusz [2019] M. Pusz,  A c++ approach to physical units,  2019\. URL: https://isocpp.org/blog/2021/09/cppcon-2019-a-cpp-approach-to-physical-units-mateusz-pusz.
  * Renda et al. [2021] M. Renda, D. Ciubotaru, C. Banu,  Betaboltz: A Monte-Carlo simulation tool for gas scattering processes,  Computer Physics Communications 267 (2021) 108057.