[PATCH/gfortran] Chop bits of subnormal numbers

[Steve Kargl]

On Sun, May 08, 2005 at 12:31:44AM +0200, Tobias Schl?ter wrote:
> 
> This won't work for negative values of p.  Apart from this, this looks fine.
> 

Damn details. ;-)

I've updated the patch to chop the bits on the absolute
value, and then set the sign if need be.

Bubblestrapped and reg. tested onn i386-*-freebsd.

OK for mainline?

-- 
Steve
-------------- next part --------------
Index: arith.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/fortran/arith.c,v
retrieving revision 1.26
diff -c -p -r1.26 arith.c
*** arith.c	14 Apr 2005 16:29:31 -0000	1.26
--- arith.c	8 May 2005 16:23:01 -0000
*************** gfc_check_integer_range (mpz_t p, int ki
*** 359,366 ****
  
    return result;
  }
- 
- 
  /* Given a real and a kind, make sure that the real lies within the
     range of the kind.  Returns ARITH_OK, ARITH_OVERFLOW or
     ARITH_UNDERFLOW.  */
--- 359,364 ----
*************** gfc_check_real_range (mpfr_t p, int kind
*** 381,389 ****
    if (mpfr_sgn (q) == 0)
      retval = ARITH_OK;
    else if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
!       retval = ARITH_OVERFLOW;
    else if (mpfr_cmp (q, gfc_real_kinds[i].subnormal) < 0)
      retval = ARITH_UNDERFLOW;
    else
      retval = ARITH_OK;
  
--- 379,420 ----
    if (mpfr_sgn (q) == 0)
      retval = ARITH_OK;
    else if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
!     retval = ARITH_OVERFLOW;
    else if (mpfr_cmp (q, gfc_real_kinds[i].subnormal) < 0)
      retval = ARITH_UNDERFLOW;
+   else if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
+     {
+       /* MPFR operates on a numbers with a given precision and enormous
+ 	exponential range.  To represent subnormal numbers the exponent is
+ 	allowed to become smaller than emin, but always retains the full
+ 	precision.  This function resets unused bits to 0 to alleviate
+ 	rounding problems.  Note, a future version of MPFR will have a
+  	mpfr_subnormalize() function, which handles this truncation in a
+ 	more efficient and robust way.  */
+ 
+       int j, k;
+       char *bin, *s;
+       mp_exp_t e;
+ 
+       bin = mpfr_get_str (NULL, &e, gfc_real_kinds[i].radix, 0, q, GMP_RNDN);
+       k = gfc_real_kinds[i].digits - (gfc_real_kinds[i].min_exponent - e);
+       for (j = k; j < gfc_real_kinds[i].digits; j++)
+ 	bin[j] = '0';
+       /* Need space for '0.', bin, 'E', and e */
+       s = (char *) gfc_getmem (strlen(bin)+10);
+       sprintf (s, "0.%sE%d", bin, (int) e);
+       mpfr_set_str (q, s, gfc_real_kinds[i].radix, GMP_RNDN);
+ 
+       if (mpfr_sgn (p) < 0)
+ 	mpfr_neg (p, q, GMP_RNDN);
+       else
+ 	mpfr_set (p, q, GMP_RNDN);
+ 
+       gfc_free (s);
+       gfc_free (bin);
+ 
+       retval = ARITH_OK;
+     }
    else
      retval = ARITH_OK;