mirror of https://github.com/flysand7/ciabatta.git
410 lines
9.3 KiB
C
410 lines
9.3 KiB
C
|
|
#include <math.h>
|
|
#include <stdint.h>
|
|
#include <fenv.h>
|
|
#include <float.h>
|
|
#include <limits.h>
|
|
|
|
#include <_compiler.h>
|
|
#if defined(_compiler_clang) || defined(_compiler_gnu)
|
|
#define just_do_it(t) __attribute__((unused)) volatile t
|
|
#endif
|
|
|
|
#define asuint64(x) ((union {double f; uint64_t i;}){x}).i
|
|
#define asdouble(x) ((union {double f; uint64_t i;}){x}).f
|
|
|
|
double nearbyint(double x) {
|
|
#pragma STDC FENV_ACCESS ON
|
|
int e = fetestexcept(FE_INEXACT);
|
|
x = rint(x);
|
|
if (!e) feclearexcept(FE_INEXACT);
|
|
return x;
|
|
}
|
|
|
|
float nearbyintf(float x) {
|
|
#pragma STDC FENV_ACCESS ON
|
|
int e = fetestexcept(FE_INEXACT);
|
|
x = rintf(x);
|
|
if (!e) feclearexcept(FE_INEXACT);
|
|
return x;
|
|
}
|
|
|
|
long double nearbyintl(long double x) {
|
|
return nearbyint(x);
|
|
}
|
|
|
|
double nextafter(double x, double y) {
|
|
union {double f; uint64_t i;} ux={x}, uy={y};
|
|
uint64_t ax, ay;
|
|
int e;
|
|
if (isnan(x) || isnan(y)) return x + y;
|
|
if (ux.i == uy.i) return y;
|
|
ax = ux.i & -1ULL/2;
|
|
ay = uy.i & -1ULL/2;
|
|
if (ax == 0) {
|
|
if (ay == 0) return y;
|
|
ux.i = (uy.i & 1ULL<<63) | 1;
|
|
} else if (ax > ay || ((ux.i ^ uy.i) & 1ULL<<63)) {
|
|
ux.i--;
|
|
}
|
|
else {
|
|
ux.i++;
|
|
}
|
|
e = ux.i >> 52 & 0x7ff;
|
|
/* raise overflow if ux.f is infinite and x is finite */
|
|
if (e == 0x7ff) just_do_it(float) _x = x+x;
|
|
/* raise underflow if ux.f is subnormal or zero */
|
|
if (e == 0) just_do_it(float) _x = x*x + ux.f*ux.f;
|
|
return ux.f;
|
|
}
|
|
|
|
float nextafterf(float x, float y) {
|
|
union {float f; uint32_t i;} ux={x}, uy={y};
|
|
uint32_t ax, ay, e;
|
|
|
|
if (isnan(x) || isnan(y)) return x + y;
|
|
if (ux.i == uy.i) return y;
|
|
ax = ux.i & 0x7fffffff;
|
|
ay = uy.i & 0x7fffffff;
|
|
if (ax == 0) {
|
|
if (ay == 0) return y;
|
|
ux.i = (uy.i & 0x80000000) | 1;
|
|
} else if (ax > ay || ((ux.i ^ uy.i) & 0x80000000)) {
|
|
ux.i--;
|
|
}
|
|
else {
|
|
ux.i++;
|
|
}
|
|
e = ux.i & 0x7f800000;
|
|
/* raise overflow if ux.f is infinite and x is finite */
|
|
if (e == 0x7f800000) just_do_it(float) _x = x+x;
|
|
/* raise underflow if ux.f is subnormal or zero */
|
|
if (e == 0) just_do_it(float) _x = x*x + ux.f*ux.f;
|
|
return ux.f;
|
|
}
|
|
|
|
long double nextafterl(long double x, long double y) {
|
|
return nextafter(x, y);
|
|
}
|
|
|
|
double nexttoward(double x, long double y) {
|
|
return nextafter(x, y);
|
|
}
|
|
|
|
float nexttowardf(float x, long double y) {
|
|
union {float f; uint32_t i;} ux = {x};
|
|
uint32_t e;
|
|
if (isnan(x) || isnan(y)) return x + y;
|
|
if (x == y) return y;
|
|
if (x == 0) {
|
|
ux.i = 1;
|
|
if (signbit(y)) ux.i |= 0x80000000;
|
|
} else if (x < y) {
|
|
if (signbit(x)) ux.i--;
|
|
else ux.i++;
|
|
} else {
|
|
if (signbit(x)) ux.i++;
|
|
else ux.i--;
|
|
}
|
|
e = ux.i & 0x7f800000;
|
|
/* raise overflow if ux.f is infinite and x is finite */
|
|
if (e == 0x7f800000) just_do_it(float) _x = x+x;
|
|
/* raise underflow if ux.f is subnormal or zero */
|
|
if (e == 0) just_do_it(float) _x = x*x + ux.f*ux.f;
|
|
return ux.f;
|
|
}
|
|
|
|
long double nexttowardl(long double x, long double y) {
|
|
return nextafterl(x, y);
|
|
}
|
|
|
|
double rint(double x) {
|
|
static const double_t toint = 1/DBL_EPSILON;
|
|
union {double f; uint64_t i;} u = {x};
|
|
int e = u.i>>52 & 0x7ff;
|
|
int s = u.i>>63;
|
|
double y;
|
|
if (e >= 0x3ff+52) return x;
|
|
if (s) y = x - toint + toint;
|
|
else y = x + toint - toint;
|
|
if (y == 0) return s ? -0.0 : +0.0;
|
|
return y;
|
|
}
|
|
|
|
float rintf(float x) {
|
|
static const float toint = 1/FLT_EPSILON;
|
|
union {float f; uint32_t i;} u = {x};
|
|
int e = u.i>>23 & 0xff;
|
|
int s = u.i>>31;
|
|
float y;
|
|
if (e >= 0x7f+23) return x;
|
|
if (s) y = x - toint + toint;
|
|
else y = x + toint - toint;
|
|
if (y == 0) return s ? -0.0f : 0.0f;
|
|
return y;
|
|
}
|
|
|
|
long double rintl(long double x) {
|
|
return rint(x);
|
|
}
|
|
|
|
#if LONG_MAX < 1U<<53 && defined(FE_INEXACT)
|
|
static long lrint_slow(double x)
|
|
{
|
|
#pragma STDC FENV_ACCESS ON
|
|
int e;
|
|
e = fetestexcept(FE_INEXACT);
|
|
x = rint(x);
|
|
if (!e && (x > LONG_MAX || x < LONG_MIN))
|
|
feclearexcept(FE_INEXACT);
|
|
/* conversion */
|
|
return x;
|
|
}
|
|
|
|
long lrint(double x)
|
|
{
|
|
uint32_t abstop = asuint64(x)>>32 & 0x7fffffff;
|
|
uint64_t sign = asuint64(x) & (1ULL << 63);
|
|
|
|
if (abstop < 0x41dfffff) {
|
|
/* |x| < 0x7ffffc00, no overflow */
|
|
double_t toint = asdouble(asuint64(1/DBL_EPSILON) | sign);
|
|
double_t y = x + toint - toint;
|
|
return (long)y;
|
|
}
|
|
return lrint_slow(x);
|
|
}
|
|
#else
|
|
long lrint(double x) {
|
|
return rint(x);
|
|
}
|
|
#endif
|
|
|
|
long lrintf(float x) {
|
|
return rintf(x);
|
|
}
|
|
|
|
long lrintl(long double x) {
|
|
return lrint(x);
|
|
}
|
|
|
|
long long llrint(double x) {
|
|
return rint(x);
|
|
}
|
|
|
|
long long llrintf(float x) {
|
|
return rintf(x);
|
|
}
|
|
|
|
long long llrintl(long double x) {
|
|
return llrint(x);
|
|
}
|
|
|
|
double round(double x) {
|
|
static const double_t toint = 1/DBL_EPSILON;
|
|
union {double f; uint64_t i;} u = {x};
|
|
int e = u.i >> 52 & 0x7ff;
|
|
double_t y;
|
|
if (e >= 0x3ff+52) return x;
|
|
if (u.i >> 63) x = -x;
|
|
if (e < 0x3ff-1) {
|
|
/* raise inexact if x!=0 */
|
|
just_do_it(float) _x = x + toint;
|
|
return 0*u.f;
|
|
}
|
|
y = x + toint - toint - x;
|
|
if (y > 0.5) y = y + x - 1;
|
|
else if (y <= -0.5) y = y + x + 1;
|
|
else y = y + x;
|
|
if (u.i >> 63) y = -y;
|
|
return y;
|
|
}
|
|
|
|
float roundf(float x) {
|
|
static const double_t toint = 1/FLT_EPSILON;
|
|
union {float f; uint32_t i;} u = {x};
|
|
int e = u.i >> 23 & 0xff;
|
|
float_t y;
|
|
if (e >= 0x7f+23) return x;
|
|
if (u.i >> 31) x = -x;
|
|
if (e < 0x7f-1) {
|
|
just_do_it(float) _x = x + toint;
|
|
return 0*u.f;
|
|
}
|
|
y = x + toint - toint - x;
|
|
if (y > 0.5f) y = y + x - 1;
|
|
else if (y <= -0.5f) y = y + x + 1;
|
|
else y = y + x;
|
|
if (u.i >> 31) y = -y;
|
|
return y;
|
|
}
|
|
|
|
long double roundl(long double x) {
|
|
return round(x);
|
|
}
|
|
|
|
long lround(double x) {
|
|
return round(x);
|
|
}
|
|
|
|
long lroundf(float x) {
|
|
return roundf(x);
|
|
}
|
|
|
|
long lroundl(long double x) {
|
|
return roundl(x);
|
|
}
|
|
|
|
long long llround(double x) {
|
|
return round(x);
|
|
}
|
|
|
|
long long llroundf(float x) {
|
|
return roundf(x);
|
|
}
|
|
|
|
long long llroundl(long double x) {
|
|
return roundl(x);
|
|
}
|
|
|
|
double ceil(double x) {
|
|
static const double_t toint = 1/DBL_EPSILON;
|
|
union {double f; uint64_t i;} u = {x};
|
|
int e = u.i >> 52 & 0x7ff;
|
|
double_t y;
|
|
|
|
if (e >= 0x3ff+52 || x == 0)
|
|
return x;
|
|
/* y = int(x) - x, where int(x) is an integer neighbor of x */
|
|
if (u.i >> 63)
|
|
y = x - toint + toint - x;
|
|
else
|
|
y = x + toint - toint - x;
|
|
/* special case because of non-nearest rounding modes */
|
|
if (e <= 0x3ff-1) {
|
|
just_do_it(double) _x = y;
|
|
return u.i >> 63 ? -0.0 : 1;
|
|
}
|
|
if (y < 0)
|
|
return x + y + 1;
|
|
return x + y;
|
|
}
|
|
|
|
float ceilf(float x) {
|
|
union {float f; uint32_t i;} u = {x};
|
|
int e = (int)(u.i >> 23 & 0xff) - 0x7f;
|
|
uint32_t m;
|
|
|
|
if (e >= 23)
|
|
return x;
|
|
if (e >= 0) {
|
|
m = 0x007fffff >> e;
|
|
if ((u.i & m) == 0)
|
|
return x;
|
|
just_do_it(float) _x = (x + 0x1p120f);
|
|
if (u.i >> 31 == 0)
|
|
u.i += m;
|
|
u.i &= ~m;
|
|
} else {
|
|
just_do_it(float) _x = (x + 0x1p120f);
|
|
if (u.i >> 31)
|
|
u.f = -0.0;
|
|
else if (u.i << 1)
|
|
u.f = 1.0;
|
|
}
|
|
return u.f;
|
|
}
|
|
|
|
long double ceill(long double x) {
|
|
return ceil(x);
|
|
}
|
|
|
|
double floor(double x) {
|
|
static const double_t toint = 1/DBL_EPSILON;
|
|
union {double f; uint64_t i;} u = {x};
|
|
int e = u.i >> 52 & 0x7ff;
|
|
double_t y;
|
|
if (e >= 0x3ff+52 || x == 0)
|
|
return x;
|
|
/* y = int(x) - x, where int(x) is an integer neighbor of x */
|
|
if (u.i >> 63)
|
|
y = x - toint + toint - x;
|
|
else
|
|
y = x + toint - toint - x;
|
|
/* special case because of non-nearest rounding modes */
|
|
if (e <= 0x3ff-1) {
|
|
just_do_it(double) _x = (y);
|
|
return u.i >> 63 ? -1 : 0;
|
|
}
|
|
if (y > 0)
|
|
return x + y - 1;
|
|
return x + y;
|
|
}
|
|
|
|
float floorf(float x) {
|
|
union {float f; uint32_t i;} u = {x};
|
|
int e = (int)(u.i >> 23 & 0xff) - 0x7f;
|
|
uint32_t m;
|
|
|
|
if (e >= 23)
|
|
return x;
|
|
if (e >= 0) {
|
|
m = 0x007fffff >> e;
|
|
if ((u.i & m) == 0)
|
|
return x;
|
|
just_do_it(float) _x = (x + 0x1p120f);
|
|
if (u.i >> 31)
|
|
u.i += m;
|
|
u.i &= ~m;
|
|
} else {
|
|
just_do_it(float) _x = (x + 0x1p120f);
|
|
if (u.i >> 31 == 0)
|
|
u.i = 0;
|
|
else if (u.i << 1)
|
|
u.f = -1.0;
|
|
}
|
|
return u.f;
|
|
}
|
|
|
|
long double floorl(long double x) {
|
|
return floor(x);
|
|
}
|
|
|
|
double trunc(double x) {
|
|
union {double f; uint64_t i;} u = {x};
|
|
int e = (int)(u.i >> 52 & 0x7ff) - 0x3ff + 12;
|
|
uint64_t m;
|
|
|
|
if (e >= 52 + 12)
|
|
return x;
|
|
if (e < 12)
|
|
e = 1;
|
|
m = -1ULL >> e;
|
|
if ((u.i & m) == 0)
|
|
return x;
|
|
just_do_it(double) _x = (x + 0x1p120f);
|
|
u.i &= ~m;
|
|
return u.f;
|
|
}
|
|
|
|
float truncf(float x) {
|
|
union {float f; uint32_t i;} u = {x};
|
|
int e = (int)(u.i >> 23 & 0xff) - 0x7f + 9;
|
|
uint32_t m;
|
|
|
|
if (e >= 23 + 9)
|
|
return x;
|
|
if (e < 9)
|
|
e = 1;
|
|
m = -1U >> e;
|
|
if ((u.i & m) == 0)
|
|
return x;
|
|
just_do_it(float) _x = (x + 0x1p120f);
|
|
u.i &= ~m;
|
|
return u.f;
|
|
}
|
|
|
|
long double truncl(long double x) {
|
|
return trunc(x);
|
|
}
|