fix

git-svn-id: http://svn2.nishi.boats/svn/milsko/trunk@555 b9cfdab3-6d41-4d17-bbe4-086880011989

src/math/default.c

@@ -173,8 +173,7 @@ static MwLLMathVTable table_u16 = {
     default_shiftLeft_u16,
     default_equal_u16,
     default_greaterThen_u16,
-    NULL
-};
+    NULL};
 
 static void default_add_u32(MwLLVec* a, MwLLVec* b, MwLLVec* out) {
 	out->un.u32.a = a->un.u32.a + b->un.u32.a;
@@ -225,8 +224,7 @@ static MwLLMathVTable table_u32 = {
     default_shiftLeft_u32,
     default_equal_u32,
     default_greaterThen_u32,
-    NULL
-};
+    NULL};
 
 static void default_add_i8(MwLLVec* a, MwLLVec* b, MwLLVec* out) {
 	out->un.i8.a = a->un.i8.a + b->un.i8.a;
@@ -311,8 +309,7 @@ static MwLLMathVTable table_i8 = {
     default_shiftLeft_u8,
     default_equal_i8,
     default_greaterThen_i8,
-    NULL
-};
+    NULL};
 static void default_add_i16(MwLLVec* a, MwLLVec* b, MwLLVec* out) {
 	out->un.i16.a = a->un.i16.a + b->un.i16.a;
 	out->un.i16.b = a->un.i16.b + b->un.i16.b;
@@ -368,8 +365,7 @@ static MwLLMathVTable table_i16 = {
     default_shiftLeft_u16,
     default_equal_i16,
     default_greaterThen_i16,
-    NULL
-};
+    NULL};
 
 static void default_add_i32(MwLLVec* a, MwLLVec* b, MwLLVec* out) {
 	out->un.i32.a = a->un.i32.a + b->un.i32.a;
@@ -412,8 +408,7 @@ static MwLLMathVTable table_i32 = {
     default_shiftLeft_u32,
     default_equal_i32,
     default_greaterThen_i32,
-    NULL
-};
+    NULL};
 
 static MwLLMathVTable* defaultMultiTable[_MwLLVecType_Max] = {
     &table_u8,	// _MwLLVecTypeU8x8

src/math/nbsd_copysign.c

@@ -0,0 +1,31 @@
+/* @(#)s_copysign.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* $Id$ */
+
+#include <Mw/BaseTypes.h>
+#include "math_internal.h"
+
+/*
+ * copysign(double x, double y)
+ * copysign(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+double
+nbsd_copysign(double x, double y) {
+	MwU32 hx, hy;
+	GET_HIGH_WORD(hx, x);
+	GET_HIGH_WORD(hy, y);
+	SET_HIGH_WORD(x, (hx & 0x7fffffff) | (hy & 0x80000000));
+	return x;
+}

src/math/nbsd_endian.c

@@ -1,5 +1,6 @@
 /* $Id$ */
-#include "nbsd_math.h"
+#include <Mw/BaseTypes.h>
+#include "math_internal.h"
 
 static char endian = 0;
 

src/math/nbsd_math.h

@@ -95,6 +95,9 @@ typedef union _ieee_double_shape_type {
 	} while(0)
 
 double nbsd_pow(double a, double b);
+double nbsd_scalbn(double x, int n);
+double nbsd_scalbln(double x, long n);
+double nbsd_copysign(double x, double y);
 char   nbsd_endian(void);
 
 #endif

src/math/nbsd_pow.c

@@ -133,7 +133,7 @@ nbsd_pow(double x, double y) {
 			k = (iy >> 20) - 0x3ff; /* exponent */
 			if(k > 20) {
 				j = ly >> (52 - k);
-				if((uint32_t)(j << (52 - k)) == ly) yisint = 2 - (j & 1);
+				if((MwU32)(j << (52 - k)) == ly) yisint = 2 - (j & 1);
 			} else if(ly == 0) {
 				j = iy >> (20 - k);
 				if((j << (20 - k)) == iy) yisint = 2 - (j & 1);
@@ -314,7 +314,7 @@ nbsd_pow(double x, double y) {
 	z  = one - (r - z);
 	GET_HIGH_WORD(j, z);
 	j += (n << 20);
-	if((j >> 20) <= 0) z = scalbn(z, n); /* subnormal output */
+	if((j >> 20) <= 0) z = nbsd_scalbn(z, n); /* subnormal output */
 	else
 		SET_HIGH_WORD(z, j);
 	return s * z;

src/math/nbsd_scalbn.c

@@ -0,0 +1,65 @@
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* $Id$ */
+
+#include <Mw/BaseTypes.h>
+#include "math_internal.h"
+
+/*
+ * scalbn (double x, int n)
+ * scalbn(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+static const double
+    two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+    twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+    hugev  = 1.0e+300,
+    tinyv  = 1.0e-300;
+
+double
+nbsd_scalbn(double x, int n) {
+	return nbsd_scalbln(x, n);
+}
+
+double
+nbsd_scalbln(double x, long n) {
+	MwI32 k, hx, lx;
+	EXTRACT_WORDS(hx, lx, x);
+	k = ((MwU32)hx & 0x7ff00000) >> 20;		    /* extract exponent */
+	if(k == 0) {					    /* 0 or subnormal x */
+		if((lx | (hx & 0x7fffffff)) == 0) return x; /* +-0 */
+		x *= two54;
+		GET_HIGH_WORD(hx, x);
+		k = (((MwU32)hx & 0x7ff00000) >> 20) - 54;
+		if(n < -50000) return tinyv * x; /*underflow*/
+	}
+	if(k == 0x7ff) return x + x; /* NaN or Inf */
+	k = k + n;
+	if(k > 0x7fe) return hugev * nbsd_copysign(hugev, x); /* overflow  */
+	if(k > 0)					      /* normal result */
+	{
+		SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
+		return x;
+	}
+	if(k <= -54) {
+		if(n > 50000)					/* in case integer overflow in n+k */
+			return hugev * nbsd_copysign(hugev, x); /*overflow*/
+		else
+			return tinyv * nbsd_copysign(tinyv, x); /*underflow*/
+	}
+	k += 54; /* subnormal result */
+	SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
+	return x * twom54;
+}