jyapayne/musl
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

long double cleanup, initial commit

Browse source 
new ldshape union, ld128 support is kept, code that used the old
ldshape union was rewritten (IEEEl2bits union of freebsd libm is
not touched yet)

ld80 __fpclassifyl no longer tries to handle invalid representation
This commit is contained in:
Szabolcs Nagy 2013-09-02 00:38:51 +00:00
commit af5f6d9556
8 changed files with 89 additions and 96 deletions
Download patch file Download diff file Expand all files Collapse all files

28
src/math/__fpclassifyl.c
View file

@ -3,28 +3,28 @@
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* invalid representations (top bit of u.i.m is wrong) are not handled */
int __fpclassifyl(long double x)
{
union ldshape u = { x };
int e = u.bits.exp;
if (!e) {
if (u.bits.m >> 63) return FP_NAN;
else if (u.bits.m) return FP_SUBNORMAL;
else return FP_ZERO;
}
union ldshape u = {x};
int e = u.i.se & 0x7fff;
if (!e)
return u.i.m ? FP_SUBNORMAL : FP_ZERO;
if (e == 0x7fff)
return u.bits.m & (uint64_t)-1>>1 ? FP_NAN : FP_INFINITE;
return u.bits.m & (uint64_t)1<<63 ? FP_NORMAL : FP_NAN;
return u.i.m << 1 ? FP_NAN : FP_INFINITE;
return FP_NORMAL;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
int __fpclassifyl(long double x)
{
union ldshape u = { x };
int e = u.bits.exp;
union ldshape u = {x};
int e = u.i.se & 0x7fff;
if (!e)
return u.bits.mlo | u.bits.mhi ? FP_SUBNORMAL : FP_ZERO;
if (e == 0x7fff)
return u.bits.mlo | u.bits.mhi ? FP_NAN : FP_INFINITE;
return u.i2.lo | u.i2.hi ? FP_SUBNORMAL : FP_ZERO;
if (e == 0x7fff) {
u.i.se = 0;
return u.i2.lo | u.i2.hi ? FP_NAN : FP_INFINITE;
}
return FP_NORMAL;
}
#endif

4
src/math/__signbitl.c
View file

@ -1,11 +1,9 @@
#include "libm.h"
// FIXME: should be a macro
#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
int __signbitl(long double x)
{
union ldshape u = {x};
return u.bits.sign;
return u.i.se >> 15;
}
#endif

6
src/math/copysignl.c
View file

@ -9,8 +9,8 @@ long double copysignl(long double x, long double y)
long double copysignl(long double x, long double y)
{
union ldshape ux = {x}, uy = {y};
ux.bits.sign = uy.bits.sign;
return ux.value;
ux.i.se &= 0x7fff;
ux.i.se |= uy.i.se & 0x8000;
return ux.f;
}
#endif

4
src/math/fabsl.c
View file

@ -9,7 +9,7 @@ long double fabsl(long double x)
{
union ldshape u = {x};
u.bits.sign = 0;
return u.value;
u.i.se &= 0x7fff;
return u.f;
}
#endif

8
src/math/ilogbl.c
View file

@ -10,12 +10,12 @@ int ilogbl(long double x)
int ilogbl(long double x)
{
union ldshape u = {x};
uint64_t m = u.bits.m;
int e = u.bits.exp;
uint64_t m = u.i.m;
int e = u.i.se & 0x7fff;
if (!e) {
if (m == 0) {
FORCE_EVAL(0/0.0f);
FORCE_EVAL(x/x);
return FP_ILOGB0;
}
/* subnormal x */
@ -25,7 +25,7 @@ int ilogbl(long double x)
if (e == 0x7fff) {
FORCE_EVAL(0/0.0f);
/* in ld80 msb is set in inf */
return m & (uint64_t)-1>>1 ? FP_ILOGBNAN : INT_MAX;
return m << 1 ? FP_ILOGBNAN : INT_MAX;
}
return e - 0x3fff;
}

81
src/math/nextafterl.c
View file

@ -6,7 +6,6 @@ long double nextafterl(long double x, long double y)
return nextafter(x, y);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
#define MSB ((uint64_t)1<<63)
long double nextafterl(long double x, long double y)
{
union ldshape ux, uy;
@ -15,32 +14,32 @@ long double nextafterl(long double x, long double y)
return x + y;
if (x == y)
return y;
ux.value = x;
ux.f = x;
if (x == 0) {
uy.value = y;
ux.bits.m = 1;
ux.bits.sign = uy.bits.sign;
} else if (x < y ^ ux.bits.sign) {
ux.bits.m++;
if ((ux.bits.m & ~MSB) == 0) {
ux.bits.m = MSB;
ux.bits.exp++;
uy.f = y;
ux.i.m = 1;
ux.i.se = uy.i.se & 0x8000;
} else if ((x < y) == !(ux.i.se & 0x8000)) {
ux.i.m++;
if (ux.i.m << 1 == 0) {
ux.i.m = 1ULL << 63;
ux.i.se++;
}
} else {
if ((ux.bits.m & ~MSB) == 0) {
ux.bits.exp--;
if (ux.bits.exp)
ux.bits.m = 0;
if (ux.i.m << 1 == 0) {
ux.i.se--;
if (ux.i.se)
ux.i.m = 0;
}
ux.bits.m--;
ux.i.m--;
}
/* raise overflow if ux.value is infinite and x is finite */
if (ux.bits.exp == 0x7fff)
/* raise overflow if ux is infinite and x is finite */
if ((ux.i.se & 0x7fff) == 0x7fff)
return x + x;
/* raise underflow if ux.value is subnormal or zero */
if (ux.bits.exp == 0)
FORCE_EVAL(x*x + ux.value*ux.value);
return ux.value;
/* raise underflow if ux is subnormal or zero */
if ((ux.i.se & 0x7fff) == 0)
FORCE_EVAL(x*x + ux.f*ux.f);
return ux.f;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
long double nextafterl(long double x, long double y)
@ -51,32 +50,26 @@ long double nextafterl(long double x, long double y)
return x + y;
if (x == y)
return y;
ux.value = x;
ux.f = x;
if (x == 0) {
uy.value = y;
ux.bits.mlo = 1;
ux.bits.sign = uy.bits.sign;
} else if (x < y ^ ux.bits.sign) {
ux.bits.mlo++;
if (ux.bits.mlo == 0) {
ux.bits.mhi++;
if (ux.bits.mhi == 0)
ux.bits.exp++;
}
uy.f = y;
ux.i.lo = 1;
ux.i.se = uy.i.se & 0x8000;
} else if ((x < y) == !(ux.i.se & 0x8000)) {
ux.i2.lo++;
if (ux.i2.lo == 0)
ux.i2.hi++;
} else {
if (ux.bits.mlo == 0) {
if (ux.bits.mhi == 0)
ux.bits.exp--;
ux.bits.mhi--;
}
ux.bits.mlo--;
if (ux.i2.lo == 0)
ux.i2.hi--;
ux.i2.lo--;
}
/* raise overflow if ux.value is infinite and x is finite */
if (ux.bits.exp == 0x7fff)
/* raise overflow if ux is infinite and x is finite */
if ((ux.i.se & 0x7fff) == 0x7fff)
return x + x;
/* raise underflow if ux.value is subnormal or zero */
if (ux.bits.exp == 0)
FORCE_EVAL(x*x + ux.value*ux.value);
return ux.value;
/* raise underflow if ux is subnormal or zero */
if ((ux.i.se & 0x7fff) == 0)
FORCE_EVAL(x*x + ux.f*ux.f);
return ux.f;
}
#endif