NetBSD Problem Report #50698

From gson@gson.org  Sat Jan 23 22:40:04 2016
Return-Path: <gson@gson.org>
Received: from mail.netbsd.org (mail.NetBSD.org [199.233.217.200])
	(using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits))
	(Client CN "mail.netbsd.org", Issuer "Postmaster NetBSD.org" (verified OK))
	by mollari.NetBSD.org (Postfix) with ESMTPS id B67C27ABF5
	for <gnats-bugs@gnats.NetBSD.org>; Sat, 23 Jan 2016 22:40:04 +0000 (UTC)
Message-Id: <20160123223955.AB49674508B@guava.gson.org>
Date: Sun, 24 Jan 2016 00:39:55 +0200 (EET)
From: gson@gson.org (Andreas Gustafsson)
Reply-To: gson@gson.org (Andreas Gustafsson)
To: gnats-bugs@gnats.NetBSD.org
Subject: hypotf() of small numbers yields infinity
X-Send-Pr-Version: 3.95

>Number:         50698
>Category:       lib
>Synopsis:       hypotf() of small numbers yields infinity
>Confidential:   no
>Severity:       serious
>Priority:       high
>Responsible:    lib-bug-people
>State:          closed
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Sat Jan 23 22:45:00 +0000 2016
>Closed-Date:    Fri Mar 11 15:16:39 +0000 2016
>Last-Modified:  Fri Mar 11 15:16:39 +0000 2016
>Originator:     Andreas Gustafsson
>Release:        NetBSD 7.0
>Organization:

>Environment:
System: NetBSD
Architecture: x86_64
Machine: amd64
>Description:

The hypotf() function in libm incorrectly returns infinity (!) when
its arguments are small numbers such as 1e-18.

I have tested this on amd64 (7.0, real hardware), i386 (-current,
qemu), and sparc (-current, qemu), and they all behave this way.

>How-To-Repeat:

$ cat test-hypotf.c
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv) {
    float a = atof(argv[1]);
    float b = atof(argv[2]);
    float c = hypotf(a, b);
    printf("%.15g\n", c);
    return 0;
}
$ cc test-hypotf.c -lm -o test-hypotf
$ ./test-hypotf 1e-18 1e-18
inf

On a non-NetBSD-machine, the test program gives the correct result:

./a.out 1e-18 1e-18
1.41421363268178e-18

>Fix:

>Release-Note:

>Audit-Trail:
From: Rin Okuyama <okuyama@flex.phys.tohoku.ac.jp>
To: gnats-bugs@NetBSD.org
Cc: 
Subject: Re: lib/50698: hypotf() of small numbers yields infinity
Date: Sun, 24 Jan 2016 19:26:06 +0900

 This is because some magic numbers in e_hypotf.c are wrong. In FreeBSD,
 this bug was fixed by commit #23397:

   https://svnweb.freebsd.org/base/head/lib/msun/src/e_hypotf.c#rev23397

 I attached the minimal patch below. With this patch, I obtain (on amd64)

   % ./test-hypotf 1e-18 1e-18
   1.41421363268178e-18
   %

 Other changes have been made for libm in FreeBSD. IMO, we should merge
 reasonable ones.

 --- src/lib/libm/src/e_hypotf.c.orig	2016-01-24 18:40:38.426417738 +0900
 +++ src/lib/libm/src/e_hypotf.c	2016-01-24 19:07:19.114100741 +0900
 @@ -43,22 +43,22 @@
  	       if(hb == 0x7f800000) w = b;
  	       return w;
  	   }
 -	   /* scale a and b by 2**-60 */
 -	   ha -= 0x5d800000; hb -= 0x5d800000;	k += 60;
 +	   /* scale a and b by 2**-68 */
 +	   ha -= 0x22000000; hb -= 0x22000000;	k += 68;
  	   SET_FLOAT_WORD(a,ha);
  	   SET_FLOAT_WORD(b,hb);
  	}
  	if(hb < 0x26800000) {	/* b < 2**-50 */
  	    if(hb <= 0x007fffff) {	/* subnormal b or 0 */
  	        if(hb==0) return a;
 -		SET_FLOAT_WORD(t1,0x3f000000);	/* t1=2^126 */
 +		SET_FLOAT_WORD(t1,0x7e800000);	/* t1=2^126 */
  		b *= t1;
  		a *= t1;
  		k -= 126;
 -	    } else {		/* scale a and b by 2^60 */
 -	        ha += 0x5d800000; 	/* a *= 2^60 */
 -		hb += 0x5d800000;	/* b *= 2^60 */
 -		k -= 60;
 +	    } else {		/* scale a and b by 2^68 */
 +	        ha += 0x22000000; 	/* a *= 2^68 */
 +		hb += 0x22000000;	/* b *= 2^68 */
 +		k -= 68;
  		SET_FLOAT_WORD(a,ha);
  		SET_FLOAT_WORD(b,hb);
  	    }

From: "Andreas Gustafsson" <gson@netbsd.org>
To: gnats-bugs@gnats.NetBSD.org
Cc: 
Subject: PR/50698 CVS commit: src/lib/libm/src
Date: Sun, 24 Jan 2016 20:34:31 +0000

 Module Name:	src
 Committed By:	gson
 Date:		Sun Jan 24 20:34:30 UTC 2016

 Modified Files:
 	src/lib/libm/src: e_hypotf.c

 Log Message:
 Fix incorrect magic numbers in scaling.  From FreeBSD commit 23397, by
 way of Rin Okuyama.  Fixes PR lib/50698.


 To generate a diff of this commit:
 cvs rdiff -u -r1.9 -r1.10 src/lib/libm/src/e_hypotf.c

 Please note that diffs are not public domain; they are subject to the
 copyright notices on the relevant files.

From: Andreas Gustafsson <gson@gson.org>
To: Rin Okuyama <okuyama@flex.phys.tohoku.ac.jp>
Cc: gnats-bugs@NetBSD.org
Subject: Re: lib/50698: hypotf() of small numbers yields infinity
Date: Mon, 25 Jan 2016 17:01:25 +0200

 Rin Okuyama wrote:
 >  This is because some magic numbers in e_hypotf.c are wrong. In FreeBSD,
 >  this bug was fixed by commit #23397:
 >  
 >    https://svnweb.freebsd.org/base/head/lib/msun/src/e_hypotf.c#rev23397
 >  
 >  I attached the minimal patch below. With this patch, I obtain (on amd64)
 >  
 >    % ./test-hypotf 1e-18 1e-18
 >    1.41421363268178e-18
 >    %

 Thank you for your analysis and patch.  The patch has now been committed.

 >  Other changes have been made for libm in FreeBSD. IMO, we should merge
 >  reasonable ones.

 Additional patches for the ones you think are reasonable are welcome.
 More test cases would be good, too.
 -- 
 Andreas Gustafsson, gson@gson.org

From: Rin Okuyama <okuyama@flex.phys.tohoku.ac.jp>
To: Andreas Gustafsson <gson@gson.org>
Cc: gnats-bugs@NetBSD.org
Subject: Re: lib/50698: hypotf() of small numbers yields infinity
Date: Mon, 1 Feb 2016 22:01:38 +0900

 On 2016/01/26 0:01, Andreas Gustafsson wrote:
 >>   Other changes have been made for libm in FreeBSD. IMO, we should merge
 >>   reasonable ones.
 >
 > Additional patches for the ones you think are reasonable are welcome.
 > More test cases would be good, too.

 Sorry for the late reply. I examined FreeBSD's libm, and found it much
 harder than I had imagined to merge their changes into ours.

 There are three kinds of differences: (1) trivial bug fixes as in this
 case, (2) long double functions missing in our libm, and (3) non-trivial
 changes to improve the accuracy or processing speed.

 We can readily merge (1) by patiently comparing the two implementations.

 For (2), we have not fully implemented long double functions. Instead,
 we have wrappers to double functions in src/lib/libm/src/ldbl_dummy.c.
 Needless to say, this is problematic for architectures with
 LDBL_MANT_DIG != DBL_MANT_DIG. They have long double functions both for
 LDBL_MANT_DIG == 64 (amd64, i386, ia64, and m68k), and 113 (aarch64,
 mips64, and sparc64). If we can merge it, we then have long double
 functions for all architectures except powerpc64 (LDBL_MANT_DIG == 106).

 However, for (2) and (3), there remains a big problem; how can we test
 them? Addition to ATF-based tests ported from us, they have their own
 test cases:

    https://svnweb.freebsd.org/base/head/lib/msun/tests/

 Function returns for some special arguments are checked in their test.
 But is it enough? Who can check hundreds of magic numbers?

State-Changed-From-To: open->pending-pullups
State-Changed-By: gson@NetBSD.org
State-Changed-When: Mon, 15 Feb 2016 16:20:46 +0000
State-Changed-Why:
pullup-7 #1120


From: "Martin Husemann" <martin@netbsd.org>
To: gnats-bugs@gnats.NetBSD.org
Cc: 
Subject: PR/50698 CVS commit: [netbsd-7] src/lib/libm/src
Date: Sun, 6 Mar 2016 18:00:08 +0000

 Module Name:	src
 Committed By:	martin
 Date:		Sun Mar  6 18:00:08 UTC 2016

 Modified Files:
 	src/lib/libm/src [netbsd-7]: e_hypotf.c

 Log Message:
 Pull up following revision(s) (requested by gson in ticket #1120):
 	lib/libm/src/e_hypotf.c: revision 1.10
 Fix incorrect magic numbers in scaling.  From FreeBSD commit 23397, by
 way of Rin Okuyama.  Fixes PR lib/50698.


 To generate a diff of this commit:
 cvs rdiff -u -r1.9 -r1.9.40.1 src/lib/libm/src/e_hypotf.c

 Please note that diffs are not public domain; they are subject to the
 copyright notices on the relevant files.

From: Andreas Gustafsson <gson@gson.org>
To: Rin Okuyama <okuyama@flex.phys.tohoku.ac.jp>
Cc: gnats-bugs@NetBSD.org
Subject: Re: lib/50698: hypotf() of small numbers yields infinity
Date: Fri, 11 Mar 2016 17:14:15 +0200

 On Feb 1, Rin Okuyama wrote:
 > On 2016/01/26 0:01, Andreas Gustafsson wrote:
 > >>   Other changes have been made for libm in FreeBSD. IMO, we should merge
 > >>   reasonable ones.
 > >
 > > Additional patches for the ones you think are reasonable are welcome.
 > > More test cases would be good, too.
 > 
 > Sorry for the late reply.

 Likewise :)

 > I examined FreeBSD's libm, and found it much
 > harder than I had imagined to merge their changes into ours.
 > 
 > There are three kinds of differences: (1) trivial bug fixes as in this
 > case, (2) long double functions missing in our libm, and (3) non-trivial
 > changes to improve the accuracy or processing speed.
 > 
 > We can readily merge (1) by patiently comparing the two implementations.
 > 
 > For (2), we have not fully implemented long double functions. Instead,
 > we have wrappers to double functions in src/lib/libm/src/ldbl_dummy.c.
 > Needless to say, this is problematic for architectures with
 > LDBL_MANT_DIG != DBL_MANT_DIG. They have long double functions both for
 > LDBL_MANT_DIG == 64 (amd64, i386, ia64, and m68k), and 113 (aarch64,
 > mips64, and sparc64). If we can merge it, we then have long double
 > functions for all architectures except powerpc64 (LDBL_MANT_DIG == 106).
 > 
 > However, for (2) and (3), there remains a big problem; how can we test
 > them? Addition to ATF-based tests ported from us, they have their own
 > test cases:
 > 
 >    https://svnweb.freebsd.org/base/head/lib/msun/tests/
 > 
 > Function returns for some special arguments are checked in their test.
 > But is it enough? Who can check hundreds of magic numbers?

 These are all good questions, but I'm afraid I don't have the answers.
 In any case, the specific bug reported in this PR has now been fixed
 and pulled up to -7, so I will close the PR.  If you have suggestions
 for additional tests or fixes, either for hypotf() or for libm in
 general, feel free to submit them as separate PRs.
 -- 
 Andreas Gustafsson, gson@gson.org

State-Changed-From-To: pending-pullups->closed
State-Changed-By: gson@NetBSD.org
State-Changed-When: Fri, 11 Mar 2016 15:16:39 +0000
State-Changed-Why:
Pullup done.


>Unformatted:

Home	PR Database Search