276 lines
7.9 KiB
C++
276 lines
7.9 KiB
C++
#include <BAN/Math.h>
|
|
|
|
#include <math.h>
|
|
#include <stddef.h>
|
|
|
|
#define BAN_FUNC1(func) \
|
|
float func##f(float a) { return BAN::Math::func<float>(a); } \
|
|
double func(double a) { return BAN::Math::func<double>(a); } \
|
|
long double func##l(long double a) { return BAN::Math::func<long double>(a); }
|
|
|
|
#define BAN_FUNC2(func) \
|
|
float func##f(float a, float b) { return BAN::Math::func<float>(a, b); } \
|
|
double func(double a, double b) { return BAN::Math::func<double>(a, b); } \
|
|
long double func##l(long double a, long double b) { return BAN::Math::func<long double>(a, b); }
|
|
|
|
#define BAN_FUNC2_PTR(func) \
|
|
float func##f(float a, float* b) { return BAN::Math::func<float>(a, b); } \
|
|
double func(double a, double* b) { return BAN::Math::func<double>(a, b); } \
|
|
long double func##l(long double a, long double* b) { return BAN::Math::func<long double>(a, b); }
|
|
|
|
#define BAN_FUNC2_TYPE(func, type) \
|
|
float func##f(float a, type b) { return BAN::Math::func<float>(a, b); } \
|
|
double func(double a, type b) { return BAN::Math::func<double>(a, b); } \
|
|
long double func##l(long double a, type b) { return BAN::Math::func<long double>(a, b); }
|
|
|
|
#define FUNC_EXPR1(func, ...) \
|
|
float func##f(float a) { return __VA_ARGS__; } \
|
|
double func(double a) { return __VA_ARGS__; } \
|
|
long double func##l(long double a) { return __VA_ARGS__; }
|
|
|
|
#define FUNC_EXPR2(func, ...) \
|
|
float func##f(float a, float b) { return __VA_ARGS__; } \
|
|
double func(double a, double b) { return __VA_ARGS__; } \
|
|
long double func##l(long double a, long double b) { return __VA_ARGS__; }
|
|
|
|
#define FUNC_EXPR3(func, ...) \
|
|
float func##f(float a, float b, float c) { return __VA_ARGS__; } \
|
|
double func(double a, double b, double c) { return __VA_ARGS__; } \
|
|
long double func##l(long double a, long double b, long double c) { return __VA_ARGS__; }
|
|
|
|
#define FUNC_EXPR1_RET(func, ret, ...) \
|
|
ret func##f(float a) { return __VA_ARGS__; } \
|
|
ret func(double a) { return __VA_ARGS__; } \
|
|
ret func##l(long double a) { return __VA_ARGS__; }
|
|
|
|
#define FUNC_EXPR2_TYPE(func, type, ...) \
|
|
float func##f(float a, type b) { return __VA_ARGS__; } \
|
|
double func(double a, type b) { return __VA_ARGS__; } \
|
|
long double func##l(long double a, type b) { return __VA_ARGS__; }
|
|
|
|
template<BAN::floating_point T> struct float_underlying;
|
|
template<> struct float_underlying<float> { using type = uint32_t; };
|
|
template<> struct float_underlying<double> { using type = uint64_t; };
|
|
template<> struct float_underlying<long double> { using type = uint64_t; };
|
|
|
|
template<BAN::integral T, size_t mantissa_bits, size_t exponent_bits, bool integral>
|
|
struct __FloatDecompose;
|
|
|
|
template<BAN::integral T, size_t mantissa_bits, size_t exponent_bits>
|
|
struct __FloatDecompose<T, mantissa_bits, exponent_bits, true>
|
|
{
|
|
T mantissa : mantissa_bits;
|
|
T : 1;
|
|
T exponent : exponent_bits;
|
|
T sign : 1;
|
|
};
|
|
|
|
template<BAN::integral T, size_t mantissa_bits, size_t exponent_bits>
|
|
struct __FloatDecompose<T, mantissa_bits, exponent_bits, false>
|
|
{
|
|
T mantissa : mantissa_bits;
|
|
T exponent : exponent_bits;
|
|
T sign : 1;
|
|
};
|
|
|
|
template<BAN::floating_point T>
|
|
struct FloatDecompose
|
|
{
|
|
using value_type = float_underlying<T>::type;
|
|
|
|
static constexpr size_t mantissa_bits
|
|
= BAN::is_same_v<T, float> ? FLT_MANT_DIG - 1
|
|
: BAN::is_same_v<T, double> ? DBL_MANT_DIG - 1
|
|
: BAN::is_same_v<T, long double> ? LDBL_MANT_DIG - 1
|
|
: 0;
|
|
|
|
static constexpr size_t exponent_bits
|
|
= BAN::is_same_v<T, float> ? 8
|
|
: BAN::is_same_v<T, double> ? 11
|
|
: BAN::is_same_v<T, long double> ? 15
|
|
: 0;
|
|
|
|
static constexpr value_type mantissa_max
|
|
= (static_cast<value_type>(1) << mantissa_bits) - 1;
|
|
|
|
static constexpr value_type exponent_max
|
|
= (static_cast<value_type>(1) << exponent_bits) - 1;
|
|
|
|
FloatDecompose(T x)
|
|
: raw(x)
|
|
{}
|
|
|
|
union
|
|
{
|
|
T raw;
|
|
__FloatDecompose<
|
|
value_type,
|
|
mantissa_bits,
|
|
exponent_bits,
|
|
BAN::is_same_v<T, long double>
|
|
> bits;
|
|
};
|
|
};
|
|
|
|
template<BAN::floating_point T1, BAN::floating_point T2>
|
|
static T1 nextafter_impl(T1 x, T2 y)
|
|
{
|
|
if (isnan(x) || isnan(y))
|
|
return BAN::numeric_limits<T1>::quiet_NaN();
|
|
|
|
if (!isfinite(x) || x == y)
|
|
return x;
|
|
|
|
FloatDecompose<T1> decompose(x);
|
|
|
|
// at zero
|
|
if ((T2)x == (T2)0.0)
|
|
{
|
|
decompose.bits.mantissa = 1;
|
|
decompose.bits.sign = (x > y);
|
|
return decompose.raw;
|
|
}
|
|
|
|
// away from zero
|
|
if ((x > y) == decompose.bits.sign)
|
|
{
|
|
decompose.bits.mantissa++;
|
|
if (decompose.bits.mantissa == 0)
|
|
{
|
|
decompose.bits.exponent++;
|
|
if (decompose.bits.exponent == decompose.exponent_max)
|
|
decompose.bits.mantissa = 0;
|
|
}
|
|
return decompose.raw;
|
|
}
|
|
|
|
// towards zero
|
|
decompose.bits.mantissa--;
|
|
if (decompose.bits.mantissa == decompose.mantissa_max)
|
|
{
|
|
if (decompose.bits.exponent == 0)
|
|
return 0.0;
|
|
decompose.bits.exponent--;
|
|
}
|
|
|
|
return decompose.raw;
|
|
}
|
|
|
|
static long double tgamma_impl(long double x)
|
|
{
|
|
constexpr long double pi = BAN::numbers::pi_v<long double>;
|
|
|
|
if (x == 0.0L)
|
|
return BAN::numeric_limits<long double>::infinity();
|
|
|
|
// reflection formula
|
|
if (x < 0.5L)
|
|
return pi / (BAN::Math::sin(pi * x) * tgamma_impl(1.0L - x));
|
|
x -= 1.0L;
|
|
|
|
// Lanczos approximation
|
|
|
|
constexpr long double g = 8.0L;
|
|
constexpr long double p[] {
|
|
0.9999999999999999298e+0L,
|
|
1.9753739023578852322e+3L,
|
|
-4.3973823927922428918e+3L,
|
|
3.4626328459862717019e+3L,
|
|
-1.1569851431631167820e+3L,
|
|
1.5453815050252775060e+2L,
|
|
-6.2536716123689161798e+0L,
|
|
3.4642762454736807441e-2L,
|
|
-7.4776171974442977377e-7L,
|
|
6.3041253821852264261e-8L,
|
|
-2.7405717035683877489e-8L,
|
|
4.0486948817567609101e-9L
|
|
};
|
|
constexpr long double sqrt_2pi = 2.5066282746310005024L;
|
|
|
|
long double A = p[0];
|
|
for (size_t i = 1; i < sizeof(p) / sizeof(*p); i++)
|
|
A += p[i] / (x + i);
|
|
|
|
const long double t = x + g + 0.5L;
|
|
return sqrt_2pi * BAN::Math::pow(t, x + 0.5L) * BAN::Math::exp(-t) * A;
|
|
}
|
|
|
|
static long double erf_impl(long double x)
|
|
{
|
|
long double sum = 0.0L;
|
|
for (size_t n = 0; n < 100; n++)
|
|
sum += x / (2.0L * n + 1.0L) * BAN::Math::ldexp(-x * x, n) / tgamma_impl(n + 1.0L);
|
|
|
|
constexpr long double sqrt_pi = 1.77245385090551602729L;
|
|
return 2.0L / sqrt_pi * sum;
|
|
}
|
|
|
|
__BEGIN_DECLS
|
|
|
|
// FIXME: add handling for nan and infinity values
|
|
|
|
BAN_FUNC1(acos)
|
|
BAN_FUNC1(acosh)
|
|
BAN_FUNC1(asin)
|
|
BAN_FUNC1(asinh)
|
|
BAN_FUNC1(atan)
|
|
BAN_FUNC2(atan2)
|
|
BAN_FUNC1(atanh)
|
|
BAN_FUNC1(cbrt)
|
|
BAN_FUNC1(ceil)
|
|
BAN_FUNC2(copysign)
|
|
BAN_FUNC1(cos)
|
|
BAN_FUNC1(cosh)
|
|
FUNC_EXPR1(erf, erf_impl(a))
|
|
FUNC_EXPR1(erfc, 1.0 - erf_impl(a))
|
|
BAN_FUNC1(exp)
|
|
BAN_FUNC1(exp2)
|
|
FUNC_EXPR1(expm1, ({ a -= 1.0; BAN::Math::exp(a); }))
|
|
FUNC_EXPR1(fabs, BAN::Math::abs(a))
|
|
FUNC_EXPR2(fdim, a > b ? a - b : 0.0)
|
|
BAN_FUNC1(floor)
|
|
FUNC_EXPR3(fma, (a * b) + c)
|
|
FUNC_EXPR2(fmax, a < b ? b : a)
|
|
FUNC_EXPR2(fmin, a < b ? a : b)
|
|
BAN_FUNC2(fmod)
|
|
BAN_FUNC2_TYPE(frexp, int*)
|
|
BAN_FUNC2(hypot)
|
|
FUNC_EXPR1_RET(ilogb, int, BAN::Math::logb(a))
|
|
// j0, j1, jn
|
|
BAN_FUNC2_TYPE(ldexp, int)
|
|
FUNC_EXPR1(lgamma, BAN::Math::log(BAN::Math::abs(tgamma_impl(a))))
|
|
FUNC_EXPR1_RET(llrint, long long, BAN::Math::rint(a))
|
|
FUNC_EXPR1_RET(llround, long long, BAN::Math::round(a))
|
|
BAN_FUNC1(log)
|
|
BAN_FUNC1(log10)
|
|
FUNC_EXPR1(log1p, ({ a += 1.0; BAN::Math::log(a); }));
|
|
BAN_FUNC1(log2)
|
|
BAN_FUNC1(logb)
|
|
FUNC_EXPR1_RET(lrint, long, BAN::Math::rint(a))
|
|
FUNC_EXPR1_RET(lround, long, BAN::Math::round(a))
|
|
BAN_FUNC2_PTR(modf)
|
|
double nan(const char* tagp) { return __builtin_nan(tagp); }
|
|
float nanf(const char* tagp) { return __builtin_nanf(tagp); }
|
|
long double nanl(const char* tagp) { return __builtin_nanl(tagp); }
|
|
FUNC_EXPR1(nearbyint, BAN::Math::rint(a))
|
|
FUNC_EXPR2(nextafter, nextafter_impl(a, b))
|
|
FUNC_EXPR2_TYPE(nexttoward, long double, nextafter_impl(a, b))
|
|
FUNC_EXPR2(pow, ({ if (isnan(a)) return a; if (isnan(b)) return b; BAN::Math::pow(a, b); }))
|
|
// remainder
|
|
// remquo
|
|
BAN_FUNC1(rint)
|
|
FUNC_EXPR1(round, ({ if (!isfinite(a)) return a; BAN::Math::round(a); }))
|
|
FUNC_EXPR2_TYPE(scalbln, long, BAN::Math::scalbn(a, b))
|
|
FUNC_EXPR2_TYPE(scalbn, int, BAN::Math::scalbn(a, b))
|
|
BAN_FUNC1(sin)
|
|
BAN_FUNC1(sinh)
|
|
BAN_FUNC1(sqrt)
|
|
BAN_FUNC1(tan)
|
|
BAN_FUNC1(tanh)
|
|
FUNC_EXPR1(tgamma, tgamma_impl(a))
|
|
BAN_FUNC1(trunc)
|
|
// y0, y1, yn
|
|
|
|
|
|
__END_DECLS
|