Add __exp*_finite optimizations and rewrite some wrappers.
2011-10-15 Ulrich Drepper <drepper@gmail.com>
+ * math/Versions [libm] (GLIBC_2.15): Add __exp_finite, __expf_finite,
+ __expl_finite.
+ * math/bits/math-finite.h: Add entries for exp.
+ * math/e_expl.c: Add __*_finite alias.
+ * sysdeps/i386/fpu/e_exp.S: Likewise.
+ * sysdeps/i386/fpu/e_expf.S: Likewise.
+ * sysdeps/i386/fpu/e_expl.c: Likewise.
+ * sysdeps/ieee754/dbl-64/e_exp.c: Likewise.
+ * sysdeps/ieee754/flt-32/e_expf.c: Likewise.
+ * sysdeps/ieee754/ldbl-128/e_expl.c: Likewise.
+ * sysdeps/ieee754/ldbl-128ibm/e_expl.c: Likewise.
+ * sysdeps/ieee754/dbl-64/w_exp.c: Complete rewrite.
+ * sysdeps/ieee754/flt-32/w_expf.c: Likewise.
+ * sysdeps/ieee754/ldbl-96/w_expl.c: Likewise.
+
* sysdeps/i386/i686/fpu/e_logf.S: No need for the fyl2xp1 use, fyl2x
is sufficient, at least on modern CPUs.
__sinh_finite; __sinhf_finite; __sinhl_finite;
__sqrt_finite; __sqrtf_finite; __sqrtl_finite;
__gamma_r_finite; __gammaf_r_finite; __gammal_r_finite;
+ __exp_finite; __expf_finite; __expl_finite;
}
}
extern long double coshl (long double) __asm__ ("__coshl_finite");
#endif
+/* exp. */
+extern double exp (double) __asm__ ("__exp_finite");
+extern float expf (float) __asm__ ("__expf_finite");
+#ifdef __MATH_DECLARE_LDOUBLE
+extern long double expl (long double) __asm__ ("__expl_finite");
+#endif
+
#ifdef __USE_GNU
/* exp10. */
extern double exp10 (double) __asm__ ("__exp10_finite");
__set_errno (ENOSYS);
return 0.0;
}
+strong_alias (__ieee754_expl, __expl_finite)
stub_warning (expl)
#include <stub-tag.h>
#include <machine/asm.h>
-RCSID("$NetBSD: e_exp.S,v 1.7 1996/07/03 17:31:28 jtc Exp $")
/* e^x = 2^(x * log2(e)) */
ENTRY(__ieee754_exp)
fldz /* Set result to 0. */
2: ret
END (__ieee754_exp)
+
+
+ENTRY(__exp_finite)
+ fldl2e
+ fmull 4(%esp) /* x * log2(e) */
+ fld %st
+ frndint /* int(x * log2(e)) */
+ fsubr %st,%st(1) /* fract(x * log2(e)) */
+ fxch
+ f2xm1 /* 2^(fract(x * log2(e))) - 1 */
+ fld1
+ faddp /* 2^(fract(x * log2(e))) */
+ fscale /* e^x */
+ fstp %st(1)
+ ret
+END(__exp_finite)
#include <machine/asm.h>
-RCSID("$NetBSD: $")
/* e^x = 2^(x * log2(e)) */
ENTRY(__ieee754_expf)
fldz /* Set result to 0. */
2: ret
END (__ieee754_expf)
+
+
+ENTRY(__expf_finite)
+ fldl2e
+ fmuls 4(%esp) /* x * log2(e) */
+ fld %st
+ frndint /* int(x * log2(e)) */
+ fsubr %st,%st(1) /* fract(x * log2(e)) */
+ fxch
+ f2xm1 /* 2^(fract(x * log2(e))) - 1 */
+ fld1
+ faddp /* 2^(fract(x * log2(e))) */
+ fscale /* e^x */
+ fstp %st(1)
+ ret
+END(__expf_finite)
"fld1\n\t" /* 4 1.0 */
"faddp\n\t" /* 3 2^(fract(x * log2(e))) */
"fstp %%st(1)\n\t" /* 2 */
- "fscale\n\t" /* 2 scale factor is st(1); e^x */
+ "fscale\n\t" /* 2 scale factor is st(1); e^x */
"fstp %%st(1)\n\t" /* 1 */
"fstp %%st(1)\n\t" /* 0 */
"jmp 2f\n\t"
: "=t" (res) : "0" (x), "m" (c0), "m" (c1) : "ax", "dx");
return res;
}
+strong_alias (__ieee754_expl, __expl_finite)
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
else return __slowexp(x);
}
}
+strong_alias (__ieee754_exp, __exp_finite)
/************************************************************************/
/* Compute e^(x+xx)(Double-Length number) .The routine also receive */
-/* @(#)w_exp.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.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: w_exp.c,v 1.6 1995/05/10 20:48:51 jtc Exp $";
-#endif
+/* Copyright (C) 2011 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+ 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.
-/*
- * wrapper exp(x)
- */
+ 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.
-#include "math.h"
-#include "math_private.h"
+ 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>
+#include <math_private.h>
-#ifdef __STDC__
static const double
-#else
-static double
-#endif
o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
u_threshold= -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */
-#ifdef __STDC__
- double __exp(double x) /* wrapper exp */
-#else
- double __exp(x) /* wrapper exp */
- double x;
-#endif
+
+/* wrapper exp */
+double
+__exp (double x)
{
-#ifdef _IEEE_LIBM
- return __ieee754_exp(x);
-#else
- double z;
- z = __ieee754_exp(x);
- if(_LIB_VERSION == _IEEE_) return z;
- if(__finite(x)) {
- if(x>o_threshold)
- return __kernel_standard(x,x,6); /* exp overflow */
- else if(x<u_threshold)
- return __kernel_standard(x,x,7); /* exp underflow */
- }
- return z;
-#endif
+ if (__builtin_expect (x > o_threshold, 0))
+ {
+ if (_LIB_VERSION != _IEEE_)
+ return __kernel_standard_f (x, x, 6);
+ }
+ else if (__builtin_expect (x < u_threshold, 0))
+ {
+ if (_LIB_VERSION != _IEEE_)
+ return __kernel_standard_f (x, x, 7);
+ }
+
+ return __ieee754_exp (x);
}
hidden_def (__exp)
weak_alias (__exp, exp)
/* Single-precision floating point e^x.
- Copyright (C) 1997, 1998, 2005, 2006 Free Software Foundation, Inc.
+ Copyright (C) 1997, 1998, 2005, 2006, 2011 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
Then e^x is approximated as
e^x = 2^n ( e^(t/512 + delta[t])
- + ( e^(t/512 + delta[t])
- * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
+ + ( e^(t/512 + delta[t])
+ * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
where
- p(x) is a polynomial approximating e(x)-1;
/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
return TWO127*x;
}
+strong_alias (__ieee754_expf, __expf_finite)
-/* w_expf.c -- float version of w_exp.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: w_expf.c,v 1.3 1995/05/10 20:48:53 jtc Exp $";
-#endif
-
-/*
- * wrapper expf(x)
- */
-
-#include "math.h"
-#include "math_private.h"
-
-#ifdef __STDC__
+/* Copyright (C) 2011 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+ 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>
+#include <math_private.h>
+
static const float
-#else
-static float
-#endif
o_threshold= 8.8722831726e+01, /* 0x42b17217 */
u_threshold= -1.0397208405e+02; /* 0xc2cff1b5 */
-#ifdef __STDC__
- float __expf(float x) /* wrapper expf */
-#else
- float __expf(x) /* wrapper expf */
- float x;
-#endif
+
+/* wrapper expf */
+float
+__expf (float x)
{
-#ifdef _IEEE_LIBM
- return __ieee754_expf(x);
-#else
- float z;
- z = __ieee754_expf(x);
- if(_LIB_VERSION == _IEEE_) return z;
- if(__finitef(x)) {
- if(x>o_threshold)
- /* exp overflow */
- return (float)__kernel_standard((double)x,(double)x,106);
- else if(x<u_threshold)
- /* exp underflow */
- return (float)__kernel_standard((double)x,(double)x,107);
- }
- return z;
-#endif
+ if (__builtin_expect (x > o_threshold, 0))
+ {
+ if (_LIB_VERSION != _IEEE_)
+ return __kernel_standard_f (x, x, 106);
+ }
+ else if (__builtin_expect (x < u_threshold, 0))
+ {
+ if (_LIB_VERSION != _IEEE_)
+ return __kernel_standard_f (x, x, 107);
+ }
+
+ return __ieee754_expf (x);
}
hidden_def (__expf)
weak_alias (__expf, expf)
/* Quad-precision floating point e^x.
- Copyright (C) 1999 Free Software Foundation, Inc.
+ Copyright (C) 1999, 2011 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Jakub Jelinek <jj@ultra.linux.cz>
Partly based on double-precision code
/* Smallest integer x for which e^x overflows. */
#define himark C[0]
11356.523406294143949491931077970765L,
-
+
/* Largest integer x for which e^x underflows. */
#define lomark C[1]
-11433.4627433362978788372438434526231L,
/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
return TWO16383*x;
}
+strong_alias (__ieee754_expl, __expl_finite)
/* Quad-precision floating point e^x.
- Copyright (C) 1999,2004,2006, 2008 Free Software Foundation, Inc.
+ Copyright (C) 1999,2004,2006, 2008, 2011 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Jakub Jelinek <jj@ultra.linux.cz>
Partly based on double-precision code
/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
return TWO1023*x;
}
+strong_alias (__ieee754_expl, __expl_finite)
-/* w_expl.c -- long double version of w_exp.c.
- * Conversion to long double by Ulrich Drepper,
- * Cygnus Support, drepper@cygnus.com.
- */
-
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: $";
-#endif
-
-/*
- * wrapper expl(x)
- */
-
-#include "math.h"
-#include "math_private.h"
-
-#ifdef __STDC__
+/* Copyright (C) 2011 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+ 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>
+#include <math_private.h>
+
static const long double
-#else
-static long double
-#endif
o_threshold= 1.135652340629414394949193107797076489134e4,
/* 0x400C, 0xB17217F7, 0xD1CF79AC */
u_threshold= -1.140019167866942050398521670162263001513e4;
/* 0x400C, 0xB220C447, 0x69C201E8 */
-#ifdef __STDC__
- long double __expl(long double x) /* wrapper exp */
-#else
- long double __expl(x) /* wrapper exp */
- long double x;
-#endif
+
+/* wrapper expl */
+long double
+__expl (long double x)
{
-#ifdef _IEEE_LIBM
- return __ieee754_expl(x);
-#else
- long double z;
- z = __ieee754_expl(x);
- if(_LIB_VERSION == _IEEE_) return z;
- if(__finitel(x)) {
- if(x>o_threshold)
- return __kernel_standard(x,x,206); /* exp overflow */
- else if(x<u_threshold)
- return __kernel_standard(x,x,207); /* exp underflow */
- }
- return z;
-#endif
+ if (__builtin_expect (x > o_threshold, 0))
+ {
+ if (_LIB_VERSION != _IEEE_)
+ return __kernel_standard (x, x, 206);
+ }
+ else if (__builtin_expect (x < u_threshold, 0))
+ {
+ if (_LIB_VERSION != _IEEE_)
+ return __kernel_standard (x, x, 207);
+ }
+
+ return __ieee754_expl (x);
}
hidden_def (__expl)
weak_alias (__expl, expl)
projects / external / glibc.git / commitdiff