#include "math.h"
#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
#include "DD.h"
#include "fp_private.h"
#include "limits.h"
extern long double sqrtl(long double);
#define REM_NAN "9"
#define kSignMask 0x80000000u
#define kExpMask 0x7ff00000u
static const Double kTZ = {{0x0,0x1}};
static const Double kUP = {{0x0,0x2}};
static const Double kDN = {{0x0,0x3}};
static const float kTwoTo52 = 4503599627370496.0f;
long double logbl( long double x )
{
DblDbl xDD;
double fpenv, sum;
xDD.f = x;
FEGETENVD(fpenv);
FESETENVD(kTZ.f); sum = xDD.d[0] + xDD.d[1]; FESETENVD(fpenv);
return logb(sum);
}
int ilogbl( long double x )
{
DblDbl xDD;
double fpenv, sum;
xDD.f = x;
FEGETENVD(fpenv);
FESETENVD(kTZ.f); sum = xDD.d[0] + xDD.d[1]; FESETENVD(fpenv);
return ilogb(sum);
}
long double scalbnl( long double x, int iscale )
{
DblDbl xDD;
double fpenv;
double Z,z;
uint32_t hibits;
FEGETENVD(fpenv);
FESETENVD(0.0);
xDD.f = x;
xDD.d[0] = scalbn(xDD.d[0],iscale); hibits = (xDD.i[0] >> 20) & 0x7ffu; if ((hibits != 0u) && (hibits != 0x7ffu)) { z = scalbn(xDD.d[1],iscale); FESETENVD(kTZ.f); Z = xDD.d[0] + z; xDD.d[1] = (xDD.d[0] - Z) + z;
xDD.d[0] = Z;
FESETENVD(fpenv);
return (xDD.f); }
Z = xDD.d[0];
if ((Z != Z) || (Z == 0.0)) xDD.d[1] = Z;
else if (hibits == 0x7ffu) xDD.d[1] = 0.0;
else xDD.d[1] = scalbn(xDD.d[1],iscale);
FESETENVD(fpenv);
return (xDD.f); }
long double scalblnl( long double x, long int n )
{
int m;
if (unlikely(n > 2097))
m = 2098;
else if (unlikely(n < -2098))
m = -2099;
else
m = (int) n;
return scalbnl(x, m);
}
long int lrintl ( long double x )
{
long double t;
long int result;
double fenv;
if (unlikely(x != x))
{
feraiseexcept(FE_INVALID);
return LONG_MIN;
}
FEGETENVD(fenv);
t = rintl ( x );
FESETENVD(fenv);
if ( t < (long double)LONG_MIN )
{
feraiseexcept(FE_INVALID);
result = LONG_MIN;
}
else if ( t > (long double)LONG_MAX )
{
feraiseexcept(FE_INVALID);
result = LONG_MAX;
}
else if (t != x)
{
feraiseexcept(FE_INEXACT);
result = (long int) t;
}
else
{
result = (long int) t;
}
return result;
}
long double hypotl( long double x, long double y )
{
DblDbl xDD, yDD;
double fpenv;
long double ldtemp;
uint32_t expx, expy;
const Float kINF = {0x7f800000};
FEGETENVD(fpenv);
FESETENVD(0.0);
xDD.f = x;
yDD.f = y;
if (xDD.i[0] >= kSignMask) {
xDD.i[0] ^= kSignMask;
xDD.i[2] ^= kSignMask;
}
if (yDD.i[0] >= kSignMask) {
yDD.i[0] ^= kSignMask;
yDD.i[2] ^= kSignMask;
}
expx = (xDD.i[0] >> 20) & 0x7ffu;
expy = (yDD.i[0] >> 20) & 0x7ffu;
if ((expx == 0x7ffu) || (expy == 0x7ffu)) {
if (xDD.d[0] == kINF.f) { FESETENVD(fpenv);
return (xDD.f);
}
else if (yDD.d[0] == kINF.f) {
FESETENVD(fpenv);
return (yDD.f);
}
else {
ldtemp = x + y;
FESETENVD(fpenv);
return ldtemp; }
}
if (yDD.f > xDD.f) { ldtemp = yDD.f;
yDD.f = xDD.f;
xDD.f = ldtemp;
}
if (yDD.d[0] == 0.0) { FESETENVD(fpenv);
return (xDD.f);
}
ldtemp = yDD.f / xDD.f; ldtemp *= ldtemp;
ldtemp = xDD.f*sqrtl(1.0L + ldtemp);
FESETENVD(fpenv);
return (ldtemp);
}
long double truncl( long double x )
{
DblDbl ldu;
double xd, result, fpenv;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = x;
if ((ldu.i[0] & kExpMask) == kExpMask) { FESETENVD(fpenv);
return (x);
}
if (fabs(ldu.d[1]) >= kTwoTo52) { FESETENVD(fpenv);
return (x);
}
FESETENVD(kTZ.f); xd = ldu.d[0] + ldu.d[1];
if (fabs(xd) < kTwoTo52) { ldu.d[1] = 0.0; if (xd >= 0.0) result = (xd + kTwoTo52) - kTwoTo52;
else
result = (xd - kTwoTo52) + kTwoTo52;
if (result == 0.0) { if (ldu.i[0] < kSignMask)
result = +0.0;
else {
result = -0.0;
ldu.d[1] = -0.0;
}
}
ldu.d[0] = result;
FESETENVD(fpenv); return (ldu.f); }
if (xd > 0.0) {
FESETENVD(kDN.f);
}
else
{
FESETENVD(kUP.f);
}
if (ldu.d[1] >= 0.0)
result = (ldu.d[1] + kTwoTo52) - kTwoTo52;
else
result = (ldu.d[1] - kTwoTo52) + kTwoTo52;
FESETENVD(0.0);
xd = ldu.d[0] + result;
ldu.d[1] = (ldu.d[0] - xd) + result;
ldu.d[0] = xd;
FESETENVD(fpenv);
return (ldu.f); }
long double floorl( long double x )
{
DblDbl ldu;
double xd, result, fpenv;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = x;
if ((ldu.i[0] & kExpMask) == kExpMask) { FESETENVD(fpenv);
return (x);
}
if (fabs(ldu.d[1]) >= kTwoTo52) { FESETENVD(fpenv);
return (x);
}
FESETENVD(kDN.f);
if (fabs(ldu.d[0]) < kTwoTo52) { xd = ldu.d[0] + ldu.d[1]; ldu.d[1] = 0.0; if (ldu.d[0] >= 0.0) result = (xd + kTwoTo52) - kTwoTo52;
else
result = (xd - kTwoTo52) + kTwoTo52;
if (result == 0.0) { if (ldu.i[0] < kSignMask)
result = 0.0;
else {
result = -0.0;
ldu.d[1] = result;
}
}
ldu.d[0] = result;
FESETENVD(fpenv); return (ldu.f); }
if (ldu.d[1] >= 0.0) result = (ldu.d[1] + kTwoTo52) - kTwoTo52;
else
result = (ldu.d[1] - kTwoTo52) + kTwoTo52;
FESETENVD(0.0);
xd = ldu.d[0] + result;
ldu.d[1] = (ldu.d[0] - xd) + result;
ldu.d[0] = xd;
FESETENVD(fpenv);
return (ldu.f); }
long double ceill( long double x )
{
DblDbl ldu;
double xd, result, fpenv;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = x;
if ((ldu.i[0] & kExpMask) == kExpMask) { FESETENVD(fpenv);
return (x);
}
if (fabs(ldu.d[1]) >= kTwoTo52) { FESETENVD(fpenv);
return (x);
}
FESETENVD(kUP.f);
if (fabs(ldu.d[0]) < kTwoTo52) { xd = ldu.d[0] + ldu.d[1]; ldu.d[1] = 0.0; if (ldu.d[0] >= 0.0) result = (xd + kTwoTo52) - kTwoTo52;
else
result = (xd - kTwoTo52) + kTwoTo52;
if (result == 0.0) { if (ldu.i[0] < kSignMask)
result = 0.0;
else {
result = -0.0;
ldu.d[1] = result;
}
}
ldu.d[0] = result;
FESETENVD(fpenv);
return (ldu.f); }
if (ldu.d[1] >= 0.0) result = (ldu.d[1] + kTwoTo52) - kTwoTo52;
else
result = (ldu.d[1] - kTwoTo52) + kTwoTo52;
FESETENVD(0.0);
xd = ldu.d[0] + result;
ldu.d[1] = (ldu.d[0] - xd) + result;
ldu.d[0] = xd;
FESETENVD(fpenv); return (ldu.f); }
long double rintl( long double x )
{
DblDbl ldu;
double fpenv;
Double fenv;
double result,xh,xt;
uint32_t rnddir;
FEGETENVD(fpenv);
FESETENVD(0.0);
fenv.f = fpenv;
rnddir = fenv.i[1] & FE_ALL_RND;
if (rnddir == FE_TONEAREST) { ldu.f = x;
if ((ldu.i[0] & kExpMask) == kExpMask) {
FESETENVD(fpenv);
return (x);
}
if (fabs(ldu.d[1]) >= kTwoTo52) {
FESETENVD(fpenv);
return (x);
}
xh = ldu.d[0]; xt = ldu.d[1];
if (fabs(xh) < kTwoTo52) { ldu.d[1] = 0.0; if (xh >= 0.0) result = (xh + kTwoTo52) - kTwoTo52;
else
result = (xh - kTwoTo52) + kTwoTo52;
if ((fabs(result - xh) == 0.5) && (xt != 0.0))
result = (xt > 0.0) ? (xh + 0.5) : (xh - 0.5);
if (result == 0.0) { if (ldu.i[0] < kSignMask) {
result = 0.0;
}
else {
result = -0.0;
ldu.d[1] = result;
}
}
ldu.d[0] = result;
FESETENVD(fpenv);
return (ldu.f); }
if (xt >= 0.0) result = (xt + kTwoTo52) - kTwoTo52;
else
result = (xt - kTwoTo52) + kTwoTo52;
ldu.d[0] = xh + result; ldu.d[1] = xh - ldu.d[0] + result;
FESETENVD(fpenv);
return (ldu.f); }
FESETENVD(fpenv);
if (rnddir == FE_TOWARDZERO) return (truncl(x));
else if (rnddir == FE_UPWARD) return (ceill(x));
else return (floorl(x));
}
long long int llrintl ( long double x )
{
long double t;
long long int result;
double fenv;
if (unlikely(x != x))
{
feraiseexcept(FE_INVALID);
return LONG_LONG_MIN;
}
FEGETENVD(fenv);
t = rintl ( x );
FESETENVD(fenv);
if ( t < (long double)LONG_LONG_MIN )
{
feraiseexcept(FE_INVALID);
result = LONG_LONG_MIN;
}
else if ( t > (long double)LONG_LONG_MAX )
{
feraiseexcept(FE_INVALID);
result = LONG_LONG_MAX;
}
else if (t != x)
{
feraiseexcept(FE_INEXACT);
result = (long long int) t;
}
else
{
result = (long long int) t;
}
return result;
}
long double nearbyintl( long double x )
{
return (rintl(x));
}
long double roundl( long double x )
{
DblDbl ldu;
double xh, xt, xd, result, fpenv;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = x;
xh = ldu.d[0];
xt = ldu.d[1];
if ((ldu.i[0] & kExpMask) == kExpMask) { FESETENVD(fpenv);
return (x);
}
if (fabs(xt) >= kTwoTo52) { FESETENVD(fpenv);
return (x);
}
FESETENVD(kTZ.f);
if (fabs(xh) < kTwoTo52) { ldu.d[1] = 0.0; if (xh >= 0.0) {
xd = (xt + 0.5) + xh;
result = (xd + kTwoTo52) - kTwoTo52;
}
else {
xd = (xt - 0.5) + xh;
result = (xd - kTwoTo52) + kTwoTo52;
}
if (result == 0.0) { if (ldu.i[0] < kSignMask)
result = 0.0;
else {
result = -0.0;
ldu.d[1] = result;
}
}
ldu.d[0] = result;
FESETENVD(fpenv); return (ldu.f); }
if (xh > 0.0) {
xd = xt + 0.5;
FESETENVD(kDN.f);
}
else {
xd = xt - 0.5;
FESETENVD(kUP.f);
}
if (xd >= 0.0)
result = (xd + kTwoTo52) - kTwoTo52;
else
result = (xd - kTwoTo52) + kTwoTo52;
FESETENVD(0.0);
xd = ldu.d[0] + result;
ldu.d[1] = (ldu.d[0] - xd) + result;
ldu.d[0] = xd;
FESETENVD(fpenv);
return (ldu.f); }
long long int llroundl ( long double x )
{
long double t;
long long int result;
double fenv;
if (unlikely(x != x))
{
feraiseexcept(FE_INVALID);
return LONG_LONG_MAX;
}
FEGETENVD(fenv);
t = roundl ( x );
FESETENVD(fenv);
if ( t < (long double)LONG_LONG_MIN )
{
feraiseexcept(FE_INVALID);
result = LONG_LONG_MIN;
}
else if ( t > (long double)LONG_LONG_MAX )
{
feraiseexcept(FE_INVALID);
result = LONG_LONG_MAX;
}
else if (t != x)
{
feraiseexcept(FE_INEXACT);
result = (long long int) t;
}
else
{
result = (long long int) t;
}
return result;
}
long double fdiml(long double x, long double y)
{
double fenv;
long double result;
FEGETENVD(fenv);
FESETENVD(0.0);
if (unlikely((x != x) || (y != y)))
{
FESETENVD(fenv);
return ( x + y );
}
else if (x > y)
result = (x - y);
else
result = 0.0L;
FESETENVD(fenv);
return result;
}
long double fmaxl( long double x, long double y )
{
DblDbl xDD,yDD;
double fpenv;
FEGETENVD(fpenv);
FESETENVD(0.0);
xDD.f = x;
yDD.f = y;
if (yDD.d[0] != yDD.d[0]) {
FESETENVD(fpenv);
return (x);
}
if (xDD.d[0] != xDD.d[0]) {
FESETENVD(fpenv);
return (y);
}
if (y > x) {
FESETENVD(fpenv);
return (y);
}
else {
FESETENVD(fpenv);
return (x);
}
}
long double fminl( long double x, long double y )
{
DblDbl xDD,yDD;
double fpenv;
FEGETENVD(fpenv);
FESETENVD(0.0);
xDD.f = x;
yDD.f = y;
if (yDD.d[0] != yDD.d[0]) {
FESETENVD(fpenv);
return (x);
}
if (xDD.d[0] != xDD.d[0]) {
FESETENVD(fpenv);
return (y);
}
if (y < x) {
FESETENVD(fpenv);
return (y);
}
else {
FESETENVD(fpenv);
return (x);
}
}
long double fmal( long double x, long double y, long double z )
{
DblDbl xDD, yDD, zDD;
double th, tm, tl, u;
xDD.f = x;
yDD.f = y;
zDD.f = z;
_d3mul( xDD.d[0], xDD.d[1], 0.0, yDD.d[0], yDD.d[1], 0.0, &th, &tm, &tl );
_d3add( th, tm, tl, zDD.d[0], zDD.d[1], 0.0, &xDD.d[0], &xDD.d[1], &u );
return xDD.f;
}
long double fabsl( long double x )
{
DblDbl xDD;
xDD.f = x;
if (xDD.i[0] < kSignMask) return (x);
else { xDD.i[0] ^= kSignMask;
xDD.i[2] ^= kSignMask;
return (xDD.f);
}
}
long double copysignl( long double x, long double y )
{
DblDbl xDD,yDD;
xDD.f = x;
yDD.f = y;
if ((xDD.i[0] & kSignMask) == (yDD.i[0] & kSignMask))
return (x);
else { xDD.i[0] ^= kSignMask;
xDD.i[2] ^= kSignMask;
return (xDD.f);
}
}
long rinttoll( long double x )
{
DblDbl ldu;
uint32_t expx, lowbit;
long temp;
double fpenv;
ldu.f = x;
expx = (ldu.i[0] >> 20) & 0x7ffu; lowbit = ldu.i[1] & 1u;
if ((expx < 0x41fu) && (expx > 0x3fdu)) { if ((lowbit == 0) && (ldu.d[1] != 0.0)) {
uint32_t signh = (ldu.i[0] >> 31) & 1u; uint32_t signt = (ldu.i[2] >> 31) & 1u; if (signh == signt) ldu.i[1] |= 1u;
else {
ldu.i[1] -= 1u; if (ldu.i[1] == 0xffffffffu)
ldu.i[0] -= 1u;
}
}
}
FEGETENVD(fpenv);
FESETENVD(0.0);
temp = rinttol(ldu.d[0]);
FESETENVD(fpenv);
return (temp);
}
long int lroundl( long double x )
{
long double t;
long int result;
double fenv;
if (unlikely(x != x))
{
feraiseexcept(FE_INVALID);
return LONG_MAX;
}
FEGETENVD(fenv);
t = roundl ( x );
FESETENVD(fenv);
if ( t < (long double)LONG_MIN )
{
feraiseexcept(FE_INVALID);
result = LONG_MIN;
}
else if ( t > (long double)LONG_MAX )
{
feraiseexcept(FE_INVALID);
result = LONG_MAX;
}
else if (t != x)
{
feraiseexcept(FE_INEXACT);
result = (long int) t;
}
else
{
result = (long int) t;
}
return result;
}
long double ldexpl(long double x, int n)
{
int m;
if (unlikely(n > 2097))
m = 2098;
else if (unlikely(n < -2098))
m = -2099;
else
m = n;
return scalbnl(x, m);
}
long double frexpl(long double x, int *exponent)
{
DblDbl ldu;
double xHead;
double fenv;
long double result;
ldu.f = x;
xHead = ldu.d[0];
FEGETENVD(fenv);
FESETENVD(0.0);
switch (__fpclassifyd(xHead))
{
case FP_NAN:
case FP_INFINITE:
FESETENVD(fenv);
return x;
default:
*exponent = (x == 0.0) ? 0 : (int)(1.0 + logbl(x));
result = scalbnl(x, - (*exponent));
FESETENVD(fenv);
return result;
}
}
long double fmodl(long double xx, long double yy)
{
int iclx,icly;
int32_t iscx,iscy,idiff;
int i;
long double absy,x1,y1,z;
long double rslt;
fenv_t OldEnv;
hexdouble OldEnvironment;
int newexc;
DblDbl xDD, yDD;
xDD.f = xx;
yDD.f = yy;
FEGETENVD( OldEnvironment.d );
FESETENVD( 0.0 );
OldEnv = OldEnvironment.i.lo;
iclx = __fpclassifyd(xDD.d[0]);
icly = __fpclassifyd(yDD.d[0]);
if (likely((iclx & icly) >= FP_NORMAL)) {
x1 = fabsl(xx);
absy = fabsl(yy);
if (absy > x1) {
rslt = xx;
goto ret;
}
else {
iscx = (int32_t) logbl(x1);
iscy = (int32_t) logbl(absy);
idiff = iscx - iscy;
if (idiff != 0) {
y1 = scalbnl(absy,-iscy);
x1 = scalbnl(x1,-iscx);
for (i = idiff; i != 0; i--) {
if ((z = x1 - y1) >= 0) {
x1 = z;
}
x1 += x1;
}
x1 = scalbnl(x1,iscy);
}
if (x1 >= absy) {
x1 -= absy;
}
}
if (xx < 0.0)
x1 = -x1;
rslt = x1;
goto ret;
}
else if ((iclx <= FP_QNAN) || (icly <= FP_QNAN)) {
rslt = xx+yy;
goto ret;
}
else if ((iclx == FP_INFINITE)||(icly == FP_ZERO)) {
rslt = nanl(REM_NAN);
OldEnvironment.i.lo |= SET_INVALID;
FESETENVD ( OldEnvironment.d );
goto ret;
}
else
rslt = xx;
ret:
FEGETENVD (OldEnvironment.d );
newexc = OldEnvironment.i.lo & FE_ALL_EXCEPT;
OldEnvironment.i.lo = OldEnv;
if ((newexc & FE_INVALID) != 0)
OldEnvironment.i.lo |= SET_INVALID;
OldEnvironment.i.lo |= newexc & ( FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW );
FESETENVD (OldEnvironment.d );
return rslt;
}
#warning remquol: cannot gaurantee exact result!
static const hexdouble Huge = HEXDOUBLE(0x7ff00000, 0x00000000);
static const hexdouble HugeHalved = HEXDOUBLE(0x7fe00000, 0x00000000);
long double remquol ( long double x, long double y, int *quo)
{
int iclx,icly;
int32_t iquo;
int32_t iscx, iscy, idiff;
int i;
long double absy,x1,y1,z;
long double rslt;
fenv_t OldEnv;
hexdouble OldEnvironment;
int newexc;
FEGETENVD ( OldEnvironment.d );
FESETENVD ( 0.0 );
__NOOP;
__NOOP;
OldEnv = OldEnvironment.i.lo;
*quo = 0;
iclx = fpclassify(x);
icly = fpclassify(y);
if (likely((iclx & icly) >= FP_NORMAL)) {
x1 = fabsl(x);
absy = fabsl(y);
iquo = 0;
iscx = (int32_t) logbl(x1);
iscy = (int32_t) logbl(absy);
idiff = iscx - iscy;
if (idiff >= 0) {
if (idiff != 0) {
y1 = scalbnl(absy,-iscy);
x1 = scalbnl(x1,-iscx);
for (i = idiff; i != 0; i--) {
if ((z = x1 - y1) >= 0) {
x1 = z;
iquo += 1;
}
iquo += iquo;
x1 += x1;
}
x1 = scalbnl(x1,iscy);
}
if (x1 >= absy) {
x1 -= absy;
iquo +=1;
}
}
if (likely( x1 < HugeHalved.d ))
z = scalbnl( x1, 1 ); else
z = Huge.d;
if ((z > absy) || ((z == absy) && ((iquo & 1) != 0))) {
x1 -= absy;
iquo += 1;
}
if (x < 0.0)
x1 = -x1;
iquo &= 0x0000007f;
if ((signbit(x) ^ signbit(y)) != 0)
iquo = -iquo;
*quo = iquo;
rslt = x1;
goto ret;
}
else if ((iclx <= FP_QNAN) || (icly <= FP_QNAN)) {
rslt = x+y;
goto ret;
}
else if ((iclx == FP_INFINITE)||(icly == FP_ZERO)) {
rslt = nan("9");
OldEnvironment.i.lo |= SET_INVALID;
FESETENVD_GRP( OldEnvironment.d );
goto ret;
}
else
rslt = x;
ret:
FEGETENVD_GRP( OldEnvironment.d );
newexc = OldEnvironment.i.lo & FE_ALL_EXCEPT;
OldEnvironment.i.lo = OldEnv;
if ((newexc & FE_INVALID) != 0)
OldEnvironment.i.lo |= SET_INVALID;
OldEnvironment.i.lo |= newexc & ( FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW );
FESETENVD_GRP( OldEnvironment.d );
return rslt;
}
long double remainderl(long double x, long double y)
{
int quo;
return (remquol(x, y, &quo));
}
long double modfl ( long double x, long double *iptr )
{
DblDbl ldu;
double xh, xt, frach, fract, inth, intt, fpenv;
long double l;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = x;
xh = ldu.d[0];
xt = ldu.d[1];
frach = modf( xh, &inth );
fract = modf( xt, &intt );
*iptr = (long double)(inth + intt);
l = ((long double)frach) + ((long double)fract);
if (l >= 1.0L)
{
l = l - 1.0L;
*iptr = *iptr + 1.0L;
}
FESETENVD(fpenv);
return l;
}
long double nextafterl(long double x, long double y)
{
static const hexdouble EPSILON = HEXDOUBLE(0x00000000, 0x00000001);
static const hexdouble PosBigHead = HEXDOUBLE(0x7fefffff, 0xffffffff);
static const hexdouble PosBigTail = HEXDOUBLE(0x7c9fffff, 0xffffffff);
DblDbl ldu;
double fpenv;
if (unlikely( ( x != x ) || ( y != y ) )) return x + y;
else if ( x == y )
{
if ( ( x == 0.0L ) && ( y == 0.0L )) return y;
else
return x;
}
else if (unlikely( isinf( x ) )) {
ldu.d[0] = PosBigHead.d;
ldu.d[1] = PosBigTail.d;
return copysignl( ldu.f, x );
}
else if (x == 0.0L)
{
ldu.d[0] = EPSILON.d;
ldu.d[1] = 0.0L;
return copysignl( ldu.f, y );
}
else if ( ( ( x < 0.0 ) && ( x < y ) ) || ( ( x > 0.0 ) && ( x > y ) ) )
{
FEGETENVD(fpenv);
FESETENVD(kTZ.f);
x = (x < 0.0L) ? x + (long double)EPSILON.d : x - (long double)EPSILON.d;
FESETENVD(fpenv);
return x;
}
else if ( ( x < 0.0 ) && ( x > y ) )
{
FEGETENVD(fpenv);
FESETENVD(kDN.f);
x = x - (long double)EPSILON.d;
FESETENVD(fpenv);
return x;
}
else {
FEGETENVD(fpenv);
FESETENVD(kUP.f);
x = x + (long double)EPSILON.d;
FESETENVD(fpenv);
return x;
}
}
float nexttowardf(float x, long double y)
{
DblDbl ldu;
double yh, yt, yd, fpenv;
if ((long double)x == y)
return (float)y;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = y;
yh = ldu.d[0];
yt = ldu.d[1];
yd = yh + yt;
FESETENVD(fpenv);
return nextafterf(x, (float)yd);
}
double nexttoward(double x, long double y)
{
DblDbl ldu;
double yh, yt, yd, fpenv;
if ((long double)x == y)
return (double)y;
FEGETENVD(fpenv);
FESETENVD(0.0);
ldu.f = y;
yh = ldu.d[0];
yt = ldu.d[1];
yd = yh + yt;
FESETENVD(fpenv);
return nextafter(x, y);
}
long double nexttowardl(long double x, long double y)
{
if (x == y)
return y; else
return nextafterl( x, y );
}
#endif