X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Flogl.c;h=de46bd357b4a944dc701a640c9cdd3daf55c0cba;hb=fba7fdb1432af9bf0abece1ec8078ed8ae3af886;hp=c4738ebc6a8cb93e9a4f96ef9adab75f66652cf7;hpb=a03ce664037ae2d0565aedc481fdd8150276b233;p=gnulib.git diff --git a/lib/logl.c b/lib/logl.c index c4738ebc6..de46bd357 100644 --- a/lib/logl.c +++ b/lib/logl.c @@ -1,4 +1,24 @@ -/* logll.c +/* Copyright 2001 by Stephen L. Moshier + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + This program 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . */ + +#include + +/* Specification. */ +#include + +/* logll.c * * Natural logarithm for 128-bit long double precision. * @@ -42,27 +62,6 @@ * */ -#include - -#include "mathl.h" - -/* Copyright 2001 by Stephen L. Moshier - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2, or (at your option) - any later version. - - This program 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 General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program; see the file COPYING. - If not, write to the Free Software Foundation, - 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ - /* log(1+x) = x - .5 x^2 + x^3 l(x) -.0078125 <= x <= +.0078125 peak relative error 1.2e-37 */ @@ -186,15 +185,19 @@ static const long double ln2b = 1.4286068203094172321214581765680755001344E-6L; long double -logl(long double x) +logl (long double x) { long double z, y, w; - long double u, t; - unsigned int m; + long double t; int k, e; /* Check for IEEE special cases. */ + /* log(NaN) = NaN. */ + if (isnanl (x)) + { + return x; + } /* log(0) = -infinity. */ if (x == 0.0L) { @@ -205,14 +208,14 @@ logl(long double x) { return (x - x) / ZERO; } - /* log (infinity or NaN) */ - if (x + x == x || x != x) + /* log (infinity) */ + if (x + x == x) { return x + x; } /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ - x = frexpl(x, &e); + x = frexpl (x, &e); if (x < 0.703125L) { x += x; @@ -236,17 +239,17 @@ logl(long double x) /* Series expansion of log(1+z). */ w = z * z; y = ((((((((((((l15 * z - + l14) * z - + l13) * z - + l12) * z - + l11) * z - + l10) * z - + l9) * z - + l8) * z - + l7) * z - + l6) * z - + l5) * z - + l4) * z + + l14) * z + + l13) * z + + l12) * z + + l11) * z + + l10) * z + + l9) * z + + l8) * z + + l7) * z + + l6) * z + + l5) * z + + l4) * z + l3) * z * w; y -= 0.5 * w; y += e * ln2b; /* Base 2 exponent offset times ln(2). */