>Number:         32861
>Category:       lib
>Synopsis:       frexp returns incorrect value for small floating-point values
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    lib-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Fri Feb 17 11:40:00 +0000 2006
>Originator:     FX Coudert
>Release:        3.0
>Organization:
>Environment:
NetBSD td143.testdrive.hp.com 3.0 NetBSD 3.0 (GENERIC.MP) #0: Mon Dec 19 01:19:07 UTC 2005  builds@works.netbsd.org:/home/builds/ab/netbsd-3-0-RELEASE/i386/200512182024Z-obj/home/builds/ab/netbsd-3-0-RELEASE/src/sys/arch/i386/compile/GENERIC.MP i386
>Description:
Calling frexp on small values (like double d = 4.94066e-324) gives incorrect results for the exponent. With the example from the NetBSD frexp man page (but with a modified value of d):

frexp(4.94066e-324, &e) = 0.5: 0.5 * 2^-1022 = 4.94066e-324

while the correct answer is (obtained on i686-linux, eg):

frexp(4.94066e-324, &e) = 0.5: 0.5 * 2^-1073 = 4.94066e-324

The non-correctness of the NetBSD result can also be checked with any scientific calculator at hand.
>How-To-Repeat:
$ cat b.c
#include <stdio.h>
#include <math.h>
#include <float.h>
int main () {
  double d = 4.94066e-324;
  int e;
  double f = frexp(d, &e);
  printf("frexp(%g, &e) = %g: %g * %d^%d = %g\n", d, f, f, FLT_RADIX, e, d);
  return 0;
}

$ gcc b.c -lm && ./a.out
frexp(4.94066e-324, &e) = 0.5: 0.5 * 2^-1022 = 4.94066e-324
>Fix: