tatami
C++ API for different matrix representations
Loading...
Searching...
No Matches
arithmetic_utils.hpp
Go to the documentation of this file.
1#ifndef TATAMI_ARITHMETIC_UTILS_HPP
2#define TATAMI_ARITHMETIC_UTILS_HPP
3
4#include <cmath>
5
12namespace tatami {
13
23enum class ArithmeticOperation : char {
24 ADD,
25 SUBTRACT,
26 MULTIPLY,
27 DIVIDE,
28 POWER,
29 MODULO,
30 INTEGER_DIVIDE
31};
32
36// We deliberately use an auto type so as to defer a decision on what the output
37// type should be; an appropriate coercion is left to the caller classes.
38template<ArithmeticOperation op_, bool right_, typename Value_, typename Scalar_>
39auto delayed_arithmetic(Value_ val, Scalar_ scalar) {
40 auto left = [&]{
41 if constexpr(right_) {
42 return val;
43 } else {
44 return scalar;
45 }
46 }();
47
48 auto right = [&]{
49 if constexpr(right_) {
50 return scalar;
51 } else {
52 return val;
53 }
54 }();
55
56 if constexpr(op_ == ArithmeticOperation::ADD) {
57 return left + right;
58
59 } else if constexpr(op_ == ArithmeticOperation::MULTIPLY) {
60 return left * right;
61
62 } else if constexpr(op_ == ArithmeticOperation::SUBTRACT) {
63 return left - right;
64
65 } else if constexpr(op_ == ArithmeticOperation::DIVIDE) {
66 // Assume that either Value_ is an IEEE-754 float, or that division by
67 // zero is impossible in this context. We don't apply manual checks
68 // here to avoid performance degradation; we also don't check that the
69 // other operations yield a value that doesn't overflow/underflow, so
70 // it would be odd to make an exception for div-by-zero errors.
71 return left / right;
72
73 } else if constexpr(op_ == ArithmeticOperation::POWER) {
74 return std::pow(left, right);
75
76 } else if constexpr(op_ == ArithmeticOperation::MODULO) {
77 // Based on a floored divide, so some work is necessary to
78 // get the right value when the sign is negative.
79 auto quo = left / right;
80 if constexpr(std::numeric_limits<decltype(quo)>::is_integer) {
81 auto rem = left % right;
82 return rem + (quo < 0 && rem != 0 ? right : 0);
83 } else {
84 auto rem = std::fmod(left, right);
85 return rem + (quo < 0 && rem != 0 ? right : 0);
86 }
87
88 } else if constexpr(op_ == ArithmeticOperation::INTEGER_DIVIDE) {
89 auto out = left / right;
90 if constexpr(std::numeric_limits<decltype(out)>::is_integer) {
91 // Using a floored divide. This little branch should be optimized
92 // away so don't worry too much about it.
93 return out - (out < 0 ? (left % right != 0) : 0);
94 } else {
95 return std::floor(out);
96 }
97 }
98}
99
100// Some of the helpers need to divide by zero to compute the fill value. The
101// only way to guarantee that division by zero is supported is with IEEE
102// floats, otherwise it would be undefined behavior at compile time, and the
103// compiler could do anything, including refusing to compile it. So, we hide
104// any '0/0' behind a constexpr to ensure that the compiler doesn't see it.
105template<ArithmeticOperation op_, bool right_, typename Value_, typename Scalar_>
106constexpr bool has_unsafe_divide_by_zero() {
107 typedef decltype(delayed_arithmetic<op_, right_>(std::declval<Value_>(), std::declval<Scalar_>())) Product;
108 if constexpr(std::numeric_limits<Product>::is_iec559) {
109 return false;
110 }
111
112 // POWER also involves division by zero for negative powers,
113 // but it returns an implementation-defined value; so it's "safe".
114 // - https://en.cppreference.com/w/cpp/numeric/math/pow
115
116 if constexpr(op_ == ArithmeticOperation::DIVIDE) {
117 return true;
118 }
119 if constexpr(op_ == ArithmeticOperation::MODULO) {
120 return true;
121 }
122 if constexpr(op_ == ArithmeticOperation::INTEGER_DIVIDE) {
123 return true;
124 }
125
126 return false;
127}
128
129// COMMENT: if the fill value is visible at a compile time, the coercion to the
130// output type is also visible. If the output type cannot hold the fill value,
131// we could potentially get more compile-time UB. (This mostly concerns NaNs or
132// Infs from divide-by-zero, but could apply to overflows or out-of-range casts
133// to integer types.) To get around this, we assume that the compiler supports
134// the IEC 559 spec, where - according to Annex F.4 of the C standard - the
135// cast from float to integer is merely unspecified, not undefined.
136
140}
141
142#endif
Flexible representations for matrix data.
Definition Extractor.hpp:15
ArithmeticOperation
Definition arithmetic_utils.hpp:23