aboutsummaryrefslogtreecommitdiffstats
path: root/util/generic/ymath.h
blob: 163e581b7a06b06847abe089d92f717f330b1733 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
#pragma once

#include <util/system/yassert.h>
#include <util/system/defaults.h>

#include <cmath>
#include <cfloat>
#include <cstdlib>

#include "typetraits.h"
#include "utility.h"

constexpr double PI = M_PI;
constexpr double M_LOG2_10 = 3.32192809488736234787; // log2(10)
constexpr double M_LN2_INV = M_LOG2E;                // 1 / ln(2) == log2(e)

/**
 * \returns                             Absolute value of the provided argument.
 */
template <class T>
constexpr T Abs(T value) {
    return std::abs(value);
}

/**
 * @returns                             Base 2 logarithm of the provided double
 *                                      precision floating point value.
 */
inline double Log2(double value) {
    return log(value) * M_LN2_INV;
}

/**
 * @returns                             Base 2 logarithm of the provided
 *                                      floating point value.
 */
inline float Log2(float value) {
    return logf(value) * static_cast<float>(M_LN2_INV);
}

/**
 * @returns                             Base 2 logarithm of the provided integral value.
 */
template <class T>
inline std::enable_if_t<std::is_integral<T>::value, double>
Log2(T value) {
    return Log2(static_cast<double>(value));
}

/** Returns 2^x */
double Exp2(double);
float Exp2f(float);

template <class T>
static constexpr T Sqr(const T t) noexcept {
    return t * t;
}

inline double Sigmoid(double x) {
    return 1.0 / (1.0 + std::exp(-x));
}

inline float Sigmoid(float x) {
    return 1.0f / (1.0f + std::exp(-x));
}

static inline bool IsFinite(double f) {
#if defined(isfinite)
    return isfinite(f);
#elif defined(_win_)
    return _finite(f) != 0;
#elif defined(_darwin_)
    return isfinite(f);
#elif defined(__GNUC__)
    return __builtin_isfinite(f);
#elif defined(_STLP_VENDOR_STD)
    return _STLP_VENDOR_STD::isfinite(f);
#else
    return std::isfinite(f);
#endif
}

static inline bool IsNan(double f) {
#if defined(_win_)
    return _isnan(f) != 0;
#else
    return std::isnan(f);
#endif
}

inline bool IsValidFloat(double f) {
    return IsFinite(f) && !IsNan(f);
}

#ifdef _MSC_VER
double Erf(double x);
#else
inline double Erf(double x) {
    return erf(x);
}
#endif

/**
 * @returns                             Natural logarithm of the absolute value
 *                                      of the gamma function of provided argument.
 */
inline double LogGamma(double x) noexcept {
#if defined(_glibc_)
    int sign;

    (void)sign;

    return lgamma_r(x, &sign);
#elif defined(__GNUC__)
    return __builtin_lgamma(x);
#elif defined(_unix_)
    return lgamma(x);
#else
    extern double LogGammaImpl(double);
    return LogGammaImpl(x);
#endif
}

/**
 * @returns                             x^n for integer n, n >= 0.
 */
template <class T, class Int>
T Power(T x, Int n) {
    static_assert(std::is_integral<Int>::value, "only integer powers are supported");
    Y_ASSERT(n >= 0);
    if (n == 0) {
        return T(1);
    }
    while ((n & 1) == 0) {
        x = x * x;
        n >>= 1;
    }
    T result = x;
    n >>= 1;
    while (n > 0) {
        x = x * x;
        if (n & 1) {
            result = result * x;
        }
        n >>= 1;
    }
    return result;
}

/**
 * Compares two floating point values and returns true if they are considered equal.
 * The two numbers are compared in a relative way, where the exactness is stronger
 * the smaller the numbers are.
 *
 * Note that comparing values where either one is 0.0 will not work.
 * The solution to this is to compare against values greater than or equal to 1.0.
 *
 * @code
 * // Instead of comparing with 0.0
 * FuzzyEquals(0.0, 1.0e-200); // This will return false
 * // Compare adding 1 to both values will fix the problem
 * FuzzyEquals(1 + 0.0, 1 + 1.0e-200); // This will return true
 * @endcode
 */
inline bool FuzzyEquals(double p1, double p2, double eps = 1.0e-13) {
    return (Abs(p1 - p2) <= eps * Min(Abs(p1), Abs(p2)));
}

/**
 * @see FuzzyEquals(double, double, double)
 */
inline bool FuzzyEquals(float p1, float p2, float eps = 1.0e-6) {
    return (Abs(p1 - p2) <= eps * Min(Abs(p1), Abs(p2)));
}

namespace NUtilMathPrivate {
    template <bool IsSigned>
    struct TCeilDivImpl {};

    template <>
    struct TCeilDivImpl<true> {
        template <class T>
        static inline constexpr T Do(T x, T y) noexcept {
            return x / y + (((x < 0) ^ (y > 0)) && (x % y));
        }
    };

    template <>
    struct TCeilDivImpl<false> {
        template <class T>
        static inline constexpr T Do(T x, T y) noexcept {
            auto quot = x / y;
            return (x % y) ? (quot + 1) : quot;
        }
    };
} // namespace NUtilMathPrivate

/**
 * @returns Equivalent to ceil((double) x / (double) y) but using only integer arithmetic operations
 */
template <class T>
inline T
#if !defined(__NVCC__)
    constexpr
#endif
    CeilDiv(T x, T y) noexcept {
    static_assert(std::is_integral<T>::value, "Integral type required.");
    Y_ASSERT(y != 0);
    return ::NUtilMathPrivate::TCeilDivImpl<std::is_signed<T>::value>::Do(x, y);
}