305 lines
11 KiB
C++
305 lines
11 KiB
C++
//! \addtogroup Math
|
||
//! \{
|
||
//! \file pimathbase.h
|
||
//! \brief
|
||
//! \~english Basic mathematical functions and defines
|
||
//! \~russian Базовые математические функции и дефайны
|
||
//! \details
|
||
//! \~english Common mathematical constants, conversion functions and utility functions
|
||
//! \~russian Общие математические константы, функции преобразования и утилиты
|
||
//! \}
|
||
/*
|
||
PIP - Platform Independent Primitives
|
||
Basic mathematical functions and defines
|
||
Ivan Pelipenko peri4ko@yandex.ru, Andrey Bychkov work.a.b@yandex.ru
|
||
|
||
This program is free software: you can redistribute it and/or modify
|
||
it under the terms of the GNU Lesser General Public License as published by
|
||
the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
|
||
|
||
You should have received a copy of the GNU Lesser General Public License
|
||
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||
*/
|
||
|
||
#ifndef PIMATHBASE_H
|
||
#define PIMATHBASE_H
|
||
|
||
#include "piinit.h"
|
||
#include "pipair.h"
|
||
#include "pivector.h"
|
||
#ifdef QNX
|
||
# undef PIP_MATH_J0
|
||
# undef PIP_MATH_J1
|
||
# undef PIP_MATH_JN
|
||
# undef PIP_MATH_Y0
|
||
# undef PIP_MATH_Y1
|
||
# undef PIP_MATH_YN
|
||
# include <math.h>
|
||
#else
|
||
# include <cmath>
|
||
#endif
|
||
|
||
#ifndef M_LN2
|
||
# define M_LN2 0.69314718055994530942
|
||
#endif
|
||
#ifndef M_LN10
|
||
# define M_LN10 2.30258509299404568402
|
||
#endif
|
||
#ifndef M_SQRT2
|
||
# define M_SQRT2 1.41421356237309514547
|
||
#endif
|
||
#ifndef M_SQRT3
|
||
# define M_SQRT3 1.73205080756887719318
|
||
#endif
|
||
#ifndef M_1_SQRT2
|
||
# define M_1_SQRT2 0.70710678118654746172
|
||
#endif
|
||
#ifndef M_1_SQRT3
|
||
# define M_1_SQRT3 0.57735026918962584208
|
||
#endif
|
||
#ifndef M_PI
|
||
# define M_PI 3.141592653589793238462643383280
|
||
#endif
|
||
#ifndef M_2PI
|
||
# define M_2PI 6.283185307179586476925286766559
|
||
#endif
|
||
#ifndef M_PI_3
|
||
# define M_PI_3 1.04719755119659774615
|
||
#endif
|
||
#ifndef M_2PI_3
|
||
# define M_2PI_3 2.0943951023931954923
|
||
#endif
|
||
#ifndef M_180_PI
|
||
# define M_180_PI 57.2957795130823208768
|
||
#endif
|
||
#ifndef M_PI_180
|
||
# define M_PI_180 1.74532925199432957692e-2
|
||
#endif
|
||
#ifndef M_SQRT_PI
|
||
# define M_SQRT_PI 1.772453850905516027298167483341
|
||
#endif
|
||
#ifndef M_E
|
||
# define M_E 2.7182818284590452353602874713527
|
||
#endif
|
||
#ifndef M_LIGHT_SPEED
|
||
# define M_LIGHT_SPEED 2.99792458e+8
|
||
#endif
|
||
#ifndef M_RELATIVE_CONST
|
||
# define M_RELATIVE_CONST -4.442807633e-10;
|
||
#endif
|
||
#ifndef M_GRAVITY_CONST
|
||
# define M_GRAVITY_CONST 398600.4418e9;
|
||
#endif
|
||
|
||
|
||
const double deg2rad = M_PI_180;
|
||
const double rad2deg = M_180_PI;
|
||
|
||
// clang-format off
|
||
inline int sign(const float & x) {return (x < 0.f) ? -1 : (x > 0.f ? 1 : 0);}
|
||
inline int sign(const double & x) {return (x < 0. ) ? -1 : (x > 0. ? 1 : 0);}
|
||
inline int sign(const ldouble & x) {return (x < 0.L) ? -1 : (x > 0.L ? 1 : 0);}
|
||
|
||
inline int pow2 (const int p ) {return 1 << p;}
|
||
inline float pow10(const float & e) {return powf(10.f, e);}
|
||
inline double pow10(const double & e) {return pow (10. , e);}
|
||
inline ldouble pow10(const ldouble & e) {return powl(10.L, e);}
|
||
// clang-format on
|
||
|
||
inline double sinc(const double & v) {
|
||
if (v == 0.) return 1.;
|
||
double t = M_PI * v;
|
||
return sin(t) / t;
|
||
}
|
||
|
||
//! Bessel function of the first kind of order 0
|
||
//! \~english Bessel function J0(x)
|
||
//! \~russian Функция Бесселя первого рода порядка 0
|
||
PIP_EXPORT double piJ0(const double & v);
|
||
//! Bessel function of the first kind of order 1
|
||
//! \~english Bessel function J1(x)
|
||
//! \~russian Функция Бесселя первого рода порядка 1
|
||
PIP_EXPORT double piJ1(const double & v);
|
||
//! Bessel function of the first kind of order n
|
||
//! \~english Bessel function Jn(n, x)
|
||
//! \~russian Функция Бесселя первого рода порядка n
|
||
PIP_EXPORT double piJn(int n, const double & v);
|
||
//! Bessel function of the second kind of order 0
|
||
//! \~english Bessel function Y0(x)
|
||
//! \~russian Функция Бесселя второго рода порядка 0
|
||
PIP_EXPORT double piY0(const double & v);
|
||
//! Bessel function of the second kind of order 1
|
||
//! \~english Bessel function Y1(x)
|
||
//! \~russian Функция Бесселя второго рода порядка 1
|
||
PIP_EXPORT double piY1(const double & v);
|
||
//! Bessel function of the second kind of order n
|
||
//! \~english Bessel function Yn(n, x)
|
||
//! \~russian Функция Бесселя второго рода порядка n
|
||
PIP_EXPORT double piYn(int n, const double & v);
|
||
|
||
// clang-format off
|
||
inline constexpr float toRad(float deg) {return deg * M_PI_180;}
|
||
inline constexpr double toRad(double deg) {return deg * M_PI_180;}
|
||
inline constexpr ldouble toRad(ldouble deg) {return deg * M_PI_180;}
|
||
inline constexpr float toDeg(float rad) {return rad * M_180_PI;}
|
||
inline constexpr double toDeg(double rad) {return rad * M_180_PI;}
|
||
inline constexpr ldouble toDeg(ldouble rad) {return rad * M_180_PI;}
|
||
// clang-format on
|
||
|
||
//! Square of a value
|
||
//! \~english Returns the square of value v (v * v)
|
||
//! \~russian Возвращает квадрат значения v (v * v)
|
||
template<typename T>
|
||
inline constexpr T sqr(const T & v) {
|
||
return v * v;
|
||
}
|
||
|
||
//! Convert linear value to decibels
|
||
//! \~english Convert linear value to decibels: 10 * log10(val)
|
||
//! \~russian Преобразовать линейное значение в децибелы: 10 * log10(val)
|
||
template<typename T>
|
||
inline constexpr T toDb(T val) {
|
||
return T(10.) * std::log10(val);
|
||
}
|
||
|
||
//! Convert decibels to linear value
|
||
//! \~english Convert decibels to linear value: 10^(val/10)
|
||
//! \~russian Преобразовать децибелы в линейное значение: 10^(val/10)
|
||
template<typename T>
|
||
inline constexpr T fromDb(T val) {
|
||
return std::pow(T(10.), val / T(10.));
|
||
}
|
||
|
||
// [-1 ; 1]
|
||
//! Generate random double in range [-1, 1]
|
||
//! \~english Returns random double in range [-1, 1]
|
||
//! \~russian Генерирует случайное число double в диапазоне [-1, 1]
|
||
PIP_EXPORT double randomd();
|
||
// [-1 ; 1] normal
|
||
//! Generate random double with normal (Gaussian) distribution
|
||
//! \~english Returns random double with normal distribution, mean=dv, stddev=sv
|
||
//! \~russian Генерирует случайное число double с нормальным распределением, среднее=dv, стандартное отклонение=sv
|
||
PIP_EXPORT double randomn(double dv = 0., double sv = 1.);
|
||
|
||
//! Absolute value of vector elements
|
||
//! \~english Returns vector with absolute values of each element
|
||
//! \~russian Возвращает вектор с абсолютными значениями каждого элемента
|
||
template<typename T>
|
||
inline PIVector<T> piAbs(const PIVector<T> & v) {
|
||
PIVector<T> result;
|
||
result.resize(v.size());
|
||
for (uint i = 0; i < v.size(); i++)
|
||
result[i] = piAbs<T>(v[i]);
|
||
return result;
|
||
}
|
||
|
||
|
||
//! Normalize angle to [0, 360) range (in-place)
|
||
//! \~english Normalizes angle to range [0, 360) degrees
|
||
//! \~russian Нормализует угол в диапазон [0, 360) градусов (на месте)
|
||
template<typename T>
|
||
void normalizeAngleDeg360(T & a) {
|
||
while (a < 0.)
|
||
a += 360.;
|
||
while (a > 360.)
|
||
a -= 360.;
|
||
}
|
||
//! Normalize angle to [0, 360) range
|
||
//! \~english Returns angle normalized to range [0, 360) degrees
|
||
//! \~russian Возвращает угол нормализованный в диапазон [0, 360) градусов
|
||
template<typename T>
|
||
double normalizedAngleDeg360(T a) {
|
||
normalizeAngleDeg360(a);
|
||
return a;
|
||
}
|
||
|
||
|
||
//! Normalize angle to (-180, 180] range (in-place)
|
||
//! \~english Normalizes angle to range (-180, 180] degrees
|
||
//! \~russian Нормализует угол в диапазон (-180, 180] градусов (на месте)
|
||
template<typename T>
|
||
void normalizeAngleDeg180(T & a) {
|
||
while (a < -180.)
|
||
a += 360.;
|
||
while (a > 180.)
|
||
a -= 360.;
|
||
}
|
||
//! Normalize angle to (-180, 180] range
|
||
//! \~english Returns angle normalized to range (-180, 180] degrees
|
||
//! \~russian Возвращает угол нормализованный в диапазон (-180, 180] градусов
|
||
template<typename T>
|
||
double normalizedAngleDeg180(T a) {
|
||
normalizeAngleDeg180(a);
|
||
return a;
|
||
}
|
||
|
||
|
||
//! Ordinary Least Squares linear regression
|
||
//! \~english Calculates linear regression coefficients using OLS method
|
||
//! \~russian Вычисляет коэффициенты линейной регрессии методом наименьших квадратов
|
||
//! \param input Vector of (x, y) pairs
|
||
//! \param out_a Output pointer for slope coefficient (a), can be nullptr
|
||
//! \param out_b Output pointer for intercept coefficient (b), can be nullptr
|
||
//! \return true on success
|
||
template<typename T>
|
||
bool OLS_Linear(const PIVector<PIPair<T, T>> & input, T * out_a, T * out_b) {
|
||
static_assert(std::is_arithmetic<T>::value, "Type must be arithmetic");
|
||
if (input.size_s() < 2) return false;
|
||
int n = input.size_s();
|
||
T a_t0 = T(), a_t1 = T(), a_t2 = T(), a_t3 = T(), a_t4 = T(), a = T(), b = T();
|
||
for (int i = 0; i < n; ++i) {
|
||
const PIPair<T, T> & cv(input[i]);
|
||
a_t0 += cv.first * cv.second;
|
||
a_t1 += cv.first;
|
||
a_t2 += cv.second;
|
||
a_t3 += cv.first * cv.first;
|
||
}
|
||
a_t4 = n * a_t3 - a_t1 * a_t1;
|
||
if (a_t4 != T()) a = (n * a_t0 - a_t1 * a_t2) / a_t4;
|
||
b = (a_t2 - a * a_t1) / n;
|
||
if (out_a != 0) *out_a = a;
|
||
if (out_b != 0) *out_b = b;
|
||
return true;
|
||
}
|
||
|
||
|
||
//! Weighted Least Squares linear regression
|
||
//! \~english Calculates linear regression coefficients using WLS method
|
||
//! \~russian Вычисляет коэффициенты линейной регрессии методом взвешенных наименьших квадратов
|
||
//! \param input Vector of (x, y) pairs
|
||
//! \param weights Vector of weights for each point
|
||
//! \param out_a Output pointer for slope coefficient (a), can be nullptr
|
||
//! \param out_b Output pointer for intercept coefficient (b), can be nullptr
|
||
//! \return true on success
|
||
template<typename T>
|
||
bool WLS_Linear(const PIVector<PIPair<T, T>> & input, const PIVector<T> & weights, T * out_a, T * out_b) {
|
||
static_assert(std::is_arithmetic<T>::value, "Type must be arithmetic");
|
||
if (input.size_s() < 2) return false;
|
||
if (input.size_s() != weights.size_s()) return false;
|
||
int n = input.size_s();
|
||
T a_t0 = T(), a_t1 = T(), a_t2 = T(), a_t3 = T(), a_t4 = T(), a_n = T(), a = T(), b = T();
|
||
for (int i = 0; i < n; ++i) {
|
||
T cp = weights[i];
|
||
const PIPair<T, T> & cv(input[i]);
|
||
a_t0 += cv.first * cv.second * cp;
|
||
a_t1 += cv.first * cp;
|
||
a_t2 += cv.second * cp;
|
||
a_t3 += cv.first * cv.first * cp;
|
||
a_n += cp;
|
||
}
|
||
a_t4 = a_n * a_t3 - a_t1 * a_t1;
|
||
if (a_t4 != T()) a = (a_n * a_t0 - a_t1 * a_t2) / a_t4;
|
||
b = (a_t2 - a * a_t1) / a_n;
|
||
if (out_a != 0) *out_a = a;
|
||
if (out_b != 0) *out_b = b;
|
||
return true;
|
||
}
|
||
|
||
#endif // PIMATHBASE_H
|