HUGE() pain in the ...

[Steve Kargl]

How I spent my Saturday.  Well, I started to go over my MPFR
patch and run various checks.  The first one involved the 
HUGE intrinsic, and I am perplexed.

Consider

  program z
    implicit none
    real(4) s, sb
    real(8) d, db
    sb = radix(s)
    s = (1.e0_4 - sb**(-digits(s)))*sb**(maxexponent(s)-1)*sb
    print *, 'huge(s) = ', huge(s), s
    db = radix(d)
    d = (1.e0_8 - db**(-digits(d)))*db**(maxexponent(d)-1)*db
    print *, 'huge(d) = ', huge(d), d
  end program z

kargl[221] gfc -o z z.f90
3.402823e38               <-- Debugging output from arith.c
1.79769313486232e308      <-- using mpfr_output_str.
kargl[222] ./z
 huge(s) =   3.4028230E+38      3.4028234E+38    
 huge(d) =                +Infinity  1.797693134862316E+308     

Okay, what do we see.  The first column should be huge() as 
computed in arith.c and the 2nd column is from the computation
  huge(x) = (1 - b**(-p)) * b**emax.
          = (1 - b**(-p)) * b**(emax-1) * b.
The single precision computation exceeds huge(s), and yet we
print 3.4028234E+38 instead of +Infinity.  Second, if we
output precision(s), we find a value of 6, and yet we are
print 8 decimal digits (mpfr_output_str gives 7 decimal digits).  
In double precision we get +Infinity printed for huge(d) and
the computation appears valid in comparison to mpfr_output_str
output.  Note, precision(d) gives 15 decimal digits, we print
16 decimal digits and mpfr_output_str yields 15 decimal digits.

I not sure if this is a failure of MPFR or our output routine.
Anyone have any insight?


-- 
Steve