Bug 43210 - frexpl() returns unnormal numbers when given denormal input
frexpl() returns unnormal numbers when given denormal input
Status: CLOSED RAWHIDE
Product: Red Hat Linux
Classification: Retired
Component: libc (Show other bugs)
7.1
i386 Linux
low Severity low
: ---
: ---
Assigned To: Jakub Jelinek
:
Depends On:
Blocks:
  Show dependency treegraph
 
Reported: 2001-06-01 16:31 EDT by Trevin Beattie
Modified: 2008-05-01 11:38 EDT (History)
0 users

See Also:
Fixed In Version:
Doc Type: Bug Fix
Doc Text:
Story Points: ---
Clone Of:
Environment:
Last Closed: 2001-06-02 15:13:37 EDT
Type: ---
Regression: ---
Mount Type: ---
Documentation: ---
CRM:
Verified Versions:
Category: ---
oVirt Team: ---
RHEL 7.3 requirements from Atomic Host:
Cloudforms Team: ---


Attachments (Terms of Use)

  None (edit)
Description Trevin Beattie 2001-06-01 16:31:34 EDT
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:
Comment 1 Jakub Jelinek 2001-06-01 17:31:12 EDT
Apparently this is wrong only in the i386 assembly frexpl, the C version
works fine. Will fix it soonish.
Comment 2 Trevin Beattie 2001-06-01 18:06:59 EDT
I've looked at the assembly code and compared it with the double and float
versions, and it looks like the fix is exceedingly simple.  All I have to do is
remove some excess code.  Here's the patch for the Gnu glibc distribution
(should be the same one Linux is using):

*** glibc-2.2.3/sysdeps/i386/fpu/s_frexpl.S.orig        Sat Jan  6 20:35:39 2001
--- glibc-2.2.3/sysdeps/i386/fpu/s_frexpl.S     Fri Jun  1 14:57:10
2001***************
*** 32,43 ****
        ASM_TYPE_DIRECTIVE(two64,@object)
  two64:        .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
        ASM_SIZE_DIRECTIVE(two64)
-       /* The following is LDBL_MAX / ldexp (1.0, 64), the largest
-          number we can handle the normal way.  */
-       ASM_TYPE_DIRECTIVE(largest,@object)
- largest:
-       .byte 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xbe, 0x7f, 0, 0
-       ASM_SIZE_DIRECTIVE(largest)
  
  #ifdef PIC
  #define MO(op) op##@GOTOFF(%edx)
--- 32,37 ----
***************
*** 66,76 ****
        cmpl    $0x7fff, %eax
        je      1f
  
!       cmpl    $0, %eax
!       je      2f
! 
!       cmpl    $0x7fbe, %eax
!       ja      4f
  
        fldt    VAL0(%esp)
  #ifdef        PIC
--- 60,67 ----
        cmpl    $0x7fff, %eax
        je      1f
  
!       cmpl    $1, %eax
!       jae     2f
  
        fldt    VAL0(%esp)
  #ifdef        PIC
***************
*** 102,116 ****
  
        LEAVE
        ret
- 
- 4:    movl    VAL2(%esp), %ecx
-       movl    %ecx, %edx
-       andl    $0x7fff, %ecx
- 
-       andl    $0x8000, %edx
-       subl    $16382, %ecx
-       orl     $0x3ffe, %edx
-       movl    %edx, VAL2(%esp)
-       jmp     1b
  END (BP_SYM (__frexpl))
  weak_alias (BP_SYM (__frexpl), BP_SYM (frexpl))
--- 93,97 ----


I've tested this in my own version of libc (which I'm porting to another o/s)
and it works great--handles all extreme cases as expected.

BTW, my test programs have uncovered a few other lesser bugs in glibc's math
functions.  Any interest?
Comment 3 Jakub Jelinek 2001-06-02 03:58:42 EDT
We are certainly interested in all bugs you find, please fill
them in.
As for this one, I'll try to check this during the weekend,
gotta leave now.
Comment 4 Jakub Jelinek 2001-06-02 15:13:33 EDT
Agreed with one minor change, see
http://sources.redhat.com/ml/libc-hacker/2001-06/msg00000.html
I've added a testcase for it too.
Comment 5 Jakub Jelinek 2001-06-07 09:17:04 EDT
Included in glibc-2.2.3-11.

Note You need to log in before you can comment on or make changes to this bug.