UTL

Collection of self-contained header-only libraries for C++17

View on GitHub

utl::math

<- to README.md

<- to implementation.hpp

math adds various helper functions that aim to reduce code verbosity when working with mathematical expressions.

Definitions

// Constants
double PI           = 3.14159265358979323846;
double PI_TWO       = 2. * PI;
double PI_HALF      = 0.5 * PI;
double E            = 2.71828182845904523536;
double GOLDEN_RATIO = 1.6180339887498948482;

// Type traits
template<typename Type>
struct is_addable_with_itself;

template<typename Type>
struct is_multipliable_by_scalar;

template <typename Type>
struct is_sized;

template <typename FuncType, typename Signature>
struct is_function_with_signature;

// Basic math functions
template<typename Type>
Type abs(Type x);

template<typename Type>
Type sign(Type x);

template<typename Type>
Type sqr(Type x); // x^2

template<typename Type>
Type cube(Type x); // x^3

template<typename Type>
Type midpoint(Type a, Type b); // (a + b) / 2

template<typename IntegerType>
int kronecker_delta(IntegerType i, IntegerType j); // (i == j) ? 1 : 0

template<typename IntegerType>
int power_of_minus_one(IntegerType power); // (-1)^power

// Degrees and radians
template<typename FloatType>
FloatType deg_to_rad(FloatType degrees);

template<typename Type>
FloatType rad_to_deg(FloatType radians);

// Meshing
struct Points {
    explicit Points(std::size_t count);
};

struct Intervals {
    explicit Intervals(std::size_t count);
    Intervals(Points points);
};

template<typename FloatType>
std::vector<FloatType> linspace(FloatType L1, FloatType L2, Intervals N);

template<typename FloatType, typename FuncType>
FloatType integrate_trapezoidal(FuncType f, FloatType L1, FloatType L2, Intervals N);

// Other utils
template<typename UintType>
UintType uint_difference(UintType a, UintType b);

template<typename SizedContainer>
int ssize(const SizedContainer& container);
        
// Branchless ternary
template<typename Type>
ArithmeticType ternary_branchless(bool condition, ArithmeticType return_if_true, ArithmeticType return_if_false);

template<typename IntType>
IntType ternary_bitselect(bool condition, IntType return_if_true, IntType return_if_false);

template<typename IntType>
IntType ternary_bitselect(bool condition, IntType return_if_true); // return_if_false == 0

// Memory units
enum class MemoryUnit { BYTE, KiB, MiB, GiB, TiB, KB, MB, GB, TB };

template<typename T, MemoryUnit units = MemoryUnit::MiB>
constexpr double memory_size(std::size_t count);

All methods have appropriate enable_if<> conditions and constexpr qualifiers, which are omitted in documentation for reduced verbosity.

Methods that deal with floating-point values require explicitly floating-point inputs for mathematical strictness.

Methods

Type traits

math::is_addable_with_itself

is_addable_with_itself<Type>::value returns at compile time whether Type objects can be added with operator+().

(see and <type_traits> and UnaryTypeTrait)

math::is_multipliable_by_scalar

is_multipliable_by_scalar<Type>::value returns at compile time whether Type objects can be multiplied by a floating point scalar with operator*().

math::is_sized

is_sized<Type>::value returns at compile time whether Type objects support .size() method.

math::is_function_with_signature<FuncType, Signature>

is_function_with_signature<FuncType, Signature>::value returns at compile time whether FuncType is a callable with signature Signature.

Useful when creating functions that accept callable type as a template argument, rather than std::function. This is usually done to avoid overhead introduced by std::function type erasure, however doing so removes explicit requirements imposed on a callable signature. Using this type trait in combination with std::enable_if allows template approach to keep explicit signature restrictions and overload method for multiple callable types.

Basic math functions

Type math::abs(Type x);
Type math::sign(Type x);
Type math::sqr(Type x);
Type math::cube(Type x);
Returns $ x $, $\mathrm{sign} (x)$, $x^2$ or $x^3$ of an appropriate type.
Type math::midpoint(Type a, Type b);

Returns $\dfrac{a + b}{2}$ of an appropriate type. Can be used with vectors or other custom types that have defined operator+() and scalar operator*().

int math::kronecker_delta(IntegerType i, IntegerType j);

Computes Kronecker delta symbol: $\delta_{ii} = \begin{cases}1, \quad i = j\0, \quad i \neq j\end{cases}$.

int math::power_of_minus_one(IntegerType power);

Computes $(-1)^{power}$ efficiently.

FloatType math::deg_to_rad(FloatType degrees);
FloatType math::rad_to_deg(FloatType radians);

Converts degrees to radians and back.

Meshing

struct Points {
    explicit Points(std::size_t count);
};

struct Intervals {
    explicit Intervals(std::size_t count);
    Intervals(Points points);
};

“Strong typedefs” for grid size in 1D meshing operations, which allows caller to express their intent more clearly.

Points implicitly converts to Intervals ($N + 1$ points corresponds to $N$ intervals), allowing most meshing functions to accept both types as an input.

std::vector<FloatType> math::linspace(FloatType L1, FloatType L2, Intervals N);

Meshes $[L_1, L_2]$ range into a regular 1D grid with $N$ intervals (which corresponds to $N + 1$ grid points). Similar to linspace() from Matlab and numpy.

FloatType math::integrate_trapezoidal(FuncType f, FloatType L1, FloatType L2, Intervals N);

Numericaly computes integral $I_h = \int\limits_{L_1}^{L_2} f(x) \mathrm{d} x$ over $N$ integration intervals using trapezoidal rule.

Other utils

UintType math::uint_difference(UintType a, UintType b);
Returns $ a - b $ for unsigned types accounting for possible integer overflow. Useful when working with small types in image processing.
int math::ssize(const SizedContainer& container);

Returns .size() value of given container casted to int.

Useful to reduce verbosity of static_cast<int>(container.size()) when working with int indexation that gets compared against container size.

ArithmeticType math::ternary_branchless(bool condition, ArithmeticType return_if_true, ArithmeticType return_if_false);

Returns condition ? return_if_true : return_if_false rewritten in a branchless way. Useful when working with GPU’s. Requires type to be arithmetic.

IntType math::ternary_bitselect(bool condition, IntType return_if_true, IntType return_if_false);
IntType math::ternary_bitselect(bool condition, IntType return_if_true);

Faster implementation of ternary_branchless() for integers. When second return is omitted, assumption return_if_false == 0 allows additional optimizations to be performed.

template<typename T, MemoryUnit units = MemoryUnit::MiB>
constexpr double memory_size(std::size_t count);

Returns size in units occupied in memory by count elements of type T. Useful to estimate memory usage of arrays, matrices and other data structures in a human-readable way.

Examples

Using math type traits

[ Run this code ]

using namespace utl;

std::cout
    << std::boolalpha
    << "are doubles addable?    -> " << math::is_addable_with_itself<double>::value              << "\n"
    << "are std::pairs addable? -> " << math::is_addable_with_itself<std::pair<int, int>>::value << "\n";

Output:

are doubles addable?    -> true
are std::pairs addable? -> false

Using basic math functions

[ Run this code ]

using namespace utl;

std::cout
    << "All methods below are constexpr and type agnostic:\n"
    << "abs(-4) = "                               << math::abs(-4)                               << "\n"
    << "sign(-4) = "                              << math::sign(-4)                              << "\n"
    << "sqr(-4) = "                               << math::sqr(-4)                               << "\n"
    << "cube(-4) = "                              << math::cube(-4)                              << "\n"
    << "deg_to_rad(180.) = "                      << math::deg_to_rad(180.)                      << "\n"
    << "rad_to_deg(PI) = "                        << math::rad_to_deg(math::PI)                  << "\n"
    << "\n"
    << "uint_difference(5u, 17u) = "              << math::uint_difference(5u, 17u)              << "\n"
    << "\n"
    << "ternary_branchless(true, 3.12, -4.17) = " << math::ternary_branchless(true, 3.12, -4.17) << "\n"
    << "ternary_bitselect(true, 15, -5) = "       << math::ternary_bitselect(true, 15, -5)       << "\n"
    << "ternary_bitselect(false, 15) = "          << math::ternary_bitselect(false, 15)          << "\n";

Output:

All methods below are constexpr and type agnostic:
abs(-4) = 4
sign(-4) = -1
sqr(-4) = 16
cube(-4) = -64
deg_to_rad(180.) = 3.14159
rad_to_deg(PI) = 180

uint_difference(5u, 17u) = 12

ternary_branchless(true, 3.12, -4.17) = 3.12
ternary_bitselect(true, 15, -5) = 15
ternary_bitselect(false, 15) = 0

Meshing and integrating

[ Run this code ]

// Mesh interval [0, PI] into 100 equal intervals => 101 linearly spaced points 
auto grid_1 = math::linspace(0., math::PI, math::Intervals(100));
auto grid_2 = math::linspace(0., math::PI, math::Points(   101)); // same as above

// Get array memory size
std::cout << "'grid_1' occupies " << math::memory_size<double, math::MemoryUnit::KB>(grid_1.size()) << " KB in memory\n\n";

// Integrate a function over an interval
auto f = [](double x){ return 4. / (1. + std::tan(x)); };
double integral = math::integrate_trapezoidal(f, 0., math::PI_HALF, math::Intervals(200));
std::cout << "Integral evaluates to: " << integral << " (should be ~PI)\n";

Output:

'grid_1' occupies 0.808 KB in memory

Integral evaluates to: 3.14159 (should be ~PI)