lmms_math.h

Go to the documentation of this file.

/*
 * lmms_math.h - defines math functions
 *
 * Copyright (c) 2004-2008 Tobias Doerffel <tobydox/at/users.sourceforge.net>
 * 
 * This file is part of LMMS - https://lmms.io
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program (see COPYING); if not, write to the
 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 * Boston, MA 02110-1301 USA.
 *
 */
 
#ifndef LMMS_MATH_H
#define LMMS_MATH_H
 
#include <QtGlobal>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <numbers>
#include <concepts>
 
#include "lmms_constants.h"
 
#ifdef __SSE2__
#include <emmintrin.h>
#endif
 
namespace lmms
{
 
 
// TODO C++23: Make constexpr since std::abs() will be constexpr
inline bool approximatelyEqual(float x, float y) noexcept
{
        return x == y || std::abs(x - y) < F_EPSILON;
}
 
// TODO C++23: Make constexpr since std::trunc() will be constexpr
inline auto fraction(std::floating_point auto x) noexcept
{
        return x - std::trunc(x);
}
 
 
// TODO C++23: Make constexpr since std::floor() will be constexpr
inline auto absFraction(std::floating_point auto x) noexcept
{
        return x - std::floor(x);
}
 

inline int fastRand() noexcept
{
        thread_local unsigned long s_next = 1;
        s_next = s_next * 1103515245 + 12345;
        return s_next / 65536 % 32768;
}
 

template<std::floating_point T>
inline T fastRand(T upper) noexcept
{
        constexpr auto FAST_RAND_RATIO = static_cast<T>(1.0 / 32768);
        return fastRand() * upper * FAST_RAND_RATIO;
}
 

template<std::integral T>
inline T fastRand(T upper) noexcept
{
        constexpr float FAST_RAND_RATIO = 1.f / 32768;
        return static_cast<T>(fastRand() * static_cast<float>(upper) * FAST_RAND_RATIO);
}
 

template<typename T> requires std::is_arithmetic_v<T>
inline auto fastRand(T from, T to) noexcept { return from + fastRand(to - from); }
 

template<std::floating_point T>
inline T fastRandInc(T upper) noexcept
{
        constexpr auto FAST_RAND_RATIO = static_cast<T>(1.0 / 32767);
        return fastRand() * upper * FAST_RAND_RATIO;
}
 

template<std::unsigned_integral T>
inline T fastRandInc(T upper) noexcept
{
        // The integer specialization of this function is kind of weird, but it is
        // necessary to prevent massive bias away from the maximum value.
        // FAST_RAND_RATIO here is 1 greater than normal, so when multiplied by, it
        // will result in a random float within [0, @p upper) instead of the usual
        // [0, @p upper].
        constexpr float FAST_RAND_RATIO = 1.f / 32768;
        // Since the random float will be in a range that does not include @upper
        // due to the above ratio, increase the upper bound by 1.
        // All values greater than @p upper get rounded down to @p upper, making the
        // chance of returning a value of @p upper the same as any other of the
        // possible values.
        // No need to copysign() unlike the signed_integral overload, since it will always be positive
        return static_cast<T>(fastRand() * (upper + 1.f) * FAST_RAND_RATIO);
}
 

template<std::signed_integral T>
inline T fastRandInc(T upper) noexcept
{
        // The integer specialization of this function is kind of weird, but it is
        // necessary to prevent massive bias away from the maximum value.
        // FAST_RAND_RATIO here is 1 greater than normal, so when multiplied by, it
        // will result in a random float within [0, @p upper) instead of the usual
        // [0, @p upper].
        constexpr float FAST_RAND_RATIO = 1.f / 32768;
        // Since the random float will be in a range that does not include @upper
        // due to the above ratio, increase the magnitude of the upper bound by 1.
        // All values greater than @p upper get rounded down to @p upper, making the
        // chance of returning a value of @p upper the same as any other of the
        // possible values.
        // HACK: Even on -O3, without this static_cast, it will convert @p upper to float twice for some reason
        const auto fupper = static_cast<float>(upper);
        const float r = fupper + std::copysign(1.f, fupper);
        // Always round towards 0 (implicit truncation occurs during static_cast).
        return static_cast<T>(fastRand() * r * FAST_RAND_RATIO);
}
 

template<typename T> requires std::is_arithmetic_v<T>
inline auto fastRandInc(T from, T to) noexcept { return from + fastRandInc(to - from); }
 

inline bool oneIn(unsigned chance) noexcept { return 0 == (fastRand() % chance); }
 

template<class T>
static void roundAt(T& value, const T& where, const T& stepSize)
{
        if (std::abs(value - where) < F_EPSILON * std::abs(stepSize))
        {
                value = where;
        }
}

inline double fastPow(double a, double b)
{
        double d;
        std::int32_t x[2];
 
        std::memcpy(x, &a, sizeof(x));
        x[1] = static_cast<std::int32_t>(b * (x[1] - 1072632447) + 1072632447);
        x[0] = 0;
 
        std::memcpy(&d, x, sizeof(d));
        return d;
}
 

template<typename T>
constexpr T sign(T val) noexcept
{ 
        return val >= 0 ? 1 : -1; 
}
 

inline float sqrt_neg(float val) 
{
        return std::sqrt(std::abs(val)) * sign(val);
}

inline float signedPowf(float v, float e)
{
        return std::pow(std::abs(v), e) * sign(v);
}
 

inline float logToLinearScale(float min, float max, float value)
{
        using namespace std::numbers;
        if (min < 0)
        {
                const float mmax = std::max(std::abs(min), std::abs(max));
                const float val = value * (max - min) + min;
                float result = signedPowf(val / mmax, e_v<float>) * mmax;
                return std::isnan(result) ? 0 : result;
        }
        float result = std::pow(value, e_v<float>) * (max - min) + min;
        return std::isnan(result) ? 0 : result;
}
 

inline float linearToLogScale(float min, float max, float value)
{
        constexpr auto inv_e = static_cast<float>(1.0 / std::numbers::e);
        const float valueLimited = std::clamp(value, min, max);
        const float val = (valueLimited - min) / (max - min);
        if (min < 0)
        {
                const float mmax = std::max(std::abs(min), std::abs(max));
                float result = signedPowf(valueLimited / mmax, inv_e) * mmax;
                return std::isnan(result) ? 0 : result;
        }
        float result = std::pow(val, inv_e) * (max - min) + min;
        return std::isnan(result) ? 0 : result;
}
 
 
// TODO C++26: Make constexpr since std::exp() will be constexpr
template<typename T> requires std::is_arithmetic_v<T>
inline auto fastPow10f(T x)
{
        using F_T = std::conditional_t<std::is_floating_point_v<T>, T, float>;
        return std::exp(std::numbers::ln10_v<F_T> * x);
}
 
 
// TODO C++26: Make constexpr since std::log() will be constexpr
template<typename T> requires std::is_arithmetic_v<T>
inline auto fastLog10f(T x)
{
        using F_T = std::conditional_t<std::is_floating_point_v<T>, T, float>;
        constexpr auto inv_ln10 = static_cast<F_T>(1.0 / std::numbers::ln10);
        return std::log(x) * inv_ln10;
}
 

inline float ampToDbfs(float amp)
{
        return fastLog10f(amp) * 20.0f;
}
 

inline float dbfsToAmp(float dbfs)
{
        return fastPow10f(dbfs * 0.05f);
}
 

inline float safeAmpToDbfs(float amp)
{
        return amp == 0.0f ? -INFINITY : ampToDbfs(amp);
}
 

inline float safeDbfsToAmp(float dbfs)
{
        return std::isinf(dbfs) ? 0.0f : dbfsToAmp(dbfs);
}
 
 
// TODO C++20: use std::formatted_size
// @brief Calculate number of digits which LcdSpinBox would show for a given number
inline int numDigitsAsInt(float f)
{
        // use rounding:
        // LcdSpinBox sometimes uses std::round(), sometimes cast rounding
        // we use rounding to be on the "safe side"
        int asInt = static_cast<int>(std::round(f));
        int digits = 1; // always at least 1
        if(asInt < 0)
        {
                ++digits;
                asInt = -asInt;
        }
        // "asInt" is positive from now
        int power = 1;
        for (int i = 1; i < 10; ++i)
        {
                power *= 10;
                if (asInt >= power) { ++digits; } // 2 digits for >=10, 3 for >=100
                else { break; }
        }
        return digits;
}
 
template <typename T>
class LinearMap
{
public:
        LinearMap(T x1, T y1, T x2, T y2)
        {
                T const dx = x2 - x1;
                assert (dx != T(0));
 
                m_a = (y2 - y1) / dx;
                m_b = y1 - m_a * x1;
        }
 
        T map(T x) const
        {
                return m_a * x + m_b;
        }
 
private:
        T m_a;
        T m_b;
};
 
#ifdef __SSE2__
// exp approximation for SSE2: https://stackoverflow.com/a/47025627/5759631
// Maximum relative error of 1.72863156e-3 on [-87.33654, 88.72283]
inline __m128 fastExp(__m128 x)
{
        __m128 f, p, r;
        __m128i t, j;
        const __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
        const __m128i m = _mm_set1_epi32 (0xff800000); /* mask for integer bits */
        const __m128 ttm23 = _mm_set1_ps (1.1920929e-7f); /* exp2(-23) */
        const __m128 c0 = _mm_set1_ps (0.3371894346f);
        const __m128 c1 = _mm_set1_ps (0.657636276f);
        const __m128 c2 = _mm_set1_ps (1.00172476f);
 
        t = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
        j = _mm_and_si128 (t, m);            /* j = (int)(floor (x/log(2))) << 23 */
        t = _mm_sub_epi32 (t, j);
        f = _mm_mul_ps (ttm23, _mm_cvtepi32_ps (t)); /* f = (x/log(2)) - floor (x/log(2)) */
        p = c0;                              /* c0 */
        p = _mm_mul_ps (p, f);               /* c0 * f */
        p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
        p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
        p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= 2^f */
        r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
        return r;
}
 
// Lost Robot's SSE2 adaptation of Kari's vectorized log approximation: https://stackoverflow.com/a/65537754/5759631
// Maximum relative error of 7.922410e-4 on [1.0279774e-38f, 3.4028235e+38f]
inline __m128 fastLog(__m128 a)
{
        __m128i aInt = _mm_castps_si128(a);
        __m128i e = _mm_sub_epi32(aInt, _mm_set1_epi32(0x3f2aaaab));
        e = _mm_and_si128(e, _mm_set1_epi32(0xff800000));
        __m128i subtr = _mm_sub_epi32(aInt, e);
        __m128 m = _mm_castsi128_ps(subtr);
        __m128 i = _mm_mul_ps(_mm_cvtepi32_ps(e), _mm_set1_ps(1.19209290e-7f));
        __m128 f = _mm_sub_ps(m, _mm_set1_ps(1.0f));
        __m128 s = _mm_mul_ps(f, f);
        __m128 r = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(0.230836749f), f), _mm_set1_ps(-0.279208571f));
        __m128 t = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(0.331826031f), f), _mm_set1_ps(-0.498910338f));
        r = _mm_add_ps(_mm_mul_ps(r, s), t);
        r = _mm_add_ps(_mm_mul_ps(r, s), f);
        r = _mm_add_ps(_mm_mul_ps(i, _mm_set1_ps(0.693147182f)), r);
        return r;
}
 
inline __m128 sse2Abs(__m128 x)
{
        return _mm_and_ps(x, _mm_castsi128_ps(_mm_set1_epi32(0x7fffffff)));// clear sign bit
}
 
inline __m128 sse2Floor(__m128 x)
{
        __m128 t = _mm_cvtepi32_ps(_mm_cvttps_epi32(x)); // trunc toward 0
        __m128 needs_correction = _mm_cmplt_ps(x, t); // checks if x < trunc
        return _mm_sub_ps(t, _mm_and_ps(needs_correction, _mm_set1_ps(1.0f)));
}
 
inline __m128 sse2Round(__m128 x)
{
        __m128 sign_mask = _mm_cmplt_ps(x, _mm_setzero_ps());// checks if x < 0
        __m128 bias_pos = _mm_set1_ps(0.5f);
        __m128 bias_neg = _mm_set1_ps(-0.5f);
        __m128 bias = _mm_or_ps(_mm_and_ps(sign_mask, bias_neg), _mm_andnot_ps(sign_mask, bias_pos));
        __m128 y = _mm_add_ps(x, bias);
        return _mm_cvtepi32_ps(_mm_cvttps_epi32(y));
}
 
#endif // __SSE2__
 
} // namespace lmms
 
#endif // LMMS_MATH_H