[PATCH/gfortran] Chop bits of subnormal numbers

[Steve Kargl]

The attached is a brute force implementation of the chop
rounding recommended by rth for subnormal numbers in the
MPFR representation.  A future version of MPFR will have
a more efficient and elegant function to perform this 
operation.

2005-05-07  Steven G. Kargl  <kargls@comcast.net>

	* arith.c (gfc_check_integer_range): zero insignificant bits
	in significands of subnormal numbers

-- 
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	7 May 2005 21:14:36 -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,414 ----
    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, p, 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 (p, s, gfc_real_kinds[i].radix, GMP_RNDN);
+       gfc_free(s);
+       gfc_free(bin);
+ 
+       retval = ARITH_OK;
+     }
    else
      retval = ARITH_OK;