diff options
Diffstat (limited to 'lib/mlibc/options/ansi/musl-generic-math/log1pl.c')
| -rw-r--r-- | lib/mlibc/options/ansi/musl-generic-math/log1pl.c | 177 |
1 files changed, 0 insertions, 177 deletions
diff --git a/lib/mlibc/options/ansi/musl-generic-math/log1pl.c b/lib/mlibc/options/ansi/musl-generic-math/log1pl.c deleted file mode 100644 index 141b5f0..0000000 --- a/lib/mlibc/options/ansi/musl-generic-math/log1pl.c +++ /dev/null @@ -1,177 +0,0 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_log1pl.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> - * - * Permission to use, copy, modify, and distribute this software for any - * purpose with or without fee is hereby granted, provided that the above - * copyright notice and this permission notice appear in all copies. - * - * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES - * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF - * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR - * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES - * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN - * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF - * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. - */ -/* - * Relative error logarithm - * Natural logarithm of 1+x, long double precision - * - * - * SYNOPSIS: - * - * long double x, y, log1pl(); - * - * y = log1pl( x ); - * - * - * DESCRIPTION: - * - * Returns the base e (2.718...) logarithm of 1+x. - * - * The argument 1+x is separated into its exponent and fractional - * parts. If the exponent is between -1 and +1, the logarithm - * of the fraction is approximated by - * - * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x). - * - * Otherwise, setting z = 2(x-1)/x+1), - * - * log(x) = z + z^3 P(z)/Q(z). - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20 - */ - -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -long double log1pl(long double x) -{ - return log1p(x); -} -#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x) - * 1/sqrt(2) <= x < sqrt(2) - * Theoretical peak relative error = 2.32e-20 - */ -static const long double P[] = { - 4.5270000862445199635215E-5L, - 4.9854102823193375972212E-1L, - 6.5787325942061044846969E0L, - 2.9911919328553073277375E1L, - 6.0949667980987787057556E1L, - 5.7112963590585538103336E1L, - 2.0039553499201281259648E1L, -}; -static const long double Q[] = { -/* 1.0000000000000000000000E0,*/ - 1.5062909083469192043167E1L, - 8.3047565967967209469434E1L, - 2.2176239823732856465394E2L, - 3.0909872225312059774938E2L, - 2.1642788614495947685003E2L, - 6.0118660497603843919306E1L, -}; - -/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), - * where z = 2(x-1)/(x+1) - * 1/sqrt(2) <= x < sqrt(2) - * Theoretical peak relative error = 6.16e-22 - */ -static const long double R[4] = { - 1.9757429581415468984296E-3L, --7.1990767473014147232598E-1L, - 1.0777257190312272158094E1L, --3.5717684488096787370998E1L, -}; -static const long double S[4] = { -/* 1.00000000000000000000E0L,*/ --2.6201045551331104417768E1L, - 1.9361891836232102174846E2L, --4.2861221385716144629696E2L, -}; -static const long double C1 = 6.9314575195312500000000E-1L; -static const long double C2 = 1.4286068203094172321215E-6L; - -#define SQRTH 0.70710678118654752440L - -long double log1pl(long double xm1) -{ - long double x, y, z; - int e; - - if (isnan(xm1)) - return xm1; - if (xm1 == INFINITY) - return xm1; - if (xm1 == 0.0) - return xm1; - - x = xm1 + 1.0; - - /* Test for domain errors. */ - if (x <= 0.0) { - if (x == 0.0) - return -1/(x*x); /* -inf with divbyzero */ - return 0/0.0f; /* nan with invalid */ - } - - /* Separate mantissa from exponent. - Use frexp so that denormal numbers will be handled properly. */ - x = frexpl(x, &e); - - /* logarithm using log(x) = z + z^3 P(z)/Q(z), - where z = 2(x-1)/x+1) */ - if (e > 2 || e < -2) { - if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ - e -= 1; - z = x - 0.5; - y = 0.5 * z + 0.5; - } else { /* 2 (x-1)/(x+1) */ - z = x - 0.5; - z -= 0.5; - y = 0.5 * x + 0.5; - } - x = z / y; - z = x*x; - z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); - z = z + e * C2; - z = z + x; - z = z + e * C1; - return z; - } - - /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ - if (x < SQRTH) { - e -= 1; - if (e != 0) - x = 2.0 * x - 1.0; - else - x = xm1; - } else { - if (e != 0) - x = x - 1.0; - else - x = xm1; - } - z = x*x; - y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); - y = y + e * C2; - z = y - 0.5 * z; - z = z + x; - z = z + e * C1; - return z; -} -#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 -// TODO: broken implementation to make things compile -long double log1pl(long double x) -{ - return log1p(x); -} -#endif |