From b9e7cd9a846cf584049a719bbd0e4cabe3c50bea Mon Sep 17 00:00:00 2001 From: Nick Clifton Date: Fri, 6 Feb 2015 16:14:04 +0000 Subject: * libc/include/complex.h (cabsl): Add prototype. (cimagl): Add prototype. (creall): Add prototype. * libc/include/ieeefp.h: Include float.h. (EXT_EXPBITS, EXT_FRACHBITS, EXT_FRACLBITS) (EXT_EXP_INFNAN. EXT_EXP_BIAS, EXT_FRACBITS): Define. (struct ieee_ext, union ieee_ext_u): New types for long double support. * libc/include/math.h (finitel): Add prototype. (hypotl): Add prototype. (sqrtl): Add prototype. * libm/common/Makefile.am (lsrc): Add sl_finite.c. * libm/common/Makefile.in: Regenerate. * libm/common/fdlibm.h (__ieee754_hypotl): Add prototype. * libm/common/hypotl.c (hypotl): Add implementation for when long double is larger than double. * libm/common/sqrtl.c (sqrtl): Likewise. * libm/common/sl_finite.c: New file. Adds implementation of the finitel function. * libm/complex/Makefile.am (lsrc): Define. (libcomplex_la_SOURCES): Add lsrc. (lib_a_SOURCES): Add lsrc. * libm/complex/Makefile.in: Regenerate. * libm/complex/cabs.c: Add documentation of cabsl function. * libm/complex/cimag.c: Add documentation of cimagl function. * libm/complex/creall.c: Add documentation of creall function. * libm/complex/cabsl.c: New file. Adds implementation of the cabsl function. * libm/complex/cimagl.c: New file. Adds implementation of the cimagl function. * libm/complex/creall.c: New file. Adds implementation of the creall function. * libm/math/Makefile.am (lsrc): Define. (libmath_la_SOURCES): Add lsrc. (lib_a_SOURCES): Add lsrc. * libm/math/Makefile.in: Regenerate. * libm/math/el_hypot.c: New file. Adds implementation of the __ieee754_hypotl function. --- newlib/libm/common/sqrtl.c | 155 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 154 insertions(+), 1 deletion(-) (limited to 'newlib/libm/common/sqrtl.c') 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 #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 +#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 */ + + -- cgit v1.2.3