Print this page
11210 libm should be cstyle(1ONBLD) clean
@@ -20,10 +20,11 @@
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
@@ -32,58 +33,63 @@
#include "libm.h"
#include "fma.h"
#include "fenv_inlines.h"
#if defined(__sparc)
-
/*
* fmaf for SPARC: 32-bit single precision, big-endian
*/
float
-__fmaf(float x, float y, float z) {
+__fmaf(float x, float y, float z)
+{
union {
unsigned i[2];
double d;
} xy, zz;
+
unsigned u, s;
int exy, ez;
/*
* the following operations can only raise the invalid exception,
* and then only if either x*y is of the form Inf*0 or one of x,
* y, or z is a signaling NaN
*/
- xy.d = (double) x * y;
- zz.d = (double) z;
+ xy.d = (double)x * y;
+ zz.d = (double)z;
/*
* if the sum xy + z will be exact, just compute it and cast the
* result to float
*/
exy = (xy.i[0] >> 20) & 0x7ff;
ez = (zz.i[0] >> 20) & 0x7ff;
+
if ((ez - exy <= 4 && exy - ez <= 28) || exy == 0x7ff || exy == 0 ||
- ez == 0x7ff || ez == 0) {
- return ((float) (xy.d + zz.d));
- }
+ ez == 0x7ff || ez == 0)
+ return ((float)(xy.d + zz.d));
/*
* collapse the tail of the smaller summand into a "sticky bit"
* so that the sum can be computed without error
*/
if (ez > exy) {
if (ez - exy < 31) {
u = xy.i[1];
s = 2 << (ez - exy);
+
if (u & (s - 1))
u |= s;
+
xy.i[1] = u & ~(s - 1);
} else if (ez - exy < 51) {
u = xy.i[0];
s = 1 << (ez - exy - 31);
+
if ((u & (s - 1)) | xy.i[1])
u |= s;
+
xy.i[0] = u & ~(s - 1);
xy.i[1] = 0;
} else {
/* collapse all of xy into a single bit */
xy.i[0] = (xy.i[0] & 0x80000000) | ((ez - 51) << 20);
@@ -91,47 +97,51 @@
}
} else {
if (exy - ez < 31) {
u = zz.i[1];
s = 2 << (exy - ez);
+
if (u & (s - 1))
u |= s;
+
zz.i[1] = u & ~(s - 1);
} else if (exy - ez < 51) {
u = zz.i[0];
s = 1 << (exy - ez - 31);
+
if ((u & (s - 1)) | zz.i[1])
u |= s;
+
zz.i[0] = u & ~(s - 1);
zz.i[1] = 0;
} else {
/* collapse all of zz into a single bit */
zz.i[0] = (zz.i[0] & 0x80000000) | ((exy - 51) << 20);
zz.i[1] = 0;
}
}
- return ((float) (xy.d + zz.d));
+ return ((float)(xy.d + zz.d));
}
-
#elif defined(__x86)
-
#if defined(__amd64)
#define NI 4
#else
#define NI 3
#endif
/*
* fmaf for x86: 32-bit single precision, little-endian
*/
float
-__fmaf(float x, float y, float z) {
+__fmaf(float x, float y, float z)
+{
union {
unsigned i[NI];
long double e;
} xy, zz;
+
unsigned u, s, cwsw, oldcwsw;
int exy, ez;
/* set rounding precision to 64 bits */
__fenv_getcwsw(&oldcwsw);
@@ -141,40 +151,44 @@
/*
* the following operations can only raise the invalid exception,
* and then only if either x*y is of the form Inf*0 or one of x,
* y, or z is a signaling NaN
*/
- xy.e = (long double) x * y;
- zz.e = (long double) z;
+ xy.e = (long double)x * y;
+ zz.e = (long double)z;
/*
* if the sum xy + z will be exact, just compute it and cast the
* result to float
*/
exy = xy.i[2] & 0x7fff;
ez = zz.i[2] & 0x7fff;
+
if ((ez - exy <= 15 && exy - ez <= 39) || exy == 0x7fff || exy == 0 ||
- ez == 0x7fff || ez == 0) {
+ ez == 0x7fff || ez == 0)
goto cont;
- }
/*
* collapse the tail of the smaller summand into a "sticky bit"
* so that the sum can be computed without error
*/
if (ez > exy) {
if (ez - exy < 31) {
u = xy.i[0];
s = 2 << (ez - exy);
+
if (u & (s - 1))
u |= s;
+
xy.i[0] = u & ~(s - 1);
} else if (ez - exy < 62) {
u = xy.i[1];
s = 1 << (ez - exy - 31);
+
if ((u & (s - 1)) | xy.i[0])
u |= s;
+
xy.i[1] = u & ~(s - 1);
xy.i[0] = 0;
} else {
/* collapse all of xy into a single bit */
xy.i[0] = 0;
@@ -183,12 +197,14 @@
}
} else {
if (exy - ez < 62) {
u = zz.i[1];
s = 1 << (exy - ez - 31);
+
if ((u & (s - 1)) | zz.i[0])
u |= s;
+
zz.i[1] = u & ~(s - 1);
zz.i[0] = 0;
} else {
/* collapse all of zz into a single bit */
zz.i[0] = 0;
@@ -203,20 +219,21 @@
/* restore the rounding precision */
__fenv_getcwsw(&cwsw);
cwsw = (cwsw & 0xfcffffff) | (oldcwsw & 0x03000000);
__fenv_setcwsw(&cwsw);
- return ((float) xy.e);
+ return ((float)xy.e);
}
#if 0
/*
* another fmaf for x86: assumes return value will be left in
* long double (80-bit double extended) precision
*/
long double
-__fmaf(float x, float y, float z) {
+__fmaf(float x, float y, float z)
+{
/*
* Note: This implementation assumes the rounding precision mode
* is set to the default, rounding to 64 bit precision. If this
* routine must work in non-default rounding precision modes, do
* the following instead:
@@ -230,12 +247,11 @@
*
* Note that the code to change rounding precision must not alter
* the exception masks or flags, since the product x * y may raise
* an invalid operation exception.
*/
- return ((long double) x * y + z);
+ return ((long double)x * y + z);
}
#endif
-
#else
#error Unknown architecture
#endif