Fixed aliasing bugs in TRUNC() by using the fdlibm macros for access

to doubles as bits.  fdlibm-1.1 had similar aliasing bugs, but these
were fixed by NetBSD or Cygnus before a modified version of fdlibm was
imported in 1994.  TRUNC() is only used by tgamma() and some
implementation-detail functions.  The aliasing bugs were detected by
compiling with gcc -O2 but don't seem to have broken tgamma() on i386's
or amd64's.  They broke my modified version of tgamma().

Moved the definition of TRUNC() to mathimpl.h so that it can be fixed
in one place, although the general version is even slower than necessary
because it has to operate on pointers to volatiles to handle its arg
sometimes being volatile.  Inefficiency of the fdlibm macros slows
down libm generally, and tgamma() is a relatively unimportant part of
libm.  The macros act as if on 32-bit words in memory, so they are
hard to optimize to direct actions on 64-bit double registers for
(non-i386) machines where this is possible.  The optimization is too
hard for gcc on amd64's, and declaring variables as volatile makes it
impossible.

3 files changed, 26 insertions, 12 deletions

diff --git a/lib/msun/bsdsrc/b_log.c b/lib/msun/bsdsrc/b_log.c
index 95d366d..b9f1473 100644
--- a/lib/msun/bsdsrc/b_log.c
+++ b/lib/msun/bsdsrc/b_log.c
@@ -77,9 +77,6 @@ __FBSDID("$FreeBSD$");
  *	+Inf	return +Inf
 */
 
-#define endian		(((*(int *) &one)) ? 1 : 0)
-#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
-
 #define N 128
 
 /* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
@@ -438,7 +435,7 @@ __log__D(x) double x;
 #endif
 {
 	int m, j;
-	double F, f, g, q, u, v, u2, one = 1.0;
+	double F, f, g, q, u, v, u2;
 	volatile double u1;
 	struct Double r;
 
diff --git a/lib/msun/bsdsrc/b_tgamma.c b/lib/msun/bsdsrc/b_tgamma.c
index babf326..a8e5755 100644
--- a/lib/msun/bsdsrc/b_tgamma.c
+++ b/lib/msun/bsdsrc/b_tgamma.c
@@ -123,20 +123,12 @@ static struct Double ratfun_gam(double, double);
 #define Pa7	-1.44705562421428915453880392761e-02
 
 static const double zero = 0., one = 1.0, tiny = 1e-300;
-static int endian;
-
-/*
- * TRUNC sets trailing bits in a floating-point number to zero.
- * is a temporary variable.
- */
-#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
 
 double
 tgamma(x)
 	double x;
 {
 	struct Double u;
-	endian = (*(int *) &one) ? 1 : 0;
 
 	if (x >= 6) {
 		if(x > 171.63)
diff --git a/lib/msun/bsdsrc/mathimpl.h b/lib/msun/bsdsrc/mathimpl.h
index 275b624..5c22e7d 100644
--- a/lib/msun/bsdsrc/mathimpl.h
+++ b/lib/msun/bsdsrc/mathimpl.h
@@ -34,9 +34,32 @@
  * $FreeBSD$
  */
 
+#ifndef _MATHIMPL_H_
+#define	_MATHIMPL_H_
+
 #include <sys/cdefs.h>
 #include <math.h>
 
+#include "../src/math_private.h"
+
+/*
+ * TRUNC() is a macro that sets the trailing 27 bits in the mantissa of an
+ * IEEE double variable to zero.  It must be expression-like for syntactic
+ * reasons, and we implement this expression using an inline function
+ * instead of a pure macro to avoid depending on the gcc feature of
+ * statement-expressions.
+ */
+#define	TRUNC(d)	(_b_trunc(&(d)))
+
+static __inline void
+_b_trunc(volatile double *_dp)
+{
+	uint32_t _lw;
+
+	GET_LOW_WORD(_lw, *_dp);
+	SET_LOW_WORD(*_dp, _lw & 0xf8000000);
+}
+
 /*
  * Functions internal to the math package, yet not static.
  */
@@ -45,3 +68,5 @@ extern double	__exp__E();
 struct Double {double a, b;};
 double __exp__D(double, double);
 struct Double __log__D(double);
+
+#endif /* !_MATHIMPL_H_ */