X-Authentication-Warning: delorie.com: mail set sender to djgpp-bounces using -f
X-Recipient: djgpp AT delorie DOT com
Message-ID: <5A74B159.4020603@gmx.de>
Date: Fri, 02 Feb 2018 19:43:37 +0100
From: "Juan Manuel Guerrero (juan DOT guerrero AT gmx DOT de) [via djgpp AT delorie DOT com]" <djgpp AT delorie DOT com>
User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; de; rv:1.9.2.13) Gecko/20101206 SUSE/3.1.7 Thunderbird/3.1.7
MIME-Version: 1.0
To: djgpp AT delorie DOT com
Subject: Re: Fixing various bugs in frexp.S
References: <5A4946EB DOT 3090500 AT gmx DOT de> <p2vmp0$14s4$1 AT gioia DOT aioe DOT org> <279ab41f-d898-4bcf-b5d4-ad9bf86e7a06 AT googlegroups DOT com> <83lgh86tza DOT fsf AT gnu DOT org>
In-Reply-To: <83lgh86tza.fsf@gnu.org>
Content-Type: text/plain; charset=ISO-8859-15; format=flowed
Content-Transfer-Encoding: 7bit
X-Provags-ID: V03:K0:2BKdukfhcIXJGkTwwHb6480QfdmIXXmFzaDU9Y0t5UxqtlPPhPn
 HGEb+IMs6nOpwSUuzyxVN2asog7W3AmVCCiFk9JYTgzOme7YbG8cfexukj3T/1xnxvv0Pia
 IfXKff7ERRnNccV14VNEK+syoCm82Bc3WAqWSjZt1W3rhEMBSMP6fKWTBzbBlxiokpAwvmO
 mY0vjx6YKo3c2yff/S6qQ==
X-UI-Out-Filterresults: notjunk:1;V01:K0:/7949yHaRn4=:DKGZVCyodborNr22fyY9xF
 Q+ZUpEeemMcTs3uP193ggWoyKZMSgIEqv35JbkGNgNms7/0SOMi/Nsudn6S5PZgFRWAus/o3O
 P2y9xRkMWmTxm+Td3Po2KhuQfminkeEQkja7nC+QBlylmw2mF5osfwQ8MjXaXYgtkV5d3ZkmK
 Osk3UNzVpwlkIiWLmIvcft1QwLPh4gni3nsB2aYAx5ll3q2kqLQZmf3CXBvpRIzAF14AMLR5z
 eIZX582cYzlo5mWkmtFkipr/HMR97HXRNW7mc0RvuNMLQvhhnGEzRzx+lTPEyc1o/87X4aMmU
 cpo3+okzAZzEu2XB0yXhhiVDsPVHIUPoMvQm67a8vKh5VUZ3nKYkjAFNiI2RmHIdCR5SNb82e
 W4rNKsc0LaCmq+1s4n+pPcXTovCLnZLetdJR6H7JNsmQ1TJdrhCKlD0SUmvsTQzgJylOUlVXM
 xRLzy/ksSip2N0c7SzuwYmYYekymD/UrOPWXBuXH/c+O79/Q4VpGzzUGQZLMaPgU3jxUPdyZq
 ia+kZLRXPQKc2hMwtEvi3fhfflqtALHLx9Eu8+iYNxgid+hmhVPwInaovEo7YkNOcz3ICnMHk
 FuBvG/pCtbmTinmVmQJwC9R42lIUxYuIGjEshAKBREZib4fz8WDQTw75DrBQSHCGUTZdm/I/m
 RVioBc4dHjhtXNVTFFQsMB+VTGICaUASB2Ngz0qy9nzCzpcbOdg3t2y1ikEdHXj6PbmvLxlO9
 U8dn1r42RiUt1KiMY3Dq7YFMbzEe7vHenv2C3WQmM1Gw/Ub5FA9dYufimiVH8BMaAeirByTQx
 1DLG804NmKJfNQNEWhKwxqPsXuPXs25usJKGMFW95K8+7MLOjY=
Reply-To: djgpp AT delorie DOT com

Am 08.01.2018 20:17, schrieb Eli Zaretskii (eliz AT gnu DOT org) [via djgpp AT delorie DOT com]:
>> Date: Mon, 8 Jan 2018 10:27:20 -0800 (PST)
>> From: "Juan Manuel Guerrero (juan DOT guerrero AT gmx DOT de) [via djgpp AT delorie DOT com]"<djgpp AT delorie DOT com>
>>
>>>>    - for NaN and infinity arguments no sign checks are done thus the
>>>>      results are always positive
>>>
>>>> All these issues are serious deviations from the ansi/posix behavior.
>>>
>>> An NaN is a NaN. It has no positive or negative value. It has no value at
>>> all. I seriously doubt ansi/posix says anything about signedness of NaNs.
>>
>> The goal is not to generate a negative NaN but to preserve the sign of the
>> imput argument.  The implementation should no convert a -NaN in a +NaN IMO.
>> This way the implementation bahaves as the implementation on my linux box.
>
> I agree with Martin that the sign of a NaN is meaningless: you cannot
> detect that sign by comparing a NaN with zero, for example.
>
> Does glibc give any explanation (in its sources, for example) for
> this behavior?

I have checked that on a BSD system, frexp[f|l] always return
a positive NaN no matter what the sign of the input NaN is.
Thus I have adjusted the frexp[f|l] functions accordingly.
They will always return a positive sign for NaN.
This means that apart from libm/math/s_frexpl.c also libm/math/s_frexp.c
and libm/math/sf_frexp.c must be adjusted because these
functions have always preserved the sign of zero, infinity and
NaN.  An inspection of the linux code seems to indicate that
they have also based their code on the same old SunPro math
code that has also been used by DJGPP.  This means that DJGPP's
frexp[f] from libm  have always preserved the sign of NaN.

libc/ansi/math/frexp.S has simple bugs. It only returns correct
values for normalized numbers and for positive zero.  Also in
case of error it sets errno to EDOM and this is wrong as can be
read in any man page of frexp.  Case of error means that the
imput value is ether infinity or NaN.  For negative zero, +/- NaN
and +/- infinity the values are simple wrong.  Below are the results
as are generated by djdev205:

IN: -0   OUT: mantissa = -0  exp = -2147483647  err = 0
IN: 0   OUT: mantissa = 0  exp = 0  err = 0
IN: -nan   OUT: mantissa = -nan  exp = 0  err = 1
IN: nan   OUT: mantissa = -nan  exp = 0  err = 1
IN: -inf   OUT: mantissa = -nan  exp = 0  err = 1
IN: inf   OUT: mantissa = -nan  exp = 0  err = 1

In other words all negative values are wrong, Inf is returned
as NaN and the negative zero has -2147483647 as exponent.
It is surprising that no one has noted this in 20 years.

If no one has really serious objections, I will commit the
patch below in a couple of days.


Regards,
Juan M. Guerrero





Index: djgpp/src/libc/ansi/math/frexp.S
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libc/ansi/math/frexp.S,v
retrieving revision 1.5
diff -U 5 -r1.5 frexp.S
--- djgpp/src/libc/ansi/math/frexp.S	3 Jan 2018 15:52:09 -0000	1.5
+++ djgpp/src/libc/ansi/math/frexp.S	2 Feb 2018 18:26:13 -0000
@@ -2,12 +2,10 @@
  /* Copyright (C) 2017 DJ Delorie, see COPYING.DJ for details */
  /* Copyright (C) 2011 DJ Delorie, see COPYING.DJ for details */
  /* Copyright (C) 1999 DJ Delorie, see COPYING.DJ for details */
  pos_NaN:
  	.long	0x7FC00000
-neg_NaN:
-	.long	0xFFC00000

  pos_inf:
  	.long	0x7F800000
  neg_inf:
  	.long	0xFF800000
@@ -90,17 +88,13 @@
  	flds	pos_NaN
  	ret

  negarg:					/* arg is - inf or - NaN */
  	cmpl	$0x80000000, %eax
-	jnz	negNaN
+	jnz	posNaN

  	movl	4(%esp), %eax
  	testl	%eax, %eax
-	jnz	negNaN
+	jnz	posNaN

  	flds	neg_inf
  	ret
-
-negNaN:
-	flds	neg_NaN
-	ret
Index: djgpp/src/libm/math/s_frexp.c
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/s_frexp.c,v
retrieving revision 1.2
diff -U 5 -r1.2 s_frexp.c
--- djgpp/src/libm/math/s_frexp.c	2 Jan 2018 22:02:08 -0000	1.2
+++ djgpp/src/libm/math/s_frexp.c	2 Feb 2018 18:26:13 -0000
@@ -106,12 +106,17 @@
  {
  	__int32_t hx, ix, lx;
  	EXTRACT_WORDS(hx,lx,x);
  	ix = 0x7fffffff&hx;
  	*eptr = 0;
-	if(ix>=0x7ff00000||((ix|lx)==0)) return x;	/* 0,inf,nan */
-	if (ix<0x00100000) {		/* subnormal */
+	if(((ix&0x000fffff)|lx)==0) return x;	/* 0,inf */
+	if(ix>=0x7ff00000) {			/* nan */
+	    hx = hx&0x7fffffff;
+	    SET_HIGH_WORD(x,hx);
+	    return x;
+	}
+	if (ix<0x00100000) {			/* subnormal */
  	    x *= two54;
  	    GET_HIGH_WORD(hx,x);
  	    ix = hx&0x7fffffff;
  	    *eptr = -54;
  	}
Index: djgpp/src/libm/math/s_frexpl.c
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/s_frexpl.c,v
retrieving revision 1.1
diff -U 5 -r1.1 s_frexpl.c
--- djgpp/src/libm/math/s_frexpl.c	2 Jan 2018 22:05:37 -0000	1.1
+++ djgpp/src/libm/math/s_frexpl.c	2 Feb 2018 18:26:13 -0000
@@ -24,14 +24,14 @@
   *	arg = x*2^exp.
   * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
   * with *exp=0.
   */

-#define IS_INF_OR_NAN(v)       ((v).ldt.exponent == 0x7FFF)
-#define IS_DENORMAL(v)         ((v).ldt.exponent == 0x0000 && (((v).ldt.mantissah & 0x80000000) == 0x00000000))
-#define IS_PSEUDO_DENORMAL(v)  ((v).ldt.exponent == 0x0000 && (((v).ldt.mantissah & 0x80000000) == 0x80000000))
-#define IS_ZERO(v)             (((v).ldt.mantissah | (v).ldt.mantissal) == 0x00000000)
+#define IS_DENORMAL(v)         (((v).ldt.exponent == 0x0000) && (((v).ldt.mantissah & 0x80000000) == 0x00000000))
+#define IS_PSEUDO_DENORMAL(v)  (((v).ldt.exponent == 0x0000) && (((v).ldt.mantissah & 0x80000000) == 0x80000000))
+#define IS_INF_OR_ZERO(v)      ((((v).ldt.mantissah & 0x7FFFFFFF) | (v).ldt.mantissal) == 0x00000000)
+#define IS_NAN(v)              (((v).ldt.exponent == 0x7FFF) && !IS_INF_OR_ZERO(v))

  #include <libc/ieee.h>
  #include "fdlibm.h"

  #ifdef __STDC__
@@ -49,11 +49,15 @@
  #endif
  {
          _longdouble_union_t value;
          value.ld = x;
  	*eptr = 0;
-	if (IS_INF_OR_NAN(value) || IS_ZERO(value)) return x;	/* 0,inf,nan */
+	if (IS_INF_OR_ZERO(value)) return x;	/* 0 or inf */
+	if (IS_NAN(value)) {			/* nan */
+	    value.ldt.sign = 0;
+	    return value.ld;
+        }
  	if (IS_DENORMAL(value)) {		/* subnormal */
  	    value.ld *= two65;
  	    *eptr = -65;
  	}
  	else if (IS_PSEUDO_DENORMAL(value)) {	/* pseudo subnormal */
Index: djgpp/src/libm/math/sf_frexp.c
===================================================================
RCS file: /cvs/djgpp/djgpp/src/libm/math/sf_frexp.c,v
retrieving revision 1.1
diff -U 5 -r1.1 sf_frexp.c
--- djgpp/src/libm/math/sf_frexp.c	7 Feb 1998 14:13:20 -0000	1.1
+++ djgpp/src/libm/math/sf_frexp.c	2 Feb 2018 18:26:13 -0000
@@ -31,11 +31,16 @@
  {
  	__int32_t hx, ix;
  	GET_FLOAT_WORD(hx,x);
  	ix = 0x7fffffff&hx;
  	*eptr = 0;
-	if(ix>=0x7f800000||(ix==0)) return x;	/* 0,inf,nan */
+	if((ix&0x007fffff)==0) return x;	/* 0,inf */
+	if(ix>=0x7f800000) {			/* nan */
+	    hx = hx&0x7fffffff;
+	    SET_HIGH_WORD(x,hx);
+	    return x;
+	}
  	if (ix<0x00800000) {		/* subnormal */
  	    x *= two25;
  	    GET_FLOAT_WORD(hx,x);
  	    ix = hx&0x7fffffff;
  	    *eptr = -25;
Index: djgpp/tests/cygnus/frexpl-t.c
===================================================================
RCS file: /cvs/djgpp/djgpp/tests/cygnus/frexpl-t.c,v
retrieving revision 1.2
diff -U 5 -r1.2 frexpl-t.c
--- djgpp/tests/cygnus/frexpl-t.c	27 Jan 2018 19:12:08 -0000	1.2
+++ djgpp/tests/cygnus/frexpl-t.c	2 Feb 2018 18:26:13 -0000
@@ -103,13 +103,13 @@
    {{.ldt = {0x0U, 0x80000000U, 0x7FFFU, 0}},  {.ldt = {0x0U, 0x80000000U, 0x7FFFU, 0}}, 0}, /* Inf */
    {{.ldt = {0x0U, 0x80000000U, 0x7FFFU, 1}},  {.ldt = {0x0U, 0x80000000U, 0x7FFFU, 1}}, 0}, /* -Inf */

    /* NaNs. */
    {{.ldt = {0x1U, 0x80000000U, 0x7FFFU, 0}},  {.ldt = {0x1U, 0x80000000U, 0x7FFFU, 0}}, 0}, /* SNaN */
-  {{.ldt = {0x1U, 0x80000000U, 0x7FFFU, 1}},  {.ldt = {0x1U, 0x80000000U, 0x7FFFU, 1}}, 0}, /* -SNaN */
+  {{.ldt = {0x1U, 0x80000000U, 0x7FFFU, 1}},  {.ldt = {0x1U, 0x80000000U, 0x7FFFU, 0}}, 0}, /* -SNaN */
    {{.ldt = {0x0U, 0xFFFFFFFFU, 0x7FFFU, 0}},  {.ldt = {0x0U, 0xFFFFFFFFU, 0x7FFFU, 0}}, 0}, /* QNaN */
-  {{.ldt = {0x0U, 0xFFFFFFFFU, 0x7FFFU, 1}},  {.ldt = {0x0U, 0xFFFFFFFFU, 0x7FFFU, 1}}, 0}, /* -QNaN */
+  {{.ldt = {0x0U, 0xFFFFFFFFU, 0x7FFFU, 1}},  {.ldt = {0x0U, 0xFFFFFFFFU, 0x7FFFU, 0}}, 0}, /* -QNaN */

    /* Number. */
    {{.ldt = {0x2168C000U, 0xC90FDAA2U, 0x3FFFU + 0x0001U, 0}},  {.ldt = {0x2168C000U, 0xC90FDAA2U, 0x3FFEU, 0}}, 2}, /* PI */
    {{.ldt = {0x2168C000U, 0xC90FDAA2U, 0x3FFFU + 0x0001U, 1}},  {.ldt = {0x2168C000U, 0xC90FDAA2U, 0x3FFEU, 1}}, 2}, /* -PI */

@@ -196,10 +196,11 @@

  static const size_t n_tests_long_double = sizeof(tests_long_double) / sizeof(tests_long_double[0]);

  #undef  IS_EQUAL
  #define IS_EQUAL(a, b) (((a).ldt.sign == (b).ldt.sign) && ((a).ldt.exponent == (b).ldt.exponent) && ((a).ldt.mantissah == (b).ldt.mantissah) && ((a).ldt.mantissal == (b).ldt.mantissal))
+#define IS_EQUAL_NAN(a, b) ((((a).ldt.sign == (b).ldt.sign) || ((a).ldt.sign != (b).ldt.sign))&& ((a).ldt.exponent == (b).ldt.exponent) && ((a).ldt.mantissah == (b).ldt.mantissah) && ((a).ldt.mantissal == (b).ldt.mantissal))


  int frexpl_test(void)
  {
    unsigned int i, counter;
@@ -210,17 +211,19 @@
      _longdouble_union_t mantissa;
      mantissa.ld = frexpl(tests_long_double[i].value.ld, &exponent);

      if (IS_EQUAL(tests_long_double[i].mantissa, mantissa) && tests_long_double[i].exponent == exponent)
        counter++;
+    else if(IS_EQUAL_NAN(tests_long_double[i].mantissa, mantissa) && tests_long_double[i].exponent == exponent)
+      counter++;
      else
        printf("frexpl test failed:  value to frexpl = %.12Lf\n"
               "mantissa result = %.12Lf  exponent result = %d\n"
               "mantissa should be = %.12Lf  exponent should be = %d\n",
               tests_long_double[i].value.ld, mantissa.ld, exponent, tests_long_double[i].mantissa.ld, tests_long_double[i].exponent);
    }

-  result = counter < n_tests_long_double || result ? 1 : 0;
+  result = counter < n_tests_long_double ? 1 : 0;
    printf("%s\n", result ? "frexpl test failed." : "frexpl test succeded.");

    return result;
  }