Re: cbrt portability

[Michael Goffioul]

On Wed, Feb 29, 2012 at 10:45 PM, Bruno Haible <address@hidden> wrote:
> This patch provides a replacement for platforms without the cbrt() function.
>
>
> 2012-02-29  Bruno Haible  <address@hidden>
>
>        cbrt: Provide replacement on MSVC and Minix.
>        * lib/math.in.h (cbrt): New declaration.
>        * lib/cbrt.c: New file.
>        * m4/cbrt.m4: New file.
>        * m4/math_h.m4 (gl_MATH_H): Test whether cbrt is declared.
>        (gl_MATH_H_DEFAULTS): Initialize GNULIB_CBRT, HAVE_CBRT.
>        * modules/math (Makefile.am): Substitute GNULIB_CBRT, HAVE_CBRT.
>        * modules/cbrt (Files): Add lib/cbrt.c, m4/cbrt.m4.
>        (Depends-on): Add dependencies.
>        (configure.ac): Arrange to compile replacement if HAVE_CBRT is 0.
>        * tests/test-math-c++.cc: Check the declaration of cbrt.
>        * doc/posix-functions/cbrt.texi: Mention that the module provides a
>        replacement.
>
> ================================= lib/cbrt.c =================================
> /* Compute cubic root of double value.
>   Copyright (C) 1997, 2012 Free Software Foundation, Inc.
>
>   Contributed by Dirk Alboth <address@hidden> and
>   Ulrich Drepper <address@hidden>, 1997.
>
>   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>
>
> /* Code based on glibc/sysdeps/ieee754/dbl-64/s_cbrt.c.  */
>
> #define CBRT2 1.2599210498948731648             /* 2^(1/3) */
> #define SQR_CBRT2 1.5874010519681994748         /* 2^(2/3) */
>
> static const double factor[5] =
> {
>  1.0 / SQR_CBRT2,
>  1.0 / CBRT2,
>  1.0,
>  CBRT2,
>  SQR_CBRT2
> };
>
>
> double
> cbrt (double x)
> {
>  if (isfinite (x) && x != 0.0)
>    {
>      double xm, ym, u, t2;
>      int xe;
>
>      /* Reduce X.  XM now is an range 1.0 to 0.5.  */
>      xm = frexp (fabs (x), &xe);
>
>      u = (0.354895765043919860
>           + ((1.50819193781584896
>               + ((-2.11499494167371287
>                   + ((2.44693122563534430
>                       + ((-1.83469277483613086
>                           + (0.784932344976639262 - 0.145263899385486377 * xm)
>                             * xm)
>                          * xm))
>                      * xm))
>                  * xm))
>              * xm));
>
>      t2 = u * u * u;
>
>      ym = u * (t2 + 2.0 * xm) / (2.0 * t2 + xm) * factor[2 + xe % 3];
>
>      return ldexp (x > 0.0 ? ym : -ym, xe / 3);
>    }
>  else
>    return x + x;
> }
> ================================= m4/cbrt.m4 =================================
> # cbrt.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_CBRT],
> [
>  AC_REQUIRE([gl_MATH_H_DEFAULTS])
>
>  dnl Determine CBRT_LIBM.
>  gl_COMMON_DOUBLE_MATHFUNC([cbrt])
>
>  dnl Test whether cbrt() exists.
>  save_LIBS="$LIBS"
>  LIBS="$LIBS $CBRT_LIBM"
>  AC_CHECK_FUNCS([cbrt])
>  LIBS="$save_LIBS"
>  if test $ac_cv_func_cbrt = no; then
>    HAVE_CBRT=0
>    dnl Find libraries needed to link lib/cbrt.c.
>    AC_REQUIRE([gl_FUNC_FABS])
>    AC_REQUIRE([gl_FUNC_FREXP])
>    AC_REQUIRE([gl_FUNC_LDEXP])
>    CBRT_LIBM=
>    dnl Append $FABS_LIBM to CBRT_LIBM, avoiding gratuitous duplicates.
>    case " $CBRT_LIBM " in
>      *" $FABS_LIBM "*) ;;
>      *) CBRT_LIBM="$CBRT_LIBM $FABS_LIBM" ;;
>    esac
>    dnl Append $FREXP_LIBM to CBRT_LIBM, avoiding gratuitous duplicates.
>    case " $CBRT_LIBM " in
>      *" $FREXP_LIBM "*) ;;
>      *) CBRT_LIBM="$CBRT_LIBM $FREXP_LIBM" ;;
>    esac
>    dnl Append $LDEXP_LIBM to CBRT_LIBM, avoiding gratuitous duplicates.
>    case " $CBRT_LIBM " in
>      *" $LDEXP_LIBM "*) ;;
>      *) CBRT_LIBM="$CBRT_LIBM $LDEXP_LIBM" ;;
>    esac
>  fi
> ])
> ==============================================================================
> --- doc/posix-functions/cbrt.texi.orig  Wed Feb 29 23:42:05 2012
> +++ doc/posix-functions/cbrt.texi       Wed Feb 29 23:08:54 2012
> @@ -8,11 +8,11 @@
>
>  Portability problems fixed by Gnulib:
> address@hidden
> address@hidden
> +This function is missing on some platforms:
> +Minix 3.1.8, MSVC 9.
> address@hidden itemize
>
>  Portability problems not fixed by Gnulib:
> address@hidden
> address@hidden
> -This function is missing on some platforms:
> -Minix 3.1.8, MSVC 9.
> address@hidden itemize
> --- lib/math.in.h.orig  Wed Feb 29 23:42:05 2012
> +++ lib/math.in.h       Wed Feb 29 23:07:01 2012
> @@ -288,6 +288,21 @@
>  #endif
>
>
> +#if @GNULIB_CBRT@
> +# if address@hidden@
> +_GL_FUNCDECL_SYS (cbrt, double, (double x));
> +# endif
> +_GL_CXXALIAS_SYS (cbrt, double, (double x));
> +_GL_CXXALIASWARN (cbrt);
> +#elif defined GNULIB_POSIXCHECK
> +# undef cbrt
> +# if HAVE_RAW_DECL_CBRT
> +_GL_WARN_ON_USE (cbrt, "cbrt is unportable - "
> +                 "use gnulib module cbrt for portability");
> +# endif
> +#endif
> +
> +
>  #if @GNULIB_CEILF@
>  # if @REPLACE_CEILF@
>  #  if !(defined __cplusplus && defined GNULIB_NAMESPACE)
> --- m4/math_h.m4.orig   Wed Feb 29 23:42:05 2012
> +++ m4/math_h.m4        Wed Feb 29 23:08:38 2012
> @@ -1,4 +1,4 @@
> -# math_h.m4 serial 75
> +# math_h.m4 serial 76
>  dnl Copyright (C) 2007-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,
> @@ -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
> -     ceilf ceill copysign copysignf copysignl cosf cosl coshf
> +     cbrt 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
> @@ -67,6 +67,7 @@
>   GNULIB_ATANF=0;      AC_SUBST([GNULIB_ATANF])
>   GNULIB_ATANL=0;      AC_SUBST([GNULIB_ATANL])
>   GNULIB_ATAN2F=0;     AC_SUBST([GNULIB_ATAN2F])
> +  GNULIB_CBRT=0;       AC_SUBST([GNULIB_CBRT])
>   GNULIB_CEIL=0;       AC_SUBST([GNULIB_CEIL])
>   GNULIB_CEILF=0;      AC_SUBST([GNULIB_CEILF])
>   GNULIB_CEILL=0;      AC_SUBST([GNULIB_CEILL])
> @@ -141,6 +142,7 @@
>   HAVE_ATANF=1;                AC_SUBST([HAVE_ATANF])
>   HAVE_ATANL=1;                AC_SUBST([HAVE_ATANL])
>   HAVE_ATAN2F=1;               AC_SUBST([HAVE_ATAN2F])
> +  HAVE_CBRT=1;                 AC_SUBST([HAVE_CBRT])
>   HAVE_COPYSIGN=1;             AC_SUBST([HAVE_COPYSIGN])
>   HAVE_COPYSIGNF=1;            AC_SUBST([HAVE_COPYSIGNF])
>   HAVE_COPYSIGNL=1;            AC_SUBST([HAVE_COPYSIGNL])
> --- modules/cbrt.orig   Wed Feb 29 23:42:05 2012
> +++ modules/cbrt        Wed Feb 29 23:40:20 2012
> @@ -2,12 +2,23 @@
>  cbrt() function: cube root.
>
>  Files:
> +lib/cbrt.c
> +m4/cbrt.m4
>  m4/mathfunc.m4
>
>  Depends-on:
> +math
> +isfinite        [test $HAVE_CBRT = 0]
> +fabs            [test $HAVE_CBRT = 0]
> +frexp           [test $HAVE_CBRT = 0]
> +ldexp           [test $HAVE_CBRT = 0]
>
>  configure.ac:
> -gl_COMMON_DOUBLE_MATHFUNC([cbrt])
> +gl_FUNC_CBRT
> +if test $HAVE_CBRT = 0; then
> +  AC_LIBOBJ([cbrt])
> +fi
> +gl_MATH_MODULE_INDICATOR([cbrt])
>
>  Makefile.am:
>
> --- modules/math.orig   Wed Feb 29 23:42:05 2012
> +++ modules/math        Wed Feb 29 23:07:43 2012
> @@ -35,6 +35,7 @@
>              -e 's/@''GNULIB_ATANF''@/$(GNULIB_ATANF)/g' \
>              -e 's/@''GNULIB_ATANL''@/$(GNULIB_ATANL)/g' \
>              -e 's/@''GNULIB_ATAN2F''@/$(GNULIB_ATAN2F)/g' \
> +             -e 's/@''GNULIB_CBRT''@/$(GNULIB_CBRT)/g' \
>              -e 's/@''GNULIB_CEIL''@/$(GNULIB_CEIL)/g' \
>              -e 's/@''GNULIB_CEILF''@/$(GNULIB_CEILF)/g' \
>              -e 's/@''GNULIB_CEILL''@/$(GNULIB_CEILL)/g' \
> @@ -109,6 +110,7 @@
>              -e 's|@''HAVE_ATANF''@|$(HAVE_ATANF)|g' \
>              -e 's|@''HAVE_ATANL''@|$(HAVE_ATANL)|g' \
>              -e 's|@''HAVE_ATAN2F''@|$(HAVE_ATAN2F)|g' \
> +             -e 's|@''HAVE_CBRT''@|$(HAVE_CBRT)|g' \
>              -e 's|@''HAVE_COPYSIGN''@|$(HAVE_COPYSIGN)|g' \
>              -e 's|@''HAVE_COPYSIGNF''@|$(HAVE_COPYSIGNF)|g' \
>              -e 's|@''HAVE_COPYSIGNL''@|$(HAVE_COPYSIGNL)|g' \
> --- tests/test-math-c++.cc.orig Wed Feb 29 23:42:05 2012
> +++ tests/test-math-c++.cc      Wed Feb 29 23:11:46 2012
> @@ -68,7 +68,9 @@
>  #endif
>  //SIGNATURE_CHECK (GNULIB_NAMESPACE::atan2, double, (double, double));
>
> -//SIGNATURE_CHECK (GNULIB_NAMESPACE::cbrt, double, (double));
> +#if GNULIB_TEST_CBRT
> +SIGNATURE_CHECK (GNULIB_NAMESPACE::cbrt, double, (double));
> +#endif
>
>  #if GNULIB_TEST_CEILF
>  SIGNATURE_CHECK (GNULIB_NAMESPACE::ceilf, float, (float));
>
>

John,

Could we make octave use it instead of having its own implementation?
This is known to fail some tests (I reported the issue some time ago).

Michael.