Commit a9401a8c authored by Elizabeth Myers's avatar Elizabeth Myers 💬
Browse files

math: implement most of glibc's __*_finite functions

These are supposed to be specialisations for speed, but these are just
faked. Some warnings were added too, if they return infinite values.

As a side effect of this change, scalbl is also now implemented.

As noted, not all functions are implemented; the big two blockers are an
implementation of j0l and y0l; I imagine Bessel functions aren't too
widely used, so I doubt that many things will want them. Someone (not
it) can implement them later.
parent db8ea145
#include <math.h> /* isfinite, isinf, isnan */
#define _GNU_SOURCE /* Extra maths functions */
#include <math.h> /* Literally everything */
#include "alias.h" /* weak_alias */
#include "alias.h" /* weak_alias */
#include "internal.h" /* GCOMPAT__assert_with_reason */
/**
* Multiplies the first argument x by FLT_RADIX (probably 2) to the power of y.
*/
long double scalbl(long double x, long double y)
{
/*
* XXX strictly not correct but:
* 1) Good Enough(TM)
* 2) scalbl is deprecated anyway
* */
return scalblnl(x, (long int)y);
}
/*
* The below require support for ynl/jnl which doesn't exist in musl and isn't
* implemented in gcompat yet
*/
#if 0
/**
* Return Bessel functions of x of the first kind of order n.
*/
long double jnl(int n, long double x)
{
/* TODO implement */
return 0;
}
/**
* Return Bessel functions of x of the first kind of order 0.
*/
long double j0l(long double n)
{
return jnl(0, n);
}
/**
* Return Bessel functions of x of the first kind of order 1.
*/
long double j1l(long double n)
{
return jnl(1, n);
}
/**
* Return Bessel functions of x of the second kind of order n.
*/
long double ynl(int n, long double x)
{
/* TODO implement */
return 0;
}
/**
* Return Bessel functions of x of the second kind of order 0.
*/
long double y0l(long double n)
{
return ynl(0, n);
}
/**
* Return Bessel functions of x of the second kind of order 1.
*/
long double y1l(long double n)
{
return ynl(1, n);
}
#endif
/**
* Test for finite value.
......@@ -73,6 +144,7 @@ int __isnan(double arg)
}
weak_alias(__isnan, isnan);
/**
* Test for a NaN.
*
......@@ -94,3 +166,1019 @@ int __isnanl(long double arg)
return isnan(arg);
}
weak_alias(__isnanl, isnanl);
/*
* Finite specialisations of functions used by glibc, that aren't supposed to
* return infinity.
*/
#define _ASSERT_FINITE(finite_fn, res) \
GCOMPAT__assert_with_reason(finite_fn(res), \
"infinite value returned in a function that returns a " \
"finite result");
#define ASSERT_FINITEF(res) _ASSERT_FINITE(isinff, res)
#define ASSERT_FINITE(res) _ASSERT_FINITE(isinf, res)
#define ASSERT_FINITEL(res) _ASSERT_FINITE(isinfl, res)
/**
* Returns the principal value of the arc cosine of x, expressed in radians.
*/
float __acosf_finite(float x)
{
float res = acosf(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the principal value of the arc cosine of x, expressed in radians.
*/
double __acos_finite(double x)
{
double res = acos(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the principal value of the arc cosine of x, expressed in radians.
*/
long double __acosl_finite(long double x)
{
long double res = acosl(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the nonnegative area hyperbolic cosine of x.
*/
double __acosh_finite(double x)
{
double res = acosh(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the nonnegative area hyperbolic cosine of x.
*/
float __acoshf_finite(float x)
{
float res = acoshf(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the nonnegative area hyperbolic cosine of x.
*/
long double __acoshl_finite(long double x)
{
long double res = acoshl(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the principal value of the arc sine of x, expressed in radians.
*/
float __asinf_finite(float x)
{
float res = asinf(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the principal value of the arc sine of x, expressed in radians.
*/
double __asin_finite(double x)
{
double res = asin(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the principal value of the arc sine of x, expressed in radians.
*/
long double __asinl_finite(long double x)
{
long double res = asinl(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the principal value of the arc tangent of x/y, expressed in radians.
*/
float __atan2f_finite(float x, float y)
{
float res = atan2f(x, y);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the principal value of the arc tangent of x/y, expressed in radians.
*/
double __atan2_finite(double x, double y)
{
double res = atan2(x, y);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the principal value of the arc tangent of x/y, expressed in radians.
*/
long double __atan2l_finite(long double x, long double y)
{
long double res = atan2l(x, y);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the area hyperbolic tangent of x.
*/
float __atanhf_finite(float x)
{
float res = atanhf(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the area hyperbolic tangent of x.
*/
double __atanh_finite(double x)
{
double res = atanh(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the area hyperbolic tangent of x.
*/
long double __atanhl_finite(long double x)
{
long double res = atanhl(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the hyperbolic cosine of x.
*/
float __coshf_finite(float x)
{
float res = coshf(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the hyperbolic cosine of x.
*/
double __cosh_finite(double x)
{
double res = cosh(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the hyperbolic cosine of x.
*/
long double __coshl_finite(long double x)
{
long double res = coshl(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Return the value of 10 raised to the power of x.
*/
float __exp10f_finite(float x)
{
float res = exp10f(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Return the value of 10 raised to the power of x.
*/
double __exp10_finite(double x)
{
double res = exp10(x);
ASSERT_FINITE(res);
return res;
}
/**
* Return the value of 10 raised to the power of x.
*/
long double __exp10l_finite(long double x)
{
long double res = exp10l(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the base-2 exponential function of x, which is 2 raised to the power x
*/
float __exp2f_finite(float x)
{
float res = exp2f(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the base-2 exponential function of x, which is 2 raised to the power x
*/
double __exp2_finite(double x)
{
double res = exp2(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the base-2 exponential function of x, which is 2 raised to the power x
*/
long double __exp2l_finite(long double x)
{
long double res = exp2l(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the base-e exponential function of x, which is e raised to the power x
*/
float __expf_finite(float x)
{
float res = expf(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the base-e exponential function of x, which is e raised to the power x
*/
double __exp_finite(double x)
{
double res = exp(x);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the base-e exponential function of x, which is e raised to the power x
*/
long double __expl_finite(long double x)
{
long double res = expl(x);
ASSERT_FINITEL(res);
return res;
}
/**
* Returns the floating-point remainder of x/y (rounded towards zero)
*/
float __fmodf_finite(float x, float y)
{
float res = fmodf(x, y);
ASSERT_FINITEF(res);
return res;
}
/**
* Returns the floating-point remainder of x/y (rounded towards zero)
*/
double __fmod_finite(double x, double y)
{
double res = fmod(x, y);
ASSERT_FINITE(res);
return res;
}
/**
* Returns the floating-point remainder of x/y (rounded towards zero)
*/
long double __fmodl_finite(long double x, long double y)
{
long double res = fmodl(x, y);
ASSERT_FINITEL(res);
return res;
}
/**
* Computes the square root of the sum of the squares of x and y, without undue
* overflow or underflow at intermediate stages of the computation.
*/
float __hypotf_finite(float x, float y)
{
float res = hypotf(x, y);
ASSERT_FINITEF(res);
return res;
}
/**
* Computes the square root of the sum of the squares of x and y, without undue
* overflow or underflow at intermediate stages of the computation.
*/
double __hypot_finite(double x, double y)
{
double res = hypot(x, y);
ASSERT_FINITE(res);
return res;
}
/**
* Computes the square root of the sum of the squares of x and y, without undue
* overflow or underflow at intermediate stages of the computation.
*/
long double __hypotl_finite(long double x, long double y)
{
long double res = hypotl(x, y);
ASSERT_FINITEL(res);
return res;
}
/**
* Return Bessel functions of x of the first kind of orders 0.
*/
float __j0f_finite(float x)
{
float res = j0f(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Return Bessel functions of x of the first kind of orders 0.
*/
double __j0_finite(double x)
{
double res = j0(x);
ASSERT_FINITE(res);
return res;
}
/* The below requires support for j0l, see above */
#if 0
/**
* Return Bessel functions of x of the first kind of orders 0.
*/
long double __j0l_finite(long double x)
{
long double res = j0l(x);
ASSERT_FINITEL(res);
return res;
}
#endif
/**
* Return Bessel functions of x of the first kind of orders 1.
*/
float __j1f_finite(float x)
{
float res = j1f(x);
ASSERT_FINITEF(res);
return res;
}
/**
* Return Bessel functions of x of the first kind of orders 1.
*/
double __j1_finite(double x)
{
double res = j1(x);
ASSERT_FINITE(res);
return res;
}
/* The below requires support for j1l, see above */
#if 0
/**
* Return Bessel functions of x of the first kind of orders 1.
*/
long double __j1l_finite(long double x)
{
long double res = j1l(x);
ASSERT_FINITEL(res);
return res;
}
#endif
/**
* Return the Bessel function of x of the first kind of order n.
*/
float __jnf_finite(int n, float x)
{
float res = jnf(n, x);
ASSERT_FINITEF(res);
return res;
}
/**
* Return the Bessel function of x of the first kind of order n.
*/
double __jn_finite(int n, double x)
{
double res = jn(n, x);
ASSERT_FINITE(res);
return res;
}
/* The below requires support for jnl, see above */
#if 0
/**
* Return the Bessel function of x of the first kind of order n.
*/
long double __jnl_finite(int n, long double x)
{
long double res = jnl(n, x);
ASSERT_FINITEL(res);
return res;
}
#endif
/**
* Returns the natural logarithm of the absolute value of the Gamma function.
*/
float __lgammaf_finite(float x)
{
float res = lgammaf(x);
ASSERT_FINITEF(res);
return res;
}
alias(__lgammaf_finite, __gammaf_finite);
/**
* Returns the natural logarithm of the absolute value of the Gamma function.
*/
double __lgamma_finite(double x)
{
double res = lgamma(x);
ASSERT_FINITE(res);
return res;
}
alias(__lgamma_finite, __gamma_finite);
/**
* Returns the natural logarithm of the absolute value of the Gamma function.
*/
long double __lgammal_finite(long double x)
{
long double res = lgammal(x);
ASSERT_FINITEL(res);