[libgfortran] PR 16137: libgfortran on non-C99 targets

[libgfortran] PR 16137: libgfortran on non-C99 targets

[Roger Sayle]

The following patch is my proposed solution to PR libgfortran/16137,
which is the inability to link libgfortran on targets without full
C99 support due to missing symbols.  The solution is to provide
replacement/stub implementations (much like libstdc++-v3 does);
checking for the required functions in configure and then providing
implementations in instrinsics/c99_functions.c when not available.

Most stub implementations are "obvious", but the one bit of cleverness
was the implementation of a portable nextafterf implementation thats
independent of the precise floating point format or its layout in
memory.  Hopefully, gfortran's use of nextafterf doesn't depend upon
its operation with denormalized numbers.  My attempts to support denorms
ran afoul of the target libm's inability to handle denormal arguments
to ldexp.  The implementation below always returns normalized values,
for example, returing FLT_MIN as the next value after zero.  I believe
this is conforming, or atleast good enough for fortran work.

The net result is that on sparc-sun-solaris2.8 there are 450 less
unexpected failures, and 978 more expected passes on the gfortran
testsuite.  On first inspection many of the remaining failures look
like LD_LIBRARY_PATH problems failing to find libgmp.so correctly.


There may be issues with people trying to link both C++ and Fortran95
code on Solaris due to duplicate symbols in their runtimes.  However,
I'd expect this to be rare, and no reason to prevent using Fortran on
its own.


The following patch has been tested on both i686-pc-linux-gnu and
sparc-sun-solaris2.8 with a full bootstrap, including gfortran, and
regression tested with a top-level "make -k check" with no new failures.

Ok for mainline?



2004-08-02  Roger Sayle  <roger@eyesopen.com>

	PR libgfortran/16137
	* configure.ac: Add tests for acosf, asinf, atan2f, atanf, ceilf,
	copysignf, cosf, coshf, expf, floorf, frexpf, hypotf, logf, log10f,
	scalbnf, sinf, sinhf, sqrtf, tanf and tanhf in libm.
	* config.h.in: Regenerate.
	* configure: Regenerate.

	* instrinsics/c99_functions.c (acosf, asinf, atan2f, atanf, ceilf,
	copysignf, cosf, coshf, expf, floorf, frexpf, hypotf, logf, log10f,
	nextafterf, scalbnf, sinf, sinhf, sqrtf, tanf, tanhf): New stub
	implementations for targets that don't support C99 float functions.


Index: configure.ac
===================================================================
RCS file: /cvs/gcc/gcc/libgfortran/configure.ac,v
retrieving revision 1.6
diff -c -3 -p -r1.6 configure.ac
*** configure.ac	18 Jul 2004 12:46:39 -0000	1.6
--- configure.ac	2 Aug 2004 02:10:19 -0000
*************** AC_CHECK_LIB([m],[csin],[need_math="no"]
*** 160,171 ****
  # Check for library functions.
  AC_CHECK_FUNCS(getrusage times)

! # Check for some C99 functions
! AC_CHECK_LIB([m],[round],[AC_DEFINE([HAVE_ROUND],[1],["c99 function"])])
! AC_CHECK_LIB([m],[roundf],[AC_DEFINE([HAVE_ROUNDF],[1],["c99 function"])])
! # And other IEEE math functions
  AC_CHECK_LIB([m],[nextafter],[AC_DEFINE([HAVE_NEXTAFTER],[1],[libm includes nextafter])])
  AC_CHECK_LIB([m],[nextafterf],[AC_DEFINE([HAVE_NEXTAFTERF],[1],[libm includes nextafterf])])

  # Let the user override this
  AC_ARG_ENABLE(cmath,
--- 160,190 ----
  # Check for library functions.
  AC_CHECK_FUNCS(getrusage times)

! # Check for C99 (and other IEEE) math functions
! AC_CHECK_LIB([m],[acosf],[AC_DEFINE([HAVE_ACOSF],[1],[libm includes acosf])])
! AC_CHECK_LIB([m],[asinf],[AC_DEFINE([HAVE_ASINF],[1],[libm includes asinf])])
! AC_CHECK_LIB([m],[atan2f],[AC_DEFINE([HAVE_ATAN2F],[1],[libm includes atan2f])])
! AC_CHECK_LIB([m],[atanf],[AC_DEFINE([HAVE_ATANF],[1],[libm includes atanf])])
! AC_CHECK_LIB([m],[ceilf],[AC_DEFINE([HAVE_CEILF],[1],[libm includes ceilf])])
! AC_CHECK_LIB([m],[copysignf],[AC_DEFINE([HAVE_COPYSIGNF],[1],[libm includes copysignf])])
! AC_CHECK_LIB([m],[cosf],[AC_DEFINE([HAVE_COSF],[1],[libm includes cosf])])
! AC_CHECK_LIB([m],[coshf],[AC_DEFINE([HAVE_COSHF],[1],[libm includes coshf])])
! AC_CHECK_LIB([m],[expf],[AC_DEFINE([HAVE_EXPF],[1],[libm includes expf])])
! AC_CHECK_LIB([m],[floorf],[AC_DEFINE([HAVE_FLOORF],[1],[libm includes floorf])])
! AC_CHECK_LIB([m],[frexpf],[AC_DEFINE([HAVE_FREXPF],[1],[libm includes frexpf])])
! AC_CHECK_LIB([m],[hypotf],[AC_DEFINE([HAVE_HYPOTF],[1],[libm includes hypotf])])
! AC_CHECK_LIB([m],[logf],[AC_DEFINE([HAVE_LOGF],[1],[libm includes logf])])
! AC_CHECK_LIB([m],[log10f],[AC_DEFINE([HAVE_LOG10F],[1],[libm includes log10f])])
  AC_CHECK_LIB([m],[nextafter],[AC_DEFINE([HAVE_NEXTAFTER],[1],[libm includes nextafter])])
  AC_CHECK_LIB([m],[nextafterf],[AC_DEFINE([HAVE_NEXTAFTERF],[1],[libm includes nextafterf])])
+ AC_CHECK_LIB([m],[round],[AC_DEFINE([HAVE_ROUND],[1],[libm includes round])])
+ AC_CHECK_LIB([m],[roundf],[AC_DEFINE([HAVE_ROUNDF],[1],[libm includes roundf])])
+ AC_CHECK_LIB([m],[scalbnf],[AC_DEFINE([HAVE_SCALBNF],[1],[libm includes scalbnf])])
+ AC_CHECK_LIB([m],[sinf],[AC_DEFINE([HAVE_SINF],[1],[libm includes sinf])])
+ AC_CHECK_LIB([m],[sinhf],[AC_DEFINE([HAVE_SINHF],[1],[libm includes sinhf])])
+ AC_CHECK_LIB([m],[sqrtf],[AC_DEFINE([HAVE_SQRTF],[1],[libm includes sqrtf])])
+ AC_CHECK_LIB([m],[tanf],[AC_DEFINE([HAVE_TANF],[1],[libm includes tanf])])
+ AC_CHECK_LIB([m],[tanhf],[AC_DEFINE([HAVE_TANHF],[1],[libm includes tanhf])])

  # Let the user override this
  AC_ARG_ENABLE(cmath,
Index: intrinsics/c99_functions.c
===================================================================
RCS file: /cvs/gcc/gcc/libgfortran/intrinsics/c99_functions.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 c99_functions.c
*** intrinsics/c99_functions.c	14 Jun 2004 17:27:20 -0000	1.2
--- intrinsics/c99_functions.c	3 Aug 2004 01:44:47 -0000
*************** Boston, MA 02111-1307, USA.  */
*** 20,29 ****
--- 20,251 ----

  #include "config.h"
  #include <sys/types.h>
+ #include <float.h>
  #include <math.h>
  #include "libgfortran.h"


+ #ifndef HAVE_ACOSF
+ float
+ acosf(float x)
+ {
+   return (float) acos(x);
+ }
+ #endif
+
+ #ifndef HAVE_ASINF
+ float
+ asinf(float x)
+ {
+   return (float) asin(x);
+ }
+ #endif
+
+ #ifndef HAVE_ATAN2F
+ float
+ atan2f(float x, float y)
+ {
+   return (float) atan2(x, y);
+ }
+ #endif
+
+ #ifndef HAVE_ATANF
+ float
+ atanf(float x)
+ {
+   return (float) atan(x);
+ }
+ #endif
+
+ #ifndef HAVE_CEILF
+ float
+ ceilf(float x)
+ {
+   return (float) ceil(x);
+ }
+ #endif
+
+ #ifndef HAVE_COPYSIGNF
+ float
+ copysignf(float x, float y)
+ {
+   return (float) copysignf(x, y);
+ }
+ #endif
+
+ #ifndef HAVE_COSF
+ float
+ cosf(float x)
+ {
+   return (float) cos(x);
+ }
+ #endif
+
+ #ifndef HAVE_COSHF
+ float
+ coshf(float x)
+ {
+   return (float) cosh(x);
+ }
+ #endif
+
+ #ifndef HAVE_EXPF
+ float
+ expf(float x)
+ {
+   return (float) exp(x);
+ }
+ #endif
+
+ #ifndef HAVE_FLOORF
+ float
+ floorf(float x)
+ {
+   return (float) floor(x);
+ }
+ #endif
+
+ #ifndef HAVE_FREXPF
+ float
+ frexpf(float x, int *exp)
+ {
+   return (float) frexp(x, exp);
+ }
+ #endif
+
+ #ifndef HAVE_HYPOTF
+ float
+ hypotf(float x, float y)
+ {
+   return (float) hypot(x, y);
+ }
+ #endif
+
+ #ifndef HAVE_LOGF
+ float
+ logf(float x)
+ {
+   return (float) log(x);
+ }
+ #endif
+
+ #ifndef HAVE_LOG10F
+ float
+ log10f(float x)
+ {
+   return (float) log10(x);
+ }
+ #endif
+
+ #ifndef HAVE_SCALBNF
+ float
+ scalbnf(float x, int y)
+ {
+   return (float) scalbn(x, y);
+ }
+ #endif
+
+ #ifndef HAVE_SINF
+ float
+ sinf(float x)
+ {
+   return (float) sin(x);
+ }
+ #endif
+
+ #ifndef HAVE_SINHF
+ float
+ sinhf(float x)
+ {
+   return (float) sinh(x);
+ }
+ #endif
+
+ #ifndef HAVE_SQRTF
+ float
+ sqrtf(float x)
+ {
+   return (float) sqrt(x);
+ }
+ #endif
+
+ #ifndef HAVE_TANF
+ float
+ tanf(float x)
+ {
+   return (float) tan(x);
+ }
+ #endif
+
+ #ifndef HAVE_TANHF
+ float
+ tanhf(float x)
+ {
+   return (float) tanh(x);
+ }
+ #endif
+
+ #ifndef HAVE_NEXTAFTERF
+ /* This is a portable implementation of nextafterf that is intended to be
+    independent of the floating point format or its in memory representation.
+    This implementation skips denormalized values, for example returning
+    FLT_MIN as the next value after zero, as many target's frexpf, scalbnf
+    and ldexpf functions don't work as expected with denormalized values.  */
+ float
+ nextafterf(float x, float y)
+ {
+   int origexp, newexp;
+
+   if (isnan(x) || isnan(y))
+     return x+y;
+   if (x == y)
+     return x;
+
+   if (x == 0.0f)
+     return y > 0.0f ? FLT_MIN : -FLT_MIN;
+
+   frexpf(x, &origexp);
+   if (x >= 0.0)
+     {
+       if (y > x)
+ 	{
+ 	  if (x < FLT_MIN)
+ 	    return FLT_MIN;
+ 	  return x + scalbnf(FLT_EPSILON, origexp-1);
+ 	}
+       else if (x > FLT_MIN)
+ 	{
+ 	  float temp = x - scalbnf(FLT_EPSILON, origexp-1);
+ 	  frexpf(temp, &newexp);
+ 	  if (newexp == origexp)
+ 	    return temp;
+ 	  return x - scalbnf(FLT_EPSILON, origexp-2);
+ 	}
+       else
+ 	return 0.0f;
+     }
+   else
+     {
+       if (y < x)
+ 	{
+ 	  if (x > -FLT_MIN)
+ 	    return -FLT_MIN;
+ 	  return x - scalbnf(FLT_EPSILON, origexp-1);
+ 	}
+       else if (x < -FLT_MIN)
+ 	{
+ 	  float temp = x + scalbnf(FLT_EPSILON, origexp-1);
+ 	  frexpf(temp, &newexp);
+ 	  if (newexp == origexp)
+ 	    return temp;
+ 	  return x + scalbnf(FLT_EPSILON, origexp-2);
+ 	}
+       else
+ 	return 0.0f;
+     }
+ }
+ #endif
+
  /* Note that if HAVE_FPCLASSIFY is not defined, then NaN is not handled */

  /* Algorithm by Steven G. Kargl.  */


Roger
--
Roger Sayle,                         E-mail: roger@eyesopen.com
OpenEye Scientific Software,         WWW: http://www.eyesopen.com/
Suite 1107, 3600 Cerrillos Road,     Tel: (+1) 505-473-7385
Santa Fe, New Mexico, 87507.         Fax: (+1) 505-473-0833

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Paul Brook]

> 2004-08-02  Roger Sayle  <roger@eyesopen.com>
>
> 	PR libgfortran/16137
> 	* configure.ac: Add tests for acosf, asinf, atan2f, atanf, ceilf,
> 	copysignf, cosf, coshf, expf, floorf, frexpf, hypotf, logf, log10f,
> 	scalbnf, sinf, sinhf, sqrtf, tanf and tanhf in libm.
> 	* config.h.in: Regenerate.
> 	* configure: Regenerate.
>
> 	* instrinsics/c99_functions.c (acosf, asinf, atan2f, atanf, ceilf,
> 	copysignf, cosf, coshf, expf, floorf, frexpf, hypotf, logf, log10f,
> 	nextafterf, scalbnf, sinf, sinhf, sqrtf, tanf, tanhf): New stub
> 	implementations for targets that don't support C99 float functions.

Ok, with the typo fix pointer out by Arnaud.

Paul

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Richard Henderson]

[-- Attachment #1: Type: text/plain, Size: 456 bytes --]

On Mon, Aug 02, 2004 at 08:34:34PM -0600, Roger Sayle wrote:
> Hopefully, gfortran's use of nextafterf doesn't depend upon
> its operation with denormalized numbers.

Well, actually, from my reading the NEAREST intrinsic *should*
in fact operate correctly with denormals, and that

	NEAREST(0.0, 1.0)

should yield the minimal denormal.

You could see if the following works for your target.  It passes
the attached tests vs the glibc implementation.


r~

[-- Attachment #2: z.c --]
[-- Type: text/plain, Size: 2061 bytes --]

#include <math.h>
#include <float.h>

float
mynextafterf(float x, float y)
{
  float absx, delta, denorm_min;
                                                                                
  if (isnan(x) || isnan(y))
    return x+y;
  if (x == y)
    return x;

  absx = fabsf (x);

  /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
  if (__FLT_DENORM_MIN__ == 0.0f)
    denorm_min = FLT_MIN;
  else
    denorm_min = __FLT_DENORM_MIN__;

  if (absx < FLT_MIN)
    delta = denorm_min;
  else
    {
      float frac;
      int exp;
                                                                                
      /* Discard the fraction from x.  */
      frac = frexpf (absx, &exp);
      delta = scalbnf (0.5f, exp);
                                                                                
      /* Scale x by the epsilon of the representation.  By rights we should
	 have been able to combine this with scalbnf, but some targets don't
	 get that correct with denormals.  */
      delta *= FLT_EPSILON;

      /* If we're going to be reducing the absolute value of X, and doing so
	 would reduce the exponent of X, then the delta to be applied is
	 one exponent smaller.  */
      if (frac == 0.5f && (y < x) == (x > 0))
	delta *= 0.5f;

      /* If that underflows to zero, then we're back to the minimum.  */
      if (delta == 0.0f)
	delta = denorm_min;
    }

  if (y < x)
    delta = -delta;

  return x + delta;
}

void test(float x, float y)
{
  float r1, r2;

  r1 = nextafterf(x, y);
  r2 = mynextafterf(x, y);

  if (r1 != r2)
    abort ();
}

int main()
{
  test (0, 1);
  test (0, -1);
  test (1, 2);
  test (1, -2);
  test (-1, 2);
  test (-1, -2);
  test (__FLT_DENORM_MIN__, 1);
  test (-__FLT_DENORM_MIN__, 1);
  test (4*__FLT_DENORM_MIN__, 1);
  test (4*__FLT_DENORM_MIN__, -1);
  test (FLT_MIN, 1);
  test (FLT_MIN, -1);
  test (4*FLT_MIN, 1);
  test (4*FLT_MIN, -1);
  test (8388608*FLT_MIN, 1);
  test (8388608*FLT_MIN, -1);
  test (16777216*FLT_MIN, 1);
  test (16777216*FLT_MIN, -1);

  return 0;
}

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Toon Moene]

Richard Henderson wrote:

> On Mon, Aug 02, 2004 at 08:34:34PM -0600, Roger Sayle wrote:
> 
>>Hopefully, gfortran's use of nextafterf doesn't depend upon
>>its operation with denormalized numbers.

> Well, actually, from my reading the NEAREST intrinsic *should*
> in fact operate correctly with denormals, and that
> 
> 	NEAREST(0.0, 1.0)
> 
> should yield the minimal denormal.

I put this question up to J3.  I think the answer might not be the 
obvious, although the most senior member in charge of floating point 
behaviour differs with my interpretation, especially with regard to the 
question whether a processor that sqaushes all denormals to zero is 
standard conforming (I think it is, he thinks it isn't).

I'll be off-line for J3 meeting 169 as of tomorrow (7th of August) 8:00 
CEST until the 15th of August, but I'll have a definite answer after 
that :-)

Cheers,

-- 
Toon Moene - mailto:toon@moene.indiv.nluug.nl - phoneto: +31 346 214290
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands
Maintainer, GNU Fortran 77: http://gcc.gnu.org/onlinedocs/g77_news.html
A maintainer of GNU Fortran 95: http://gcc.gnu.org/fortran/

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Roger Sayle]

On Tue, 3 Aug 2004, Richard Henderson wrote:
> On Mon, Aug 02, 2004 at 08:34:34PM -0600, Roger Sayle wrote:
> > Hopefully, gfortran's use of nextafterf doesn't depend upon
> > its operation with denormalized numbers.
>
> Well, actually, from my reading the NEAREST intrinsic *should*
> in fact operate correctly with denormals, and that
>
> 	NEAREST(0.0, 1.0)
>
> should yield the minimal denormal.
>
> You could see if the following works for your target.  It passes
> the attached tests vs the glibc implementation.

Sorry for taking so long to reply, but it's taken me a few days to
find the time to realize that the aborts that I'm seeing on x86
when compiling with optimization are caused by the "excess precision"
problem.  Doh!  Indeed compiling your function with -ffloat-store or
-O0 allows the tests to pass on i686-pc-linux-gnu.

The following tweak to your function resolves the problem at the
source :>


  float
  mynextafterf(float x, float y)
  {
-   float absx, delta, denorm_min;
+   float absx, denorm_min;
+   volatile float delta;

    if (isnan(x) || isnan(y))
      return x+y;


Would this be acceptable on all platforms, or would it be better to
#ifdef the "volatile" for just IA-32?

Roger
--

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Paul Brook]

On Friday 06 August 2004 21:52, Roger Sayle wrote:
> Sorry for taking so long to reply, but it's taken me a few days to
> find the time to realize that the aborts that I'm seeing on x86
> when compiling with optimization are caused by the "excess precision"
> problem.  Doh!  Indeed compiling your function with -ffloat-store or
> -O0 allows the tests to pass on i686-pc-linux-gnu.
>
> The following tweak to your function resolves the problem at the
> source :>
>
>
>   float
>   mynextafterf(float x, float y)
>   {
> -   float absx, delta, denorm_min;
> +   float absx, denorm_min;
> +   volatile float delta;
>
>     if (isnan(x) || isnan(y))
>       return x+y;
>
>
> Would this be acceptable on all platforms, or would it be better to
> #ifdef the "volatile" for just IA-32?

Ok with a comment saying why it's marked volatile. I've done something similar 
elsewhere.

Paul

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Richard Henderson]

On Fri, Aug 06, 2004 at 10:30:25PM +0100, Paul Brook wrote:
> > Would this be acceptable on all platforms, or would it be better to
> > #ifdef the "volatile" for just IA-32?
> 
> Ok with a comment saying why it's marked volatile. I've done something
> similar elsewhere.

A better solution is

static inline float
force_round_float (float x)
{
  if (!__builtin_constant_p (x) && FLT_EVAL_METHOD != 0)
    asm ("" : "=m"(x) : "m"(x));
  return x;
}

used as needed.  Volatile on "delta" doesn't actually address the
situation at the correct point, I don't think.

We also don't pay the penalty on well-behaved targets.



r~

Re: [libgfortran] PR 16137: libgfortran on non-C99 targets

[Richard Henderson]

On Fri, Aug 06, 2004 at 10:04:47PM +0200, Toon Moene wrote:
> I put this question up to J3.  I think the answer might not be the 
> obvious, although the most senior member in charge of floating point 
> behaviour differs with my interpretation, especially with regard to the 
> question whether a processor that sqaushes all denormals to zero is 
> standard conforming (I think it is, he thinks it isn't).

Seems to me that if the system reliably squashes all denormals, that
that is a consistent machine model for the arithmetic, and so is 
conforming.

The problem would come when you've got a bit to flip at runtime,
and it only eliminates *some* of the denormals.  E.g. this NEAREST
thing, where one might do the computation in integer modes, and so
generate a number that you couldn't generate with actual fp 
arithmetic.


r~