Re: [gfortran] Fix NEAREST constant folder

[Tobias Schlüter <tobias dot schlueter at physik dot uni-muenchen dot de>]

Richard Henderson wrote:
> On Tue, Apr 12, 2005 at 01:56:51PM -0700, Steve Kargl wrote:
> 
>>I suspect we will need to rip MPFR out of gfortran, and
>>use gcc/real.c.
> 
> I'm not thrilled about this, since there had been a vague push
> to move real.c to use MPFR as well.  What about creating some
> wrapper functions that allow us to fake subnormal numbers?

Here's a very first go at such a wrapper.  This is just a simple
self-contained program which outlines the possible interface, and I don't
think it will look much like this on inclusion, but I have a few questions and
thought it would be better to accompany them with some code (even if its debug
messages are in German and it may be off-by-one or contain similar bugs):

- mpfr returns the direction in which it has rounded.  Do we need this
information anywhere, what should fp_normalize return?  Currently it returns
one of a few magic numbers.
- the way it is this will probably slow arithmetic down considerably, because
I have to always extract both the mantissa and the exponent, even if we're not
dealing with a subnormal.  Is there a way to get the exponent without the
mantissa?
- other suggestions?

- Tobi

#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>

/* mpfr doesn't support denormals.  Since most targets use denormals
   in their arithmetic, we have to provide a means to calculate
   correctly rounded results even in those cases.  Hence we have to
   teach mpfr about denormals.

   What we do is very simple: the global variable emin contains the
   minimum exponent, prec the number of binary digits for the real
   kind that is currently operated with.  If the exponent exp of the
   result of a calculation is smaller than emin, we do away with the
   lower (emin - exp) bytes of the results mantissa.  

   In order to make sure that this happens consistently, wrappers for
   all mpfr functions are provided.  */

int prec;
int emin;

static int
fp_normalize (mpfr_t out, mpfr_t in)
{
  mpz_t mantissa;
  mp_exp_t exp;
  int shift, i;

  /* Special case zero.  */
  if (mpfr_sgn (in) == 0)
    {
      mpfr_set (out, in, GMP_RNDN);
      return 0;
    }	

  mpz_init (mantissa);
  exp = mpfr_get_z_exp (mantissa, in);

  gmp_printf ("Mantisse Start: 0x%Zx E%d\n", mantissa, exp);

  /* The returned exponent is (normalized exponent - prec), adjust for
     this.  */
  exp = exp + prec;
  if (exp >= emin)
    {
      mpz_clear (mantissa);
      mpfr_set (out, in, GMP_RNDN);
      return 2;
    }

  /* We have a subnormal number.  First determine if there are no
     digits to be kept.  */
  if (exp < emin - prec)
    {
      printf ("exp = %d\n", exp);
      mpfr_set_ui (out, 0, GMP_RNDN);
      mpz_clear (mantissa);

      return -1;
    }
  
  /* There are some digits that should remain, clear lowest bits.  */
  for (i=0; i < emin - exp; i++)
    mpz_clrbit (mantissa, i);
  
  gmp_printf ("Mantisse Ende: 0x%Zx\n", mantissa);

  mpfr_set_z (out, mantissa, GMP_RNDN);
  mpfr_print_binary (out);
  puts(" Binärergebnis vor shift");
  mpfr_mul_2si (out, out, exp - prec, GMP_RNDN);
  mpfr_print_binary (out);
  puts(" Binärergebnis nach shift");
  
  mpz_clear (mantissa);

  return 1;
}

#define fpfunc2(op, op2_t) \
   int \
   fp_ ## op (mpfr_t res, mpfr_t op1, op2_t op2, mp_rnd_t rnd) \
   { \
     mpfr_ ## op (res, op1, op2, rnd); \
     mpfr_print_binary (res); puts (""); \
     return fp_normalize (res, res); \
   }

fpfunc2(add, mpfr_t);
fpfunc2(add_ui, unsigned long int);
fpfunc2(sub, mpfr_t);
fpfunc2(mul, mpfr_t);
fpfunc2(div, mpfr_t);
fpfunc2(pow_ui, unsigned long int);
fpfunc2(pow, mpfr_t);
fpfunc2(mul_2ui, unsigned long int);
fpfunc2(mul_2si, signed long int);
/* etc.  */

int
main ()
{
  mpfr_t x, y;

  prec = 5;
  mpfr_set_default_prec (prec);
  emin = -1;

  mpfr_init_set_str (x, ".125", 10, GMP_RNDN);

  mpfr_init_set_str (y, ".375", 10, GMP_RNDN);
  mpfr_div_ui (y, y, 16, GMP_RNDN);

  printf ("%d\n", fp_add (x, x, y, GMP_RNDN));

  mpfr_print_binary (x);
  puts("");

  return 0;
}