fmaf.c revision 1ec68d336ba97cd53f46053ac10401d16014d075
/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License (the "License").
* You may not use this file except in compliance with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#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) {
union {
unsigned i[2];
double d;
unsigned u, s;
/*
* 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;
/*
* if the sum xy + z will be exact, just compute it and cast the
* result to float
*/
}
/*
* collapse the tail of the smaller summand into a "sticky bit"
* so that the sum can be computed without error
*/
u = xy.i[1];
if (u & (s - 1))
u |= s;
u = xy.i[0];
u |= s;
xy.i[0] = u & ~(s - 1);
xy.i[1] = 0;
} else {
/* collapse all of xy into a single bit */
xy.i[1] = 0;
}
} else {
u = zz.i[1];
if (u & (s - 1))
u |= s;
u = zz.i[0];
u |= s;
zz.i[0] = u & ~(s - 1);
zz.i[1] = 0;
} else {
/* collapse all of zz into a single bit */
zz.i[1] = 0;
}
}
}
#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) {
union {
unsigned i[NI];
long double e;
/* set rounding precision to 64 bits */
/*
* 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;
/*
* if the sum xy + z will be exact, just compute it and cast the
* result to float
*/
goto cont;
}
/*
* collapse the tail of the smaller summand into a "sticky bit"
* so that the sum can be computed without error
*/
u = xy.i[0];
if (u & (s - 1))
u |= s;
xy.i[0] = u & ~(s - 1);
u = xy.i[1];
if ((u & (s - 1)) | xy.i[0])
u |= s;
xy.i[0] = 0;
} else {
/* collapse all of xy into a single bit */
xy.i[0] = 0;
}
} else {
u = zz.i[1];
if ((u & (s - 1)) | zz.i[0])
u |= s;
zz.i[0] = 0;
} else {
/* collapse all of zz into a single bit */
zz.i[0] = 0;
}
}
cont:
/* restore the rounding precision */
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) {
/*
* 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:
*
* long double t;
*
* <set rp mode to round to 64 bit precision>
* t = x * y;
* <restore rp mode>
* return t + z;
*
* 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);
}
#endif
#else
#endif