Get the ld80 routines from OpenBSD to build on mac and linux.

Bump version number and SO major version, since we have
introduced new long double APIs.

Make.inc

@@ -3,8 +3,8 @@
 OS := $(shell uname)
 # Do not forget to bump SOMINOR when changing VERSION,
 # and SOMAJOR when breaking ABI in a backward-incompatible way
-VERSION = 0.4
-SOMAJOR = 1
+VERSION = 0.5
+SOMAJOR = 2
 SOMINOR = 0
 DESTDIR =
 prefix = /usr/local
@@ -64,7 +64,7 @@ endif
 
 ifeq ($(ARCH),x86_64)
 override ARCH := amd64
-endif 
+endif
 
 ifneq (,$(findstring MINGW,$(OS)))
 override OS=WINNT

ld80/Make.files

@@ -1,10 +1,12 @@
 $(CUR_SRCS) += 	invtrig.c \
-            e_acoshl.c     e_hypotl.c     e_powl.c       k_tanl.c       s_exp2l.c      s_nanl.c \
+            e_acoshl.c     e_powl.c       k_tanl.c       s_exp2l.c \
             e_atanhl.c     e_lgammal.c    e_sinhl.c      s_asinhl.c     s_expm1l.c \
-            e_coshl.c      e_log10l.c     e_tgammal.c    s_floorl.c     s_nextafterl.c \
-            e_expl.c       e_log2l.c      k_cosl.c       s_ceill.c      s_log1pl.c     s_tanhl.c \
-            e_fmodl.c      e_logl.c       k_sinl.c       s_erfl.c       s_modfl.c      s_truncl.c
-#           s_remquol.c
+            e_coshl.c      e_log10l.c     e_tgammal.c \
+            e_expl.c       e_log2l.c      k_cosl.c       s_log1pl.c     s_tanhl.c \
+            e_logl.c       k_sinl.c       s_erfl.c
+
+#           s_remquol.c    e_fmodl.c      s_truncl.c
+#           e_hypotl.c     s_floorl.c     s_nextafterl.c s_ceill.c      s_modfl.c
 
 ifneq ($(OS), WINNT)
 $(CUR_SRCS) += s_nanl.c

ld80/e_acoshl.c

@@ -24,7 +24,7 @@
  *	acoshl(NaN) is NaN without signal.
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_atanhl.c

@@ -28,7 +28,7 @@
  *
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_coshl.c

@@ -31,7 +31,7 @@
  *	only coshl(0)=1 is exact for finite x.
  */
 
-#include "math.h"
+#include "openlibm.h"
 #include "math_private.h"
 
 static const long double one = 1.0, half=0.5, huge = 1.0e4900L;

ld80/e_expl.c

@@ -72,7 +72,7 @@
 
 /*	Exponential function	*/
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_hypotl.c

@@ -42,7 +42,7 @@
  * 	than 1 ulps (units in the last place)
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_lgammal.c

@@ -86,9 +86,10 @@
  *
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
+extern int signgam;
 
 static const long double
   half = 0.5L,

ld80/e_log10l.c

@@ -63,7 +63,7 @@
  * log domain:       x < 0; returns MINLOG
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_log2l.c

@@ -63,7 +63,7 @@
  * log domain:       x < 0; returns NAN
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_logl.c

@@ -63,7 +63,7 @@
  * log domain:       x < 0; returns NAN
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_powl.c

@@ -73,7 +73,7 @@
  */
 
 #include <float.h>
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_rem_pio2l.h

@@ -6,7 +6,7 @@
  *
  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  *
@@ -17,8 +17,8 @@
 //__FBSDID("$FreeBSD: src/lib/msun/ld80/e_rem_pio2l.h,v 1.3 2011/06/18 13:56:33 benl Exp $");
 
 /* ld80 version of __ieee754_rem_pio2l(x,y)
- * 
- * return the remainder of x rem pi/2 in y[0]+y[1] 
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
  * use __kernel_rem_pio2()
  */
 
@@ -26,7 +26,7 @@
 
 #include "openlibm.h"
 #include "math_private.h"
-#include "fpmath.h"
+#include "openlibm.h"
 
 #define	BIAS	(LDBL_MAX_EXP - 1)
 
@@ -102,24 +102,24 @@ __ieee754_rem_pio2l(long double x, long double *y)
 		union IEEEl2bits u2;
 	        int ex1;
 	        j  = ex;
-	        y[0] = r-w; 
+	        y[0] = r-w;
 		u2.e = y[0];
 		ex1 = u2.xbits.expsign & 0x7fff;
 	        i = j-ex1;
 	        if(i>22) {  /* 2nd iteration needed, good to 141 */
 		    t  = r;
-		    w  = fn*pio2_2;	
+		    w  = fn*pio2_2;
 		    r  = t-w;
-		    w  = fn*pio2_2t-((t-r)-w);	
+		    w  = fn*pio2_2t-((t-r)-w);
 		    y[0] = r-w;
 		    u2.e = y[0];
 		    ex1 = u2.xbits.expsign & 0x7fff;
 		    i = j-ex1;
 		    if(i>61) {	/* 3rd iteration need, 180 bits acc */
 		    	t  = r;	/* will cover all possible cases */
-		    	w  = fn*pio2_3;	
+		    	w  = fn*pio2_3;
 		    	r  = t-w;
-		    	w  = fn*pio2_3t-((t-r)-w);	
+		    	w  = fn*pio2_3t-((t-r)-w);
 		    	y[0] = r-w;
 		    }
 		}
@@ -127,7 +127,7 @@ __ieee754_rem_pio2l(long double x, long double *y)
 	    y[1] = (r-y[0])-w;
 	    return n;
 	}
-    /* 
+    /*
      * all other (large) arguments
      */
 	if(ex==0x7fff) {		/* x is inf or NaN */

ld80/e_sinhl.c

@@ -28,7 +28,7 @@
  *	only sinhl(0)=0 is exact for finite x.
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/e_tgammal.c

@@ -58,9 +58,10 @@
  */
 
 #include <float.h>
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
+extern int signgam;
 
 /*
 tgamma(x+2)  = tgamma(x+2) P(x)/Q(x)

ld80/invtrig.h

@@ -28,7 +28,7 @@
 
 #include <float.h>
 
-#include "fpmath.h"
+#include <openlibm.h>
 
 #define	BIAS		(LDBL_MAX_EXP - 1)
 #define	MANH_SIZE	LDBL_MANH_SIZE

ld80/s_asinhl.c

@@ -21,7 +21,7 @@
  *		  := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_ceill.c

@@ -19,7 +19,7 @@
  *	Inexact flag raised if x not equal to ceil(x).
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_erfl.c

@@ -99,7 +99,7 @@
  */
 
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_exp2l.c

@@ -32,7 +32,7 @@
 
 #include "amd64/bsd_ieeefp.h"
 
-#include "fpmath.h"
+#include "openlibm.h"
 #include "openlibm.h"
 #include "math_private.h"
 

ld80/s_expm1l.c

@@ -57,7 +57,7 @@
  *
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 static const long double MAXLOGL = 1.1356523406294143949492E4L;
 

ld80/s_floorl.c

@@ -19,7 +19,7 @@
  *	Inexact flag raised if x not equal to floor(x).
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_log1pl.c

@@ -59,7 +59,7 @@
  * log domain:       x-1 < 0; returns NAN
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_modfl.c

@@ -20,7 +20,7 @@
  *	No exception.
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_nanl.c

@@ -28,7 +28,7 @@
 
 #include "openlibm.h"
 
-#include "fpmath.h"
+#include "openlibm.h"
 #include "math_private.h"
 
 DLLEXPORT long double

ld80/s_nextafterl.c

@@ -17,7 +17,7 @@
  *   Special cases:
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 
@@ -87,4 +87,4 @@ nextafterl(long double x, long double y)
 	return x;
 }
 
-__strong_alias(nexttowardl, nextafterl);
+//__strong_alias(nexttowardl, nextafterl);

ld80/s_nexttoward.c

@@ -17,7 +17,7 @@
  *   Special cases:
  */
 
-#include <math.h>
+#include <openlibm.h>
 #include <float.h>
 
 #include "math_private.h"

ld80/s_nexttowardf.c

@@ -10,7 +10,7 @@
  * ====================================================
  */
 
-#include <math.h>
+#include <openlibm.h>
 #include <float.h>
 
 #include "math_private.h"

ld80/s_remquol.c

@@ -14,7 +14,7 @@
 #include <machine/ieee.h>
 
 #include <float.h>
-#include <math.h>
+#include <openlibm.h>
 #include <stdint.h>
 
 #include "math_private.h"

ld80/s_tanhl.c

@@ -34,7 +34,7 @@
  *	only tanhl(0)=0 is exact for finite argument.
  */
 
-#include <math.h>
+#include <openlibm.h>
 
 #include "math_private.h"
 

ld80/s_truncl.c

@@ -21,10 +21,10 @@
  */
 
 #include <sys/types.h>
-#include <machine/ieee.h>
+//#include <machine/ieee.h>
 
 #include <float.h>
-#include <math.h>
+#include "openlibm.h"
 #include <stdint.h>
 
 #include "math_private.h"

src/Make.files

@@ -50,7 +50,7 @@ $(CUR_SRCS) +=	e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
 	s_frexpl.c s_logbl.c s_nexttoward.c \
 	s_remquol.c \
 	s_sinl.c s_sincosl.c s_tanl.c s_truncl.c w_cabsl.c \
-	s_nextafterl.c s_rintl.c s_scalbnl.c
+	s_nextafterl.c s_rintl.c s_scalbnl.c polevll.c
 #	s_cbrtl.c
 endif
 

src/math_private.h

@@ -22,6 +22,7 @@
 #include "fpmath.h"
 #include <complex.h>
 #include <stdint.h>
+#include "math_private_openbsd.h"
 
 /*
  * The original fdlibm code used statements like:
@@ -228,7 +229,7 @@ do {								\
  * Common routine to process the arguments to nan(), nanf(), and nanl().
  */
 void _scan_nan(u_int32_t *__words, int __num_words, const char *__s);
- 
+
 #ifdef __GNUCLIKE_ASM
 
 /* Asm versions of some functions. */

src/math_private_openbsd.h

@@ -0,0 +1,218 @@
+/*	$OpenBSD: math_private.h,v 1.17 2014/06/02 19:31:17 kettenis Exp $	*/
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ */
+
+#ifndef _MATH_PRIVATE_OPENBSD_H_
+#define _MATH_PRIVATE_OPENBSD_H_
+
+#if _IEEE_WORD_ORDER == _BIG_ENDIAN
+
+typedef union
+{
+  long double value;
+  struct {
+    u_int32_t mswhi;
+    u_int32_t mswlo;
+    u_int32_t lswhi;
+    u_int32_t lswlo;
+  } parts32;
+  struct {
+    u_int64_t msw;
+    u_int64_t lsw;
+  } parts64;
+} ieee_quad_shape_type;
+
+#endif
+
+#if _IEEE_WORD_ORDER == _LITTLE_ENDIAN
+
+typedef union
+{
+  long double value;
+  struct {
+    u_int32_t lswlo;
+    u_int32_t lswhi;
+    u_int32_t mswlo;
+    u_int32_t mswhi;
+  } parts32;
+  struct {
+    u_int64_t lsw;
+    u_int64_t msw;
+  } parts64;
+} ieee_quad_shape_type;
+
+#endif
+
+/* Get two 64 bit ints from a long double.  */
+
+#define GET_LDOUBLE_WORDS64(ix0,ix1,d)				\
+do {								\
+  ieee_quad_shape_type qw_u;					\
+  qw_u.value = (d);						\
+  (ix0) = qw_u.parts64.msw;					\
+  (ix1) = qw_u.parts64.lsw;					\
+} while (0)
+
+/* Set a long double from two 64 bit ints.  */
+
+#define SET_LDOUBLE_WORDS64(d,ix0,ix1)				\
+do {								\
+  ieee_quad_shape_type qw_u;					\
+  qw_u.parts64.msw = (ix0);					\
+  qw_u.parts64.lsw = (ix1);					\
+  (d) = qw_u.value;						\
+} while (0)
+
+/* Get the more significant 64 bits of a long double mantissa.  */
+
+#define GET_LDOUBLE_MSW64(v,d)					\
+do {								\
+  ieee_quad_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  (v) = sh_u.parts64.msw;					\
+} while (0)
+
+/* Set the more significant 64 bits of a long double mantissa from an int.  */
+
+#define SET_LDOUBLE_MSW64(d,v)					\
+do {								\
+  ieee_quad_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  sh_u.parts64.msw = (v);					\
+  (d) = sh_u.value;						\
+} while (0)
+
+/* Get the least significant 64 bits of a long double mantissa.  */
+
+#define GET_LDOUBLE_LSW64(v,d)					\
+do {								\
+  ieee_quad_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  (v) = sh_u.parts64.lsw;					\
+} while (0)
+
+/* A union which permits us to convert between a long double and
+   three 32 bit ints.  */
+
+#if _IEEE_WORD_ORDER == _BIG_ENDIAN
+
+typedef union
+{
+  long double value;
+  struct {
+#ifdef __LP64__
+    int padh:32;
+#endif
+    int exp:16;
+    int padl:16;
+    u_int32_t msw;
+    u_int32_t lsw;
+  } parts;
+} ieee_extended_shape_type;
+
+#endif
+
+#if _IEEE_WORD_ORDER == _LITTLE_ENDIAN
+
+typedef union
+{
+  long double value;
+  struct {
+    u_int32_t lsw;
+    u_int32_t msw;
+    int exp:16;
+    int padl:16;
+#ifdef __LP64__
+    int padh:32;
+#endif
+  } parts;
+} ieee_extended_shape_type;
+
+#endif
+
+/* Get three 32 bit ints from a double.  */
+
+#define GET_LDOUBLE_WORDS(se,ix0,ix1,d)				\
+do {								\
+  ieee_extended_shape_type ew_u;				\
+  ew_u.value = (d);						\
+  (se) = ew_u.parts.exp;					\
+  (ix0) = ew_u.parts.msw;					\
+  (ix1) = ew_u.parts.lsw;					\
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+
+#define SET_LDOUBLE_WORDS(d,se,ix0,ix1)				\
+do {								\
+  ieee_extended_shape_type iw_u;				\
+  iw_u.parts.exp = (se);					\
+  iw_u.parts.msw = (ix0);					\
+  iw_u.parts.lsw = (ix1);					\
+  (d) = iw_u.value;						\
+} while (0)
+
+/* Get the more significant 32 bits of a long double mantissa.  */
+
+#define GET_LDOUBLE_MSW(v,d)					\
+do {								\
+  ieee_extended_shape_type sh_u;				\
+  sh_u.value = (d);						\
+  (v) = sh_u.parts.msw;						\
+} while (0)
+
+/* Set the more significant 32 bits of a long double mantissa from an int.  */
+
+#define SET_LDOUBLE_MSW(d,v)					\
+do {								\
+  ieee_extended_shape_type sh_u;				\
+  sh_u.value = (d);						\
+  sh_u.parts.msw = (v);						\
+  (d) = sh_u.value;						\
+} while (0)
+
+/* Get int from the exponent of a long double.  */
+
+#define GET_LDOUBLE_EXP(se,d)					\
+do {								\
+  ieee_extended_shape_type ge_u;				\
+  ge_u.value = (d);						\
+  (se) = ge_u.parts.exp;					\
+} while (0)
+
+/* Set exponent of a long double from an int.  */
+
+#define SET_LDOUBLE_EXP(d,se)					\
+do {								\
+  ieee_extended_shape_type se_u;				\
+  se_u.value = (d);						\
+  se_u.parts.exp = (se);					\
+  (d) = se_u.value;						\
+} while (0)
+
+/*
+ * Common routine to process the arguments to nan(), nanf(), and nanl().
+ */
+void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
+
+/*
+ * Functions internal to the math package, yet not static.
+ */
+double __exp__D(double, double);
+struct Double __log__D(double);
+long double __p1evll(long double, void *, int);
+long double __polevll(long double, void *, int);
+
+#endif /* _MATH_PRIVATE_OPENBSD_H_ */

src/polevll.c

@@ -0,0 +1,104 @@
+/*	$OpenBSD: polevll.c,v 1.2 2013/11/12 20:35:09 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <[email protected]>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							polevll.c
+ *							p1evll.c
+ *
+ *	Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * long double x, y, coef[N+1], polevl[];
+ *
+ * y = polevll( x, coef, N );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates polynomial of degree N:
+ *
+ *                     2          N
+ * y  =  C  + C x + C x  +...+ C x
+ *        0    1     2          N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C  , ..., coef[N] = C  .
+ *            N                   0
+ *
+ *  The function p1evll() assumes that coef[N] = 1.0 and is
+ * omitted from the array.  Its calling arguments are
+ * otherwise the same as polevll().
+ *
+ *
+ * SPEED:
+ *
+ * In the interest of speed, there are no checks for out
+ * of bounds arithmetic.  This routine is used by most of
+ * the functions in the library.  Depending on available
+ * equipment features, the user may wish to rewrite the
+ * program in microcode or assembly language.
+ *
+ */
+
+#include "openlibm.h"
+
+#include "math_private.h"
+
+/*
+ * Polynomial evaluator:
+ *  P[0] x^n  +  P[1] x^(n-1)  +  ...  +  P[n]
+ */
+long double
+__polevll(long double x, void *PP, int n)
+{
+	long double y;
+	long double *P;
+
+	P = (long double *)PP;
+	y = *P++;
+	do {
+		y = y * x + *P++;
+	} while (--n);
+
+	return (y);
+}
+
+/*
+ * Polynomial evaluator:
+ *  x^n  +  P[0] x^(n-1)  +  P[1] x^(n-2)  +  ...  +  P[n]
+ */
+long double
+__p1evll(long double x, void *PP, int n)
+{
+	long double y;
+	long double *P;
+
+	P = (long double *)PP;
+	n -= 1;
+	y = x + *P++;
+	do {
+		y = y * x + *P++;
+	} while (--n);
+
+	return (y);
+}