diff options
Diffstat (limited to 'newlib/libm/common')
| -rw-r--r-- | newlib/libm/common/Makefile.am | 4 | ||||
| -rw-r--r-- | newlib/libm/common/Makefile.in | 13 | ||||
| -rw-r--r-- | newlib/libm/common/fdlibm.h | 2 | ||||
| -rw-r--r-- | newlib/libm/common/hypotl.c | 55 | ||||
| -rw-r--r-- | newlib/libm/common/sl_finite.c | 25 | ||||
| -rw-r--r-- | newlib/libm/common/sqrtl.c | 155 |
6 files changed, 244 insertions, 10 deletions
diff --git a/newlib/libm/common/Makefile.am b/newlib/libm/common/Makefile.am index b370b2e09..d5e0ef959 100644 --- a/newlib/libm/common/Makefile.am +++ b/newlib/libm/common/Makefile.am @@ -31,8 +31,8 @@ lsrc = atanl.c cosl.c sinl.c tanl.c tanhl.c frexpl.c modfl.c ceill.c fabsl.c \ scalbnl.c exp2l.c scalblnl.c tgammal.c nearbyintl.c lrintl.c llrintl.c \ roundl.c lroundl.c llroundl.c truncl.c remquol.c fdiml.c fmaxl.c fminl.c \ fmal.c acoshl.c atanhl.c remainderl.c lgammal.c erfl.c erfcl.c \ - logbl.c nexttowardf.c nexttoward.c nexttowardl.c log2l.c - + logbl.c nexttowardf.c nexttoward.c nexttowardl.c log2l.c \ + sl_finite.c libcommon_la_LDFLAGS = -Xcompiler -nostdlib diff --git a/newlib/libm/common/Makefile.in b/newlib/libm/common/Makefile.in index af3c88d71..ae5949a48 100644 --- a/newlib/libm/common/Makefile.in +++ b/newlib/libm/common/Makefile.in @@ -139,7 +139,7 @@ am__objects_3 = lib_a-atanl.$(OBJEXT) lib_a-cosl.$(OBJEXT) \ lib_a-erfl.$(OBJEXT) lib_a-erfcl.$(OBJEXT) \ lib_a-logbl.$(OBJEXT) lib_a-nexttowardf.$(OBJEXT) \ lib_a-nexttoward.$(OBJEXT) lib_a-nexttowardl.$(OBJEXT) \ - lib_a-log2l.$(OBJEXT) + lib_a-log2l.$(OBJEXT) lib_a-sl_finite.$(OBJEXT) @HAVE_LONG_DOUBLE_TRUE@@USE_LIBTOOL_FALSE@am__objects_4 = \ @HAVE_LONG_DOUBLE_TRUE@@USE_LIBTOOL_FALSE@ $(am__objects_3) @USE_LIBTOOL_FALSE@am_lib_a_OBJECTS = $(am__objects_1) \ @@ -173,7 +173,7 @@ am__objects_7 = atanl.lo cosl.lo sinl.lo tanl.lo tanhl.lo frexpl.lo \ lroundl.lo llroundl.lo truncl.lo remquol.lo fdiml.lo fmaxl.lo \ fminl.lo fmal.lo acoshl.lo atanhl.lo remainderl.lo lgammal.lo \ erfl.lo erfcl.lo logbl.lo nexttowardf.lo nexttoward.lo \ - nexttowardl.lo log2l.lo + nexttowardl.lo log2l.lo sl_finite.lo @HAVE_LONG_DOUBLE_TRUE@@USE_LIBTOOL_TRUE@am__objects_8 = \ @HAVE_LONG_DOUBLE_TRUE@@USE_LIBTOOL_TRUE@ $(am__objects_7) @USE_LIBTOOL_TRUE@am_libcommon_la_OBJECTS = $(am__objects_5) \ @@ -360,7 +360,8 @@ lsrc = atanl.c cosl.c sinl.c tanl.c tanhl.c frexpl.c modfl.c ceill.c fabsl.c \ scalbnl.c exp2l.c scalblnl.c tgammal.c nearbyintl.c lrintl.c llrintl.c \ roundl.c lroundl.c llroundl.c truncl.c remquol.c fdiml.c fmaxl.c fminl.c \ fmal.c acoshl.c atanhl.c remainderl.c lgammal.c erfl.c erfcl.c \ - logbl.c nexttowardf.c nexttoward.c nexttowardl.c log2l.c + logbl.c nexttowardf.c nexttoward.c nexttowardl.c log2l.c \ + sl_finite.c libcommon_la_LDFLAGS = -Xcompiler -nostdlib @USE_LIBTOOL_TRUE@noinst_LTLIBRARIES = libcommon.la @@ -1238,6 +1239,12 @@ lib_a-log2l.o: log2l.c lib_a-log2l.obj: log2l.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-log2l.obj `if test -f 'log2l.c'; then $(CYGPATH_W) 'log2l.c'; else $(CYGPATH_W) '$(srcdir)/log2l.c'; fi` +lib_a-sl_finite.o: sl_finite.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sl_finite.o `test -f 'sl_finite.c' || echo '$(srcdir)/'`sl_finite.c + +lib_a-sl_finite.obj: sl_finite.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sl_finite.obj `if test -f 'sl_finite.c'; then $(CYGPATH_W) 'sl_finite.c'; else $(CYGPATH_W) '$(srcdir)/sl_finite.c'; fi` + mostlyclean-libtool: -rm -f *.lo diff --git a/newlib/libm/common/fdlibm.h b/newlib/libm/common/fdlibm.h index a4b7fffe7..821e4dedb 100644 --- a/newlib/libm/common/fdlibm.h +++ b/newlib/libm/common/fdlibm.h @@ -146,6 +146,8 @@ extern double scalb __P((double, double)); #endif extern double significand __P((double)); +extern long double __ieee754_hypotl __P((long double, long double)); + /* ieee style elementary functions */ extern double __ieee754_sqrt __P((double)); extern double __ieee754_acos __P((double)); diff --git a/newlib/libm/common/hypotl.c b/newlib/libm/common/hypotl.c index 3934b8ecc..cf67ccf58 100644 --- a/newlib/libm/common/hypotl.c +++ b/newlib/libm/common/hypotl.c @@ -29,14 +29,61 @@ POSSIBILITY OF SUCH DAMAGE. */ #include <math.h> -#include "local.h" +#include <errno.h> +#include "fdlibm.h" -/* On platforms where long double is as wide as double. */ -#ifdef _LDBL_EQ_DBL long double hypotl (long double x, long double y) { +#ifdef _LDBL_EQ_DBL + + /* On platforms where long double is as wide as double. */ return hypot(x, y); -} + +#else + + long double z; + + z = __ieee754_hypotl (x, y); + + if (_LIB_VERSION == _IEEE_) + return z; + + if ((! finitel (z)) && finitel (x) && finitel (y)) + { + /* hypot (finite, finite) overflow. */ + struct exception exc; + + exc.type = OVERFLOW; + exc.name = "hypotl"; + exc.err = 0; + exc.arg1 = x; + exc.arg2 = y; + + if (_LIB_VERSION == _SVID_) + exc.retval = HUGE; + else + { +#ifndef HUGE_VAL +#define HUGE_VAL inf + double inf = 0.0; + + SET_HIGH_WORD (inf, 0x7ff00000); /* Set inf to infinite. */ #endif + exc.retval = HUGE_VAL; + } + if (_LIB_VERSION == _POSIX_) + errno = ERANGE; + else if (! matherr (& exc)) + errno = ERANGE; + + if (exc.err != 0) + errno = exc.err; + + return (long double) exc.retval; + } + + return z; +#endif /* ! _LDBL_EQ_DBL */ +} diff --git a/newlib/libm/common/sl_finite.c b/newlib/libm/common/sl_finite.c new file mode 100644 index 000000000..0c27b4fae --- /dev/null +++ b/newlib/libm/common/sl_finite.c @@ -0,0 +1,25 @@ +/* Copyright (C) 2015 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. + */ + +/* finitel(x) returns 1 is x is finite, else 0; */ + +#include <math.h> + +int +finitel (long double x) +{ +#ifdef _LDBL_EQ_DBL + return finite (x); +#else + /* Let the compiler do this for us. + Note - we do not know how many bits there are in a long double. + Some architectures for example have an 80-bit long double whereas + others use 128-bits. We use macros and comiler builtin functions + to avoid specific knowledge of the long double format. */ + return __builtin_isinf_sign (x) == 0; +#endif +} + diff --git a/newlib/libm/common/sqrtl.c b/newlib/libm/common/sqrtl.c index 1e7d7c1d6..8cd717e19 100644 --- a/newlib/libm/common/sqrtl.c +++ b/newlib/libm/common/sqrtl.c @@ -31,12 +31,165 @@ POSSIBILITY OF SUCH DAMAGE. #include <math.h> #include "local.h" -/* On platforms where long double is as wide as double. */ #ifdef _LDBL_EQ_DBL +/* On platforms where long double is as wide as double. */ long double sqrtl (long double x) { return sqrt(x); } + +#else + + /* This code is based upon the version in the BSD math's library. + That code is... + * + * Copyright (c) 2007 Steven G. Kargl + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice unmodified, this list of conditions, and the following + * disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include <float.h> +#include "ieeefp.h" + +#ifndef LDBL_NBIT +#define LDBL_NBIT 0 #endif +#ifndef LDBL_MAX_EXP +#define LDBL_MAX_EXP DBL_MAX_EXP +#endif + +/* Return (x + ulp) for normal positive x. Assumes no overflow. */ + +static inline long double +inc (long double x) +{ + union ieee_ext_u ux = { .extu_ld = x, }; + + if (++ux.extu_ext.ext_fracl == 0) + { + if (++ux.extu_ext.ext_frach == 0) + { + ux.extu_ext.ext_exp++; + ux.extu_ext.ext_frach |= LDBL_NBIT; + } + } + + return ux.extu_ld; +} + +/* Return (x - ulp) for normal positive x. Assumes no underflow. */ + +static inline long double +dec (long double x) +{ + union ieee_ext_u ux = { .extu_ld = x, }; + + if (ux.extu_ext.ext_fracl-- == 0) + { + if (ux.extu_ext.ext_frach-- == LDBL_NBIT) + { + ux.extu_ext.ext_exp--; + ux.extu_ext.ext_frach |= LDBL_NBIT; + } + } + + return ux.extu_ld; +} + +/* This is slow, but simple and portable. */ + +long double +sqrtl (long double x) +{ + union ieee_ext_u ux = { .extu_ld = x, }; + int k, r; + long double lo, xn; + + /* If x = NaN, then sqrt(x) = NaN. */ + /* If x = Inf, then sqrt(x) = Inf. */ + /* If x = -Inf, then sqrt(x) = NaN. */ + if (ux.extu_ext.ext_exp == LDBL_MAX_EXP * 2 - 1) + return (x * x + x); + + /* If x = +-0, then sqrt(x) = +-0. */ + if (x == 0.0L || x == -0.0L) + return x; + + /* If x < 0, then raise invalid and return NaN. */ + if (ux.extu_ext.ext_sign) + return ((x - x) / (x - x)); + + if (ux.extu_ext.ext_exp == 0) + { + /* Adjust subnormal numbers. */ + ux.extu_ld *= 0x1.0p514; + k = -514; + } + else + k = 0; + + /* ux.extu_ld is a normal number, so break it into ux.extu_ld = e*2^n where + ux.extu_ld = (2*e)*2^2k for odd n and ux.extu_ld = (4*e)*2^2k for even n. */ + + if ((ux.extu_ext.ext_exp - EXT_EXP_BIAS) & 1) + { + /* n is even. */ + k += ux.extu_ext.ext_exp - EXT_EXP_BIAS - 1; /* 2k = n - 2. */ + ux.extu_ext.ext_exp = EXT_EXP_BIAS + 1; /* ux.extu_ld in [2,4). */ + } + else + { + k += ux.extu_ext.ext_exp - EXT_EXP_BIAS; /* 2k = n - 1. */ + ux.extu_ext.ext_exp = EXT_EXP_BIAS; /* ux.extu_ld in [1,2). */ + } + + /* Newton's iteration. + Split ux.extu_ld into a high and low part to achieve additional precision. */ + + xn = sqrt ((double) ux.extu_ld); /* 53-bit estimate of sqrtl(x). */ + +#if LDBL_MANT_DIG > 100 + xn = (xn + (ux.extu_ld / xn)) * 0.5; /* 106-bit estimate. */ +#endif + + lo = ux.extu_ld; + ux.extu_ext.ext_fracl = 0; /* Zero out lower bits. */ + lo = (lo - ux.extu_ld) / xn; /* Low bits divided by xn. */ + xn = xn + (ux.extu_ld / xn); /* High portion of estimate. */ + ux.extu_ld = xn + lo; /* Combine everything. */ + ux.extu_ext.ext_exp += (k >> 1) - 1; + + xn = x / ux.extu_ld; /* Chopped quotient (inexact?). */ + + /* For simplicity we round to nearest. */ + xn = inc (xn); /* xn = xn + ulp. */ + + ux.extu_ld = ux.extu_ld + xn; /* Chopped sum. */ + ux.extu_ext.ext_exp--; + + return ux.extu_ld; +} +#endif /* ! _LDBL_EQ_DBL */ + + |