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