merge from upstream, updated to 0.7.5

Make.inc

@@ -10,41 +10,40 @@ VERSION = 0.7.0
 SOMAJOR = 3
 SOMINOR = 0
 DESTDIR =
-prefix = /usr/local
-bindir = $(prefix)/bin
-libdir = $(prefix)/lib
-includedir = $(prefix)/include
+prefix ?= /usr/local
+bindir ?= $(prefix)/bin
+libdir ?= $(prefix)/lib
+includedir ?= $(prefix)/include
 
 ifeq ($(OS), FreeBSD)
-pkgconfigdir = $(prefix)/libdata/pkgconfig
+pkgconfigdir ?= $(prefix)/libdata/pkgconfig
 else
-pkgconfigdir = $(libdir)/pkgconfig
+pkgconfigdir ?= $(libdir)/pkgconfig
 endif
 
-USEGCC = 1
-USECLANG = 0
+USEGCC ?= 1
+USECLANG ?= 0
 
 ifneq (,$(findstring $(OS),Darwin FreeBSD OpenBSD))
-USEGCC = 0
-USECLANG = 1
+USEGCC ?= 0
+USECLANG ?= 1
 endif
 
-AR = $(TOOLPREFIX)ar
-
 ifeq ($(ARCH),wasm32)
-CC = clang-8
-USEGCC = 0
-CFLAGS_add += -fno-builtin -fno-strict-aliasing
+USECLANG = 1
+TOOLPREFIX = llvm-
 endif
 
+AR ?= $(TOOLPREFIX)ar
+
 ifeq ($(USECLANG),1)
-USEGCC = 0
+USEGCC ?= 0
 CC = clang
 CFLAGS_add += -fno-builtin -fno-strict-aliasing
 endif
 
 ifeq ($(USEGCC),1)
-CC = $(TOOLPREFIX)gcc
+CC ?= $(TOOLPREFIX)gcc
 CFLAGS_add += -fno-gnu89-inline -fno-builtin
 endif
 

aarch64/fenv.c

@@ -26,7 +26,6 @@
  * $FreeBSD: src/lib/msun/arm/fenv.c,v 1.3 2011/10/16 05:37:56 das Exp $
  */
 
-#define	__fenv_static
 #include <openlibm_fenv.h>
 
 #ifdef __GNUC_GNU_INLINE__

amd64/fenv.c

@@ -29,7 +29,7 @@
 #include "bsd_fpu.h"
 #include "math_private.h"
 
-#ifndef OPENLIBM_USE_HOST_FENV_H
+#ifdef _WIN32
 #define __fenv_static OLM_DLLEXPORT
 #endif
 #include <openlibm_fenv.h>

include/openlibm_math.h

@@ -304,6 +304,7 @@ OLM_DLLEXPORT double	trunc(double);
  * BSD math library entry points
  */
 #if __BSD_VISIBLE
+OLM_DLLEXPORT int	isinff(float) __pure2;
 OLM_DLLEXPORT int	isnanf(float) __pure2;
 
 /*

powerpc/fenv.c

@@ -26,7 +26,6 @@
  * $FreeBSD$
  */
 
-#define	__fenv_static
 #include <openlibm_fenv.h>
 
 #ifdef __GNUC_GNU_INLINE__

src/e_hypotl.c

@@ -82,7 +82,7 @@ hypotl(long double x, long double y)
 	        man_t manh, manl;
 		GET_LDBL_MAN(manh,manl,b);
 		if((manh|manl)==0) return a;
-		t1=0;
+		t1=1;
 		SET_HIGH_WORD(t1,ESW(MAX_EXP-2));	/* t1=2^(MAX_EXP-2) */
 		b *= t1;
 		a *= t1;

src/e_powf.c

@@ -25,6 +25,9 @@ bp[] = {1.0, 1.5,},
 dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
 dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
 zero    =  0.0,
+half    =  0.5,
+qrtr    =  0.25,
+thrd    =  3.33333343e-01, /* 0x3eaaaaab */
 one	=  1.0,
 two	=  2.0,
 two24	=  16777216.0,	/* 0x4b800000 */
@@ -74,7 +77,7 @@ __ieee754_powf(float x, float y)
     /* y!=zero: result is NaN if either arg is NaN */
 	if(ix > 0x7f800000 ||
 	   iy > 0x7f800000)
-		return (x+0.0F)+(y+0.0F);
+	    return nan_mix(x, y);
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
@@ -103,15 +106,10 @@ __ieee754_powf(float x, float y)
 	if(iy==0x3f800000) {	/* y is  +-1 */
 	    if(hy<0) return one/x; else return x;
 	}
-        if(hy==0x40000000) return x*x;   /* y is  2 */
-        if(hy==0x40400000) return x*x*x; /* y is  3 */
-        if(hy==0x40800000) {             /* y is  4 */
-            u = x*x;
-            return u*u;
-        }
-        if(hy==0x3f000000) {             /* y is  0.5 */
+	if(hy==0x40000000) return x*x; /* y is  2 */
+	if(hy==0x3f000000) {	/* y is  0.5 */
 	    if(hx>=0)	/* x >= +0 */
-                return __ieee754_sqrtf(x);
+	    return __ieee754_sqrtf(x);
 	}
 
 	ax   = fabsf(x);
@@ -139,12 +137,12 @@ __ieee754_powf(float x, float y)
     /* |y| is huge */
 	if(iy>0x4d000000) { /* if |y| > 2**27 */
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+	    if(ix<0x3f7ffff7) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
 	    if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
 	/* now |1-x| is tiny <= 2**-20, suffice to compute
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
 	    t = ax-1;		/* t has 20 trailing zeros */
-	    w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
+	    w = (t*t)*(half-t*(thrd-t*qrtr));
 	    u = ivln2_h*t;	/* ivln2_h has 16 sig. bits */
 	    v = t*ivln2_l-w*ivln2;
 	    t1 = u+v;
@@ -183,10 +181,10 @@ __ieee754_powf(float x, float y)
 	    r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
 	    r += s_l*(s_h+s);
 	    s2  = s_h*s_h;
-	    t_h = (float)3.0+s2+r;
+	    t_h = 3+s2+r;
 	    GET_FLOAT_WORD(is,t_h);
 	    SET_FLOAT_WORD(t_h,is&0xfffff000);
-	    t_l = r-((t_h-(float)3.0)-s2);
+	    t_l = r-((t_h-3)-s2);
 	/* u+v = s*(1+...) */
 	    u = s_h*t_h;
 	    v = s_l*t_h+t_l*s;
@@ -198,7 +196,7 @@ __ieee754_powf(float x, float y)
 	    z_h = cp_h*p_h;		/* cp_h+cp_l = 2/(3*log2) */
 	    z_l = cp_l*p_h+p_l*cp+dp_l[k];
 	/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
-	    t = (float)n;
+	    t = n;
 	    t1 = (((z_h+z_l)+dp_h[k])+t);
 	    GET_FLOAT_WORD(is,t1);
 	    SET_FLOAT_WORD(t1,is&0xfffff000);

src/math_private.h

@@ -203,10 +203,9 @@ do {								\
 } while (0)
 
 
-//VBS
+#ifndef __FreeBSD__
 #define	STRICT_ASSIGN(type, lval, rval)	((lval) = (rval))
-
-/* VBS
+#else
 #ifdef FLT_EVAL_METHOD
 // Attempt to get strict C99 semantics for assignment with non-C99 compilers.
 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
@@ -215,7 +214,7 @@ do {								\
 #define	STRICT_ASSIGN(type, lval, rval) do {	\
 	volatile type __lval;			\
 						\
-	if (sizeof(type) >= sizeof(double))	\
+	if (sizeof(type) >= sizeof(long double))	\
 		(lval) = (rval);		\
 	else {					\
 		__lval = (rval);		\
@@ -224,13 +223,31 @@ do {								\
 } while (0)
 #endif
 #endif
-*/
+#endif
 
 /*
  * Common routine to process the arguments to nan(), nanf(), and nanl().
  */
 void __scan_nan(u_int32_t *__words, int __num_words, const char *__s);
 
+/*
+ * Mix 1 or 2 NaNs.  First add 0 to each arg.  This normally just turns
+ * signaling NaNs into quiet NaNs by setting a quiet bit.  We do this
+ * because we want to never return a signaling NaN, and also because we
+ * don't want the quiet bit to affect the result.  Then mix the converted
+ * args using addition.  The result is typically the arg whose mantissa
+ * bits (considered as in integer) are largest.
+ *
+ * Technical complications: the result in bits might depend on the precision
+ * and/or on compiler optimizations, especially when different register sets
+ * are used for different precisions.  Try to make the result not depend on
+ * at least the precision by always doing the main mixing step in long double
+ * precision.  Try to reduce dependencies on optimizations by adding the
+ * the 0's in different precisions (unless everything is in long double
+ * precision).
+ */
+#define	nan_mix(x, y)	(((x) + 0.0L) + ((y) + 0))
+
 #ifdef __GNUCLIKE_ASM
 
 /* Asm versions of some functions. */

test/test-211.c

@@ -0,0 +1,12 @@
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+
+int
+main()
+{
+  float x = 0xd.65874p-4f;
+  float y = 4.0f;
+  float z = powf (x, y);
+  assert(z==0x1.f74424p-2);
+}