new module 'cbrtl'

[Bruno Haible]

Here is a new module that provides the cbrtl() function.


2012-02-29  Bruno Haible  <address@hidden>

        New module 'cbrtl'.
        * lib/math.in.h (cbrtl): New declaration.
        * lib/cbrtl.c: New file.
        * m4/cbrtl.m4: New file.
        * m4/math_h.m4 (gl_MATH_H): Test whether cbrtl is declared.
        (gl_MATH_H_DEFAULTS): Initialize GNULIB_CBRTL, HAVE_CBRTL,
        HAVE_DECL_CBRTL.
        * modules/math (Makefile.am): Substitute GNULIB_CBRTL, HAVE_CBRTL,
        HAVE_DECL_CBRTL.
        * modules/cbrtl: New file.
        * tests/test-math-c++.cc: Check the declaration of cbrtl.
        * doc/posix-functions/cbrtl.texi: Mention the new module.

================================= lib/cbrtl.c =================================
/* Compute cubic root of long double value.
   Copyright (C) 2012 Free Software Foundation, Inc.
   Cephes Math Library Release 2.2: January, 1991
   Copyright 1984, 1991 by Stephen L. Moshier
   Adapted for glibc October, 2001.

   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
cbrtl (long double x)
{
  return cbrt (x);
}

#else

/* Code based on glibc/sysdeps/ieee754/ldbl-128/s_cbrtl.c.  */

/*                                                      cbrtl.c
 *
 *      Cube root, long double precision
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, cbrtl();
 *
 * y = cbrtl( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the cube root of the argument, which may be negative.
 *
 * Range reduction involves determining the power of 2 of
 * the argument.  A polynomial of degree 2 applied to the
 * mantissa, and multiplication by the cube root of 1, 2, or 4
 * approximates the root to within about 0.1%.  Then Newton's
 * iteration is used three times to converge to an accurate
 * result.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE       -8,8       100000      1.3e-34     3.9e-35
 *    IEEE    exp(+-707)    100000      1.3e-34     4.3e-35
 *
 */

static const long double CBRT2 = 1.259921049894873164767210607278228350570251L;
static const long double CBRT4 = 1.587401051968199474751705639272308260391493L;
static const long double CBRT2I = 
0.7937005259840997373758528196361541301957467L;
static const long double CBRT4I = 
0.6299605249474365823836053036391141752851257L;

long double
cbrtl (long double x)
{
  if (isfinite (x) && x != 0.0L)
    {
      int e, rem, sign;
      long double z;

      if (x > 0)
        sign = 1;
      else
        {
          sign = -1;
          x = -x;
        }

      z = x;
      /* extract power of 2, leaving mantissa between 0.5 and 1  */
      x = frexpl (x, &e);

      /* Approximate cube root of number between .5 and 1,
         peak relative error = 1.2e-6  */
      x = ((((1.3584464340920900529734e-1L * x
              - 6.3986917220457538402318e-1L) * x
             + 1.2875551670318751538055e0L) * x
            - 1.4897083391357284957891e0L) * x
           + 1.3304961236013647092521e0L) * x + 3.7568280825958912391243e-1L;

      /* exponent divided by 3 */
      if (e >= 0)
        {
          rem = e;
          e /= 3;
          rem -= 3 * e;
          if (rem == 1)
            x *= CBRT2;
          else if (rem == 2)
            x *= CBRT4;
        }
      else
        {                           /* argument less than 1 */
          e = -e;
          rem = e;
          e /= 3;
          rem -= 3 * e;
          if (rem == 1)
            x *= CBRT2I;
          else if (rem == 2)
            x *= CBRT4I;
          e = -e;
        }

      /* multiply by power of 2 */
      x = ldexpl (x, e);

      /* Newton iteration */
      x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
      x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
      x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;

      if (sign < 0)
        x = -x;
      return x;
    }
  else
    return x + x;
}

#endif
================================= m4/cbrtl.m4 =================================
# cbrtl.m4 serial 1
dnl Copyright (C) 2012 Free Software Foundation, Inc.
dnl This file is free software; the Free Software Foundation
dnl gives unlimited permission to copy and/or distribute it,
dnl with or without modifications, as long as this notice is preserved.

AC_DEFUN([gl_FUNC_CBRTL],
[
  AC_REQUIRE([gl_MATH_H_DEFAULTS])
  AC_REQUIRE([gl_LONG_DOUBLE_VS_DOUBLE])
  AC_REQUIRE([gl_FUNC_CBRT])

  dnl Test whether cbrtl() exists. Assume that cbrtl(), if it exists, is
  dnl defined in the same library as cbrt().
  save_LIBS="$LIBS"
  LIBS="$LIBS $CBRT_LIBM"
  AC_CHECK_FUNCS([cbrtl])
  LIBS="$save_LIBS"
  if test $ac_cv_func_cbrtl = yes; then
    CBRTL_LIBM="$CBRT_LIBM"
    dnl Also check whether it's declared.
    dnl IRIX 6.5 has cbrtl() in libm but doesn't declare it in <math.h>.
    AC_CHECK_DECL([cbrtl], , [HAVE_DECL_CBRTL=0], [[#include <math.h>]])
  else
    HAVE_CBRTL=0
    HAVE_DECL_CBRTL=0
    dnl Find libraries needed to link lib/cbrtl.c.
    if test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 1; then
      CBRTL_LIBM="$CBRT_LIBM"
    else
      AC_REQUIRE([gl_FUNC_FREXPL])
      AC_REQUIRE([gl_FUNC_LDEXPL])
      CBRTL_LIBM=
      dnl Append $FREXPL_LIBM to CBRTL_LIBM, avoiding gratuitous duplicates.
      case " $CBRTL_LIBM " in
        *" $FREXPL_LIBM "*) ;;
        *) CBRTL_LIBM="$CBRTL_LIBM $FREXPL_LIBM" ;;
      esac
      dnl Append $LDEXPL_LIBM to CBRTL_LIBM, avoiding gratuitous duplicates.
      case " $CBRTL_LIBM " in
        *" $LDEXPL_LIBM "*) ;;
        *) CBRTL_LIBM="$CBRTL_LIBM $LDEXPL_LIBM" ;;
      esac
    fi
  fi
  AC_SUBST([CBRTL_LIBM])
])
================================ modules/cbrtl ================================
Description:
cbrtl() function: cube root.

Files:
lib/cbrtl.c
m4/cbrtl.m4
m4/cbrt.m4
m4/mathfunc.m4

Depends-on:
math
cbrt            [test $HAVE_CBRTL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 
= 1]
isfinite        [test $HAVE_CBRTL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 
= 0]
frexpl          [test $HAVE_CBRTL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 
= 0]
ldexpl          [test $HAVE_CBRTL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 
= 0]

configure.ac:
gl_FUNC_CBRTL
if test $HAVE_CBRTL = 0; then
  AC_LIBOBJ([cbrtl])
fi
gl_MATH_MODULE_INDICATOR([cbrtl])

Makefile.am:

Include:
<math.h>

Link:
$(CBRTL_LIBM)

License:
LGPL

Maintainer:
Bruno Haible
===============================================================================
--- doc/posix-functions/cbrtl.texi.orig Thu Mar  1 02:46:16 2012
+++ doc/posix-functions/cbrtl.texi      Thu Mar  1 02:32:37 2012
@@ -4,15 +4,18 @@
 
 POSIX specification:@* 
@url{http://www.opengroup.org/onlinepubs/9699919799/functions/cbrtl.html}
 
-Gnulib module: ---
+Gnulib module: cbrtl
 
 Portability problems fixed by Gnulib:
 @itemize
address@hidden
+This function is missing on some platforms:
+FreeBSD 6.0, NetBSD 5.0, OpenBSD 3.8, Minix 3.1.8, AIX 5.1, HP-UX 11, Solaris 
9, Cygwin, MSVC 9, Interix 3.5, BeOS.
address@hidden
+This function is not declared on some platforms:
+IRIX 6.5.
 @end itemize
 
 Portability problems not fixed by Gnulib:
 @itemize
address@hidden
-This function is missing on some platforms:
-FreeBSD 6.0, NetBSD 5.0, OpenBSD 3.8, Minix 3.1.8, AIX 5.1, HP-UX 11, IRIX 
6.5, Solaris 9, Cygwin, MSVC 9, Interix 3.5, BeOS.
 @end itemize
--- lib/math.in.h.orig  Thu Mar  1 02:46:16 2012
+++ lib/math.in.h       Thu Mar  1 02:31:26 2012
@@ -316,6 +316,20 @@
 # endif
 #endif
 
+#if @GNULIB_CBRTL@
+# if address@hidden@
+_GL_FUNCDECL_SYS (cbrtl, long double, (long double x));
+# endif
+_GL_CXXALIAS_SYS (cbrtl, long double, (long double x));
+_GL_CXXALIASWARN (cbrtl);
+#elif defined GNULIB_POSIXCHECK
+# undef cbrtl
+# if HAVE_RAW_DECL_CBRTL
+_GL_WARN_ON_USE (cbrtl, "cbrtl is unportable - "
+                 "use gnulib module cbrtl for portability");
+# endif
+#endif
+
 
 #if @GNULIB_CEILF@
 # if @REPLACE_CEILF@
--- m4/math_h.m4.orig   Thu Mar  1 02:46:16 2012
+++ m4/math_h.m4        Thu Mar  1 02:31:39 2012
@@ -40,7 +40,7 @@
   dnl corresponding gnulib module is not in use.
   gl_WARN_ON_USE_PREPARE([[#include <math.h>]],
     [acosf acosl asinf asinl atanf atanl
-     cbrt cbrtf ceilf ceill copysign copysignf copysignl cosf cosl coshf
+     cbrt cbrtf cbrtl ceilf ceill copysign copysignf copysignl cosf cosl coshf
      expf expl fabsf fabsl floorf floorl fma fmaf fmal
      fmod fmodf fmodl frexpf frexpl hypotf hypotl
      ldexpf ldexpl logb logf logl log10f log10l modf modff modfl powf
@@ -69,6 +69,7 @@
   GNULIB_ATAN2F=0;     AC_SUBST([GNULIB_ATAN2F])
   GNULIB_CBRT=0;       AC_SUBST([GNULIB_CBRT])
   GNULIB_CBRTF=0;      AC_SUBST([GNULIB_CBRTF])
+  GNULIB_CBRTL=0;      AC_SUBST([GNULIB_CBRTL])
   GNULIB_CEIL=0;       AC_SUBST([GNULIB_CEIL])
   GNULIB_CEILF=0;      AC_SUBST([GNULIB_CEILF])
   GNULIB_CEILL=0;      AC_SUBST([GNULIB_CEILL])
@@ -145,6 +146,7 @@
   HAVE_ATAN2F=1;               AC_SUBST([HAVE_ATAN2F])
   HAVE_CBRT=1;                 AC_SUBST([HAVE_CBRT])
   HAVE_CBRTF=1;                AC_SUBST([HAVE_CBRTF])
+  HAVE_CBRTL=1;                AC_SUBST([HAVE_CBRTL])
   HAVE_COPYSIGN=1;             AC_SUBST([HAVE_COPYSIGN])
   HAVE_COPYSIGNF=1;            AC_SUBST([HAVE_COPYSIGNF])
   HAVE_COPYSIGNL=1;            AC_SUBST([HAVE_COPYSIGNL])
@@ -192,6 +194,7 @@
   HAVE_DECL_ASINL=1;           AC_SUBST([HAVE_DECL_ASINL])
   HAVE_DECL_ATANL=1;           AC_SUBST([HAVE_DECL_ATANL])
   HAVE_DECL_CBRTF=1;           AC_SUBST([HAVE_DECL_CBRTF])
+  HAVE_DECL_CBRTL=1;           AC_SUBST([HAVE_DECL_CBRTL])
   HAVE_DECL_CEILF=1;           AC_SUBST([HAVE_DECL_CEILF])
   HAVE_DECL_CEILL=1;           AC_SUBST([HAVE_DECL_CEILL])
   HAVE_DECL_COSL=1;            AC_SUBST([HAVE_DECL_COSL])
--- modules/math.orig   Thu Mar  1 02:46:16 2012
+++ modules/math        Thu Mar  1 02:31:58 2012
@@ -37,6 +37,7 @@
              -e 's/@''GNULIB_ATAN2F''@/$(GNULIB_ATAN2F)/g' \
              -e 's/@''GNULIB_CBRT''@/$(GNULIB_CBRT)/g' \
              -e 's/@''GNULIB_CBRTF''@/$(GNULIB_CBRTF)/g' \
+             -e 's/@''GNULIB_CBRTL''@/$(GNULIB_CBRTL)/g' \
              -e 's/@''GNULIB_CEIL''@/$(GNULIB_CEIL)/g' \
              -e 's/@''GNULIB_CEILF''@/$(GNULIB_CEILF)/g' \
              -e 's/@''GNULIB_CEILL''@/$(GNULIB_CEILL)/g' \
@@ -113,6 +114,7 @@
              -e 's|@''HAVE_ATAN2F''@|$(HAVE_ATAN2F)|g' \
              -e 's|@''HAVE_CBRT''@|$(HAVE_CBRT)|g' \
              -e 's|@''HAVE_CBRTF''@|$(HAVE_CBRTF)|g' \
+             -e 's|@''HAVE_CBRTL''@|$(HAVE_CBRTL)|g' \
              -e 's|@''HAVE_COPYSIGN''@|$(HAVE_COPYSIGN)|g' \
              -e 's|@''HAVE_COPYSIGNF''@|$(HAVE_COPYSIGNF)|g' \
              -e 's|@''HAVE_COPYSIGNL''@|$(HAVE_COPYSIGNL)|g' \
@@ -160,6 +162,7 @@
              -e 's|@''HAVE_DECL_ASINL''@|$(HAVE_DECL_ASINL)|g' \
              -e 's|@''HAVE_DECL_ATANL''@|$(HAVE_DECL_ATANL)|g' \
              -e 's|@''HAVE_DECL_CBRTF''@|$(HAVE_DECL_CBRTF)|g' \
+             -e 's|@''HAVE_DECL_CBRTL''@|$(HAVE_DECL_CBRTL)|g' \
              -e 's|@''HAVE_DECL_CEILF''@|$(HAVE_DECL_CEILF)|g' \
              -e 's|@''HAVE_DECL_CEILL''@|$(HAVE_DECL_CEILL)|g' \
              -e 's|@''HAVE_DECL_COSL''@|$(HAVE_DECL_COSL)|g' \
--- tests/test-math-c++.cc.orig Thu Mar  1 02:46:16 2012
+++ tests/test-math-c++.cc      Thu Mar  1 02:13:12 2012
@@ -74,6 +74,9 @@
 #if GNULIB_TEST_CBRT
 SIGNATURE_CHECK (GNULIB_NAMESPACE::cbrt, double, (double));
 #endif
+#if GNULIB_TEST_CBRTL
+SIGNATURE_CHECK (GNULIB_NAMESPACE::cbrtl, long double, (long double));
+#endif
 
 #if GNULIB_TEST_CEILF
 SIGNATURE_CHECK (GNULIB_NAMESPACE::ceilf, float, (float));