 @@ -15,33 +15,34 @@
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  */
 #include <sys/cdefs.h>
-__RCSID("$NetBSD: n_exp2.c,v 1.2 2014/03/12 19:42:18 martin Exp $");
+__RCSID("$NetBSD: n_exp2.c,v 1.3 2014/03/23 15:26:47 martin Exp $");
 #ifdef __FBSDID
 __FBSDID("$FreeBSD: src/lib/msun/src/s_exp2.c,v 1.7 2008/02/22 02:27:34 das Exp $");
 #endif
 #include <stdint.h>
 #include <float.h>
 #include <string.h>
 #include "math.h"
 #define	TBLBITS	8
 #define	TBLSIZE	(1 << TBLBITS)
 static const double
     huge     = 0x1p126,
     redux    = 0x1.8p52 / TBLSIZE,
     P1	     = 0x1.62e42fefa39efp-1,
     P2	     = 0x1.ebfbdff82c575p-3,
     P3	     = 0x1.c6b08d704a0a6p-5,
     P4	     = 0x1.3b2ab88f70400p-7,
 @@ -334,59 +335,61 @@ static const double tbl[TBLSIZE * 2] = {
  *   Note that the range of i is +-TBLSIZE/2, so we actually index the tables
  *   by i0 = i + TBLSIZE/2.  For cache efficiency, exp2t[] and eps[] are
  *   virtual tables, interleaved in the real table tbl[].
+ *
  *   This method is due to Gal, with many details due to Gal and Bachelis:
+ *
  *	Gal, S. and Bachelis, B.  An Accurate Elementary Mathematical Library
  *	for the IEEE Floating Point Standard.  TOMS 17(1), 26-46 (1991).
  */
 double
 exp2(double x)
+{
 	double r, t, twopk, z;
-	uint32_t hx, ix, i0;
+	uint32_t hx, ix, i0, temp;
 	int k, big;
 	/* Filter out exceptional cases. */
-	hx = ((uint32_t*)&x)[0];
+	memcpy(&hx, &x, sizeof(hx));
 	ix = hx & 0x7fffffff;		/* high word of |x| */
 	if(ix >= 0x40900000) {			/* |x| >= 1024 */
 		if(x >= 0x1.0p10)
 			return (huge * huge); /* overflow */
 		if(x <= -0x1.0ccp10)
 			return (twom1000 * twom1000); /* underflow */
 	} else if (ix < 0x3c900000) {		/* |x| < 0x1p-54 */
 		return (1.0 + x);
+	}
 	/* Reduce x, computing z, i0, and k. */
 	(*(volatile double*)&t) = x + redux;
 	i0 = ((uint32_t*)&t)[1];
 	i0 += TBLSIZE / 2;
 	k = (i0 >> TBLBITS) << 20;
 	i0 = (i0 & (TBLSIZE - 1)) << 1;
 	t -= redux;
 	z = x - t;
 	/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
 	t = tbl[i0];		/* exp2t[i0] */
 	z -= tbl[i0 + 1];	/* eps[i0]   */
 	big = k >= -1021 << 20;
 	if (big) {
-		((uint32_t*)&twopk)[0] = 0x3ff00000+k;
+		temp = 0x3ff00000+k;
-		((uint32_t*)&twopk)[1] = 0;
+		twopk = 0.0;
 		memcpy(&twopk, &temp, sizeof(temp));
 	} else {
-		((uint32_t*)&twopk)[0] = 0x3ff00000+k + (1000 << 20);
+		temp =  0x3ff00000+k + (1000 << 20);
-		((uint32_t*)&twopk)[1] = 0;
+		twopk = 0.0;
 		memcpy(&twopk, &temp, sizeof(temp));
+	}
 	r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
 	/* Scale by 2**(k>>20). */
 	if (big) {
 		/*if (k == 1024 << 20)
 			return (r * 2.0 * 0x1p1023); */
 		return (r * twopk);
 	} else {
 		return (r * twopk * twom1000);
+	}
+}

 @@ -15,33 +15,34 @@
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  */
 #include <sys/cdefs.h>
-__RCSID("$NetBSD: n_exp2f.c,v 1.1 2014/03/06 10:55:57 martin Exp $");
+__RCSID("$NetBSD: n_exp2f.c,v 1.2 2014/03/23 15:26:47 martin Exp $");
 #ifdef __FBSDID
 __FBSDID("$FreeBSD: src/lib/msun/src/s_exp2f.c,v 1.9 2008/02/22 02:27:34 das Exp $");
 #endif
 #include <stdint.h>
 #include <float.h>
 #include <string.h>
 #include "math.h"
 #define	TBLBITS	4
 #define	TBLSIZE	(1 << TBLBITS)
 static const float
     huge    = 0x1p100f,
     redux   = 0x1.8p23f / TBLSIZE,
     P1	    = 0x1.62e430p-1f,
     P2	    = 0x1.ebfbe0p-3f,
     P3	    = 0x1.c6b348p-5f,
     P4	    = 0x1.3b2c9cp-7f;
 @@ -89,46 +90,47 @@ static const double exp2ft[TBLSIZE] = {
  *   roundoff error insignificant and simplifies the scaling step.
+ *
  *   This method is due to Tang, but I do not use his suggested parameters:
+ *
  *	Tang, P.  Table-driven Implementation of the Exponential Function
  *	in IEEE Floating-Point Arithmetic.  TOMS 15(2), 144-157 (1989).
  */
 float
 exp2f(float x)
+{
 	double tv, twopk, u, z;
 	float t;
 	uint32_t hx, ix, i0;
-	int32_t k;
+	int32_t k, temp;
 	/* Filter out exceptional cases. */
-	hx = *((uint32_t*)&x);
+	memcpy(&hx, &x, sizeof(hx));
 	ix = hx & 0x7fffffff;		/* high word of |x| */
 	if(ix >= 0x43000000) {			/* |x| >= 128 */
 		if(x >= 0x1.0p7f)
 			return (huge * huge);	/* overflow */
 		if(x <= -0x1.2cp7f)
 			return (twom100 * twom100); /* underflow */
 	} else if (ix <= 0x33000000) {		/* |x| <= 0x1p-25 */
 		return (1.0f + x);
+	}
 	/* Reduce x, computing z, i0, and k. */
-	*((volatile float*)&t) = x + redux;
+	i0 = x + redux;
-	i0 = *((uint32_t*)&t);
+	memcpy(&t, &i0, sizeof(t));
 	i0 += TBLSIZE / 2;
 	k = (i0 >> TBLBITS) << 20;
 	i0 &= TBLSIZE - 1;
 	t -= redux;
 	z = x - t;
-	((uint32_t*)&twopk)[0] = 0x3ff00000+k;
+	temp = 0x3ff00000+k;
-	((uint32_t*)&twopk)[1] = 0;
+	twopk = 0.0;
 	memcpy(&twopk, &temp, sizeof(temp));
 	/* Compute r = exp2(y) = exp2ft[i0] * p(z). */
 	tv = exp2ft[i0];
 	u = tv * z;
 	tv = tv + u * (P1 + z * P2) + u * (z * z) * (P3 + z * P4);
 	/* Scale by 2**(k>>20). */
 	return (tv * twopk);
+}