Re: [gfortran] I/O of large integer and real kinds, round 2

[Richard Henderson <rth at redhat dot com>]

On Thu, Jun 09, 2005 at 12:50:04AM +0200, FX Coudert wrote:
> +      while (x > 1e100L)
> +      {
> +        result += 100;
> +        x /= 1e100L;
> +      }
> +      while (x > 1e10L)
> +      {
> +        result += 10;
> +        x /= 1e10L;
> +      }

You don't have to work this hard.  IEEE double has a dynamic range of
1e308, or so.  Since we're falling back to log10 for evaluation of the
significand, there's no point in reducing the value below DBL_MAX.

Now, real(10) doesn't have *any* more dynamic range, just extra
precision.  So in that case you might as well pass on to log10
immediately and do nothing else.

For real(16), the dynamic range is 1e4932 or so.  In that case,
reducing by 1e100 is going to take quite a long time, and be remarkably
inaccurate.  While you could adjust your scheme to reduce by 1e1000,
it seems like 

long double log10l (long double x)
{
#if LDBL_MAX_EXP > DBL_MAX_EXP
  if (x > DBL_MAX)
    {
      int p2_result = 0;
      if (x >= 0x1p16384L) { p2_result += 16384; x /= 0x1p16384L; }
      if (x >= 0x1p8192L) { p2_result += 8192; x /= 0x1p8192L; }
      if (x >= 0x1p4096L) { p2_result += 4096; x /= 0x1p4096L; }
      if (x >= 0x1p2048L) { p2_result += 2048; x /= 0x1p2048L; }
      if (x >= 0x1p1024L) { p2_result += 1024; x /= 0x1p1024L; }
      return log10 (x) + p2_result * .30102999566398119521373889472449302L;
    }
  else
#endif
    return log10 (x);
}

would be a whole lot faster and be more accurate.  This takes into
account knowledge that real(16) has 4 bits more exponent than real(8),
and specifically targets those four bits for reduction.  The last
reduction (0x1p1024L) is so that something just a bit larger than
DBL_MAX (e.g. DBL_MAX+0x1p912L) doesn't get rounded up to +Inf.


r~