diff options
Diffstat (limited to 'newlib/libm/common/sf_log2.c')
| -rw-r--r-- | newlib/libm/common/sf_log2.c | 131 |
1 files changed, 91 insertions, 40 deletions
diff --git a/newlib/libm/common/sf_log2.c b/newlib/libm/common/sf_log2.c index 9fbaaca72..f52cd385d 100644 --- a/newlib/libm/common/sf_log2.c +++ b/newlib/libm/common/sf_log2.c @@ -1,48 +1,99 @@ -/* sf_log2.c -- float version of s_log2.c. - * Modification of sf_exp10.c by Yaakov Selkowitz 2009. - */ +/* Single-precision log2 function. + Copyright (c) 2017 ARM Ltd. All rights reserved. -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * wrapper log2f(x) - */ + 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, 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. + 3. The name of the company may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY ARM LTD ``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 ARM LTD 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 "fdlibm.h" -#include <errno.h> +#if !__OBSOLETE_MATH + #include <math.h> -#undef log2 -#undef log2f - -#ifdef __STDC__ - float log2f(float x) /* wrapper log2f */ -#else - float log2f(x) /* wrapper log2f */ - float x; -#endif -{ - return (logf(x) / (float_t) M_LN2); -} +#include <stdint.h> +#include "math_config.h" -#ifdef _DOUBLE_IS_32BITS +/* +LOG2F_TABLE_BITS = 4 +LOG2F_POLY_ORDER = 4 -#ifdef __STDC__ - double log2(double x) -#else - double log2(x) - double x; -#endif +ULP error: 0.752 (nearest rounding.) +Relative error: 1.9 * 2^-26 (before rounding.) +*/ + +#define N (1 << LOG2F_TABLE_BITS) +#define T __log2f_data.tab +#define A __log2f_data.poly +#define OFF 0x3f330000 + +float +log2f (float x) { - return (double) log2f((float) x); -} + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t z, r, r2, p, y, y0, invc, logc; + uint32_t ix, iz, top, tmp; + int k, i; + + ix = asuint (x); +#if WANT_ROUNDING + /* Fix sign of zero with downward rounding when x==1. */ + if (__builtin_expect (ix == 0x3f800000, 0)) + return 0; +#endif + if (__builtin_expect (ix - 0x00800000 >= 0x7f800000 - 0x00800000, 0)) + { + /* x < 0x1p-126 or inf or nan. */ + if (ix * 2 == 0) + return __math_divzerof (1); + if (ix == 0x7f800000) /* log2(inf) == inf. */ + return x; + if ((ix & 0x80000000) || ix * 2 >= 0xff000000) + return __math_invalidf (x); + /* x is subnormal, normalize it. */ + ix = asuint (x * 0x1p23f); + ix -= 23 << 23; + } + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - LOG2F_TABLE_BITS)) % N; + top = tmp & 0xff800000; + iz = ix - top; + k = (int32_t) tmp >> 23; /* arithmetic shift */ + invc = T[i].invc; + logc = T[i].logc; + z = (double_t) asfloat (iz); -#endif /* defined(_DOUBLE_IS_32BITS) */ + /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ + r = z * invc - 1; + y0 = logc + (double_t) k; + + /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ + r2 = r * r; + y = A[1] * r + A[2]; + y = A[0] * r2 + y; + p = A[3] * r + y0; + y = y * r2 + p; + return (float) y; +} +#endif /* !__OBSOLETE_MATH */ |