From Bugzilla Helper:
User-Agent: Mozilla/4.77 [en] (X11; U; Linux 2.4.2-2 i586)

Description of problem:
If you pass a denormal number to frexpl() (i.e., a 'long double' in the
range 2^-16383 to 2^-16445), frexpl() will convert the number to an
unnormal format by scaling the exponent without changing the mantissa, and
the stored exponent will always be -16382.

An "unnormal" format is one in which the exponent is non-zero, but the
integer bit is zero.  The FPU never generates unnormal numbers.  If an
unnormal number is used as the operand in a floating-point instruction, the
FPU will treat it as Not-a-Number.


How reproducible:
Always

Steps to Reproduce:
Compile and run the following code:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

main ()
{
  union {
    long double ld;
    struct {
      unsigned long long f:63;
      unsigned long i:1;
      unsigned short e:15;
      unsigned short s:1;
    } part;
  } x, y;
  int m, n;

  srand (time (NULL));

  /* For comparison, start with a normalized number. */
  m = 31 - rand () % 63;
  x.ld = 1.0L;
  y.ld = 1.0L;
  x.part.e += m;	/* Scale X by m */
  printf ("2^%d = %Lg = [ %c %d.%016llx *2^%d ]\n", m, x.ld,
	  x.part.s ? '-' : '+', x.part.i, x.part.f * 2,
	  (signed int) x.part.e - 0x3fff);

  y.ld = frexpl (x.ld, &n);
  printf ("frexpl(%Lg) = { %Lg [ %c %d.%016llx *2^%d ], %d }\n", x.ld,
	  y.ld, y.part.s ? '-' : '+', y.part.i, y.part.f * 2,
	  (signed int) y.part.e - 0x3fff, n);

  m = -16383 - rand () % 63;
  x.ld = 1.0L;
  y.ld = 1.0L;
  x.part.e += -8192;	/* Scale X by -8192 */
  y.part.e += m + 8192; /* Scale Y by the rest of m */
  x.ld = x.ld * y.ld;	/* Scale X down to -16383..-16445 (denormal) */
  printf ("2^%d = %Lg = [ %c %d.%016llx *2^%d ]\n", m, x.ld,
	  x.part.s ? '-' : '+', x.part.i, x.part.f * 2,
	  (signed int) x.part.e - 0x3fff);

  y.ld = frexpl (x.ld, &n);
  printf ("frexpl(%Lg) = { %Lg [ %c %d.%016llx *2^%d ], %d }\n", x.ld,
	  y.ld, y.part.s ? '-' : '+', y.part.i, y.part.f * 2,
	  (signed int) y.part.e - 0x3fff, n);
}


Actual Results:  2^-14 = 6.10352e-05 = [ + 1.0000000000000000 *2^-14 ]
frexpl(6.10352e-05) = { 0.5 [ + 1.0000000000000000 *2^-1 ], -13 }
2^-16424 = 7.64454e-4945 = [ + 0.0000000000400000 *2^-16383 ]
frexpl(7.64454e-4945) = { 0.000000000000113687 [ + 0.0000000000400000 *2^-1
], -16382 }


Expected Results:  The last line should show:
frexpl(7.64454e-4945) = { 0.5 [ + 1.0000000000000000 *2^-1 ], -16423 }

According to WG14/N869 7.12.6.4, frexp must return a normalized value X and
an integer exponent EXP, such that X has a magnitude in the interval [1/2,
1) or zero, and X*2^EXP equals the original value.


Additional info: