- add some more math functions (patch from gentoo/solar)

10 files changed, 313 insertions, 1 deletions

diff --git a/libm/Makefile.in b/libm/Makefile.in
index 01c63d15a..0d06ec718 100644
--- a/libm/Makefile.in
+++ b/libm/Makefile.in
@@ -71,7 +71,9 @@ libm_CSRC := \
 	w_log.c w_log10.c w_pow.c w_remainder.c w_scalb.c w_sinh.c \
 	w_sqrt.c nan.c carg.c s_llrint.c \
 	s_fpclassify.c s_fpclassifyf.c s_signbit.c s_signbitf.c \
-	s_isnan.c s_isnanf.c s_isinf.c s_isinff.c
+	s_isnan.c s_isnanf.c s_isinf.c s_isinff.c \
+	s_fdim.c s_fma.c s_fmax.c s_fmin.c s_nearbyint.c \
+	s_remquo.c s_scalbln.c w_exp2.c w_tgamma.c
 FL_MOBJ := \
 	acosf.o acoshf.o asinf.o asinhf.o atan2f.o atanf.o atanhf.o cbrtf.o \
 	ceilf.o copysignf.o cosf.o coshf.o erfcf.o erff.o exp2f.o expf.o \
diff --git a/libm/s_fdim.c b/libm/s_fdim.c
new file mode 100644
index 000000000..6fdf3c2d6
--- /dev/null
+++ b/libm/s_fdim.c
@@ -0,0 +1,25 @@
+/* Copyright (C) 2002 by  Red Hat, Incorporated. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software
+ * is freely granted, provided that this notice is preserved.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(fdim)
+#ifdef __STDC__
+	double fdim(double x, double y)
+#else
+	double fdim(x,y)
+	double x;
+	double y;
+#endif
+{
+  int c = __fpclassify(x);
+  if (c == FP_NAN || c == FP_INFINITE)
+    return HUGE_VAL;
+
+  return x > y ? x - y : 0.0;
+}
+libm_hidden_def(fdim)
diff --git a/libm/s_fma.c b/libm/s_fma.c
new file mode 100644
index 000000000..e5ff5a722
--- /dev/null
+++ b/libm/s_fma.c
@@ -0,0 +1,33 @@
+/* Compute x * y + z as ternary operation.
+   Copyright (C) 1997, 2001 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+
+double
+__fma (double x, double y, double z)
+{
+  return (x * y) + z;
+}
+weak_alias (__fma, fma)
+
+#ifdef NO_LONG_DOUBLE
+strong_alias (__fma, __fmal)
+weak_alias (__fmal, fmal)
+#endif
diff --git a/libm/s_fmax.c b/libm/s_fmax.c
new file mode 100644
index 000000000..61ec12d81
--- /dev/null
+++ b/libm/s_fmax.c
@@ -0,0 +1,26 @@
+/* Copyright (C) 2002 by  Red Hat, Incorporated. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software
+ * is freely granted, provided that this notice is preserved.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(fmax)
+#ifdef __STDC__
+	double fmax(double x, double y)
+#else
+	double fmax(x,y)
+	double x;
+	double y;
+#endif
+{
+  if (__fpclassify(x) == FP_NAN)
+    return x;
+  if (__fpclassify(y) == FP_NAN)
+    return y;
+
+  return x > y ? x : y;
+}
+libm_hidden_def(fmax)
diff --git a/libm/s_fmin.c b/libm/s_fmin.c
new file mode 100644
index 000000000..e30808d85
--- /dev/null
+++ b/libm/s_fmin.c
@@ -0,0 +1,26 @@
+/* Copyright (C) 2002 by  Red Hat, Incorporated. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software
+ * is freely granted, provided that this notice is preserved.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(fmin)
+#ifdef __STDC__
+	double fmin(double x, double y)
+#else
+	double fmin(x,y)
+	double x;
+	double y;
+#endif
+{
+  if (__fpclassify(x) == FP_NAN)
+    return x;
+  if (__fpclassify(y) == FP_NAN)
+    return y;
+
+  return x < y ? x : y;
+}
+libm_hidden_def(fmin)
diff --git a/libm/s_nearbyint.c b/libm/s_nearbyint.c
new file mode 100644
index 000000000..1bd8d4054
--- /dev/null
+++ b/libm/s_nearbyint.c
@@ -0,0 +1,25 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(nearbyint)
+#ifdef __STDC__
+	double nearbyint(double x)
+#else
+	double nearbyint(x)
+	double x;
+#endif
+{
+  return rint(x);
+}
+libm_hidden_def(nearbyint)
diff --git a/libm/s_remquo.c b/libm/s_remquo.c
new file mode 100644
index 000000000..18853e211
--- /dev/null
+++ b/libm/s_remquo.c
@@ -0,0 +1,38 @@
+/* Copyright (C) 2002 by  Red Hat, Incorporated. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software
+ * is freely granted, provided that this notice is preserved.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(remquo)
+#ifdef __STDC__
+	double remquo(double x, double y, int *quo)	/* wrapper remquo */
+#else
+	double remquo(x,y,quo)			/* wrapper remquo */
+	double x,y;
+        int *quo;
+#endif
+{
+        int signx, signy, signres;
+        int mswx;
+        int mswy;
+        double x_over_y;
+
+        GET_HIGH_WORD(mswx, x);
+        GET_HIGH_WORD(mswy, y);
+
+        signx = (mswx & 0x80000000) >> 31;
+        signy = (mswy & 0x80000000) >> 31;
+
+        signres = (signx ^ signy) ? -1 : 1;
+
+        x_over_y = fabs(x / y);
+
+        *quo = signres * (lrint(x_over_y) & 0x7f);
+
+        return remainder(x,y);
+}
+libm_hidden_def(remquo)
diff --git a/libm/s_scalbln.c b/libm/s_scalbln.c
new file mode 100644
index 000000000..a780f3710
--- /dev/null
+++ b/libm/s_scalbln.c
@@ -0,0 +1,63 @@
+/* @(#)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.
+ * ====================================================
+ */
+
+/*
+ * 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.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(scalbln)
+#ifdef __STDC__
+static const double
+#else
+static double
+#endif
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
+
+#ifdef __STDC__
+	double scalbln (double x, long int n)
+#else
+	double scalbln (x,n)
+	double x; long int n;
+#endif
+{
+	int32_t k,hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+        k = (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 = ((hx&0x7ff00000)>>20) - 54;
+	    }
+        if (k==0x7ff) return x+x;		/* NaN or Inf */
+        k = k+n;
+        if (n> 50000 || k >  0x7fe)
+	  return huge*copysign(huge,x); /* overflow  */
+	if (n< -50000) return tiny*copysign(tiny,x); /*underflow*/
+        if (k > 0)				/* normal result */
+	    {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
+        if (k <= -54)
+	  return tiny*copysign(tiny,x);	/*underflow*/
+        k += 54;				/* subnormal result */
+	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
+        return x*twom54;
+}
+libm_hidden_def(scalbln)
diff --git a/libm/w_exp2.c b/libm/w_exp2.c
new file mode 100644
index 000000000..20a9cf35b
--- /dev/null
+++ b/libm/w_exp2.c
@@ -0,0 +1,27 @@
+
+/* @(#)w_exp2.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.
+ * ====================================================
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(exp2)
+#ifdef __STDC__
+	double exp2(double x)
+#else
+	double exp2(x)
+	double x;
+#endif
+{
+  return pow(2.0, x);
+}
+libm_hidden_def(exp2)
diff --git a/libm/w_tgamma.c b/libm/w_tgamma.c
new file mode 100644
index 000000000..62b3c6275
--- /dev/null
+++ b/libm/w_tgamma.c
@@ -0,0 +1,47 @@
+/* @(#)w_gamma.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.
+ * ====================================================
+ */
+
+/* double gamma(double x)
+ * Return  the logarithm of the Gamma function of x or the Gamma function of x,
+ * depending on the library mode.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(tgamma)
+#ifdef __STDC__
+	double tgamma(double x)
+#else
+	double tgamma(x)
+	double x;
+#endif
+{
+        double y;
+	int local_signgam;
+	y = __ieee754_gamma_r(x,&local_signgam);
+	if (local_signgam < 0) y = -y;
+#ifdef _IEEE_LIBM
+	return y;
+#else
+	if(_LIB_VERSION == _IEEE_) return y;
+
+	if(!finite(y)&&finite(x)) {
+	  if(floor(x)==x&&x<=0.0)
+	    return __kernel_standard(x,x,41); /* tgamma pole */
+	  else
+	    return __kernel_standard(x,x,40); /* tgamma overflow */
+	}
+	return y;
+#endif
+}
+libm_hidden_def(tgamma)