X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Flogl.c;h=0c00a2b22857aaedc34d8f5a196a859461b9452b;hb=7ef6c64e210ac0979d7e8ac69bc5b5208c2405ab;hp=374b731fddf962fc66374ea54af65441443b23e1;hpb=267a39bafd249d7eb9c37df06dc6defcf41cb343;p=gnulib.git

diff --git a/lib/logl.c b/lib/logl.c
index 374b731fd..0c00a2b22 100644
--- a/lib/logl.c
+++ b/lib/logl.c
@@ -1,4 +1,49 @@
-/*							logll.c
+/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
+
+   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 <http://www.gnu.org/licenses/>.  */
+
+#include <config.h>
+
+/* Specification.  */
+#include <math.h>
+
+#if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
+
+long double
+logl (long double x)
+{
+  return log (x);
+}
+
+#elif HAVE_LOGL
+
+long double
+logl (long double x)
+# undef logl
+{
+  /* Work around the OSF/1 5.1 bug.  */
+  if (x == 0.0L)
+    /* Return -Infinity.  */
+    return -1.0L / 0.0L;
+  return logl (x);
+}
+
+#else
+
+/* Code based on glibc/sysdeps/ieee754/ldbl-128/e_logl.c.  */
+
+/*                                                      logll.c
  *
  * Natural logarithm for 128-bit long double precision.
  *
@@ -42,27 +87,6 @@
  *
  */
 
-#include <math.h>
-
-#include "mathl.h"
-
-/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
-
-   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,
-   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
-
 /* log(1+x) = x - .5 x^2 + x^3 l(x)
    -.0078125 <= x <= +.0078125
    peak relative error 1.2e-37 */
@@ -123,7 +147,7 @@ static const long double logtbl[92] = {
 -2.7902661731604211834685052867305795169688E-4L,
 -1.2335696813916860754951146082826952093496E-4L,
 -3.0677461025892873184042490943581654591817E-5L,
-#define ZERO logtbl[38]
+# define ZERO logtbl[38]
  0.0000000000000000000000000000000000000000E0L,
 -3.0359557945051052537099938863236321874198E-5L,
 -1.2081346403474584914595395755316412213151E-4L,
@@ -186,15 +210,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 +233,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 +264,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).  */
@@ -256,3 +284,5 @@ logl(long double x)
   y += e * ln2a;
   return y;
 }
+
+#endif