Re: Confused with frexp and atangent (mathfp branch)

[Jeff Johnston <jjohnstn at redhat dot com>]

Roman Belenov wrote:
> I'm working with a toolchain that uses newlib with mathfp branch used for
> libm. I've encountered some problems - it seems that newlib functions give
> incorrect results.
>
> First, frexp(0., &exp) writes -1022 to exp (instead of 0 required by standard)
> as a result of calculating
>     *exp = ((hd & 0x7ff00000) >> 20) - 1022;
> (since 0. is represented as all-zero-bits-value). Probably explicit test for
> zero (together with NAN and INF) is required.

Yes, a fix is required.  The check for special values is done on the 
result instead of the input.

> Second, there is confusing usage of branch variable in atangent() function,
> resulting in incorrect results from atan2. As a result of the following code
>
>   if (arctan2)
>     {
>       if (u < 0.0 || branch == 2)
>         res = __PI - res;
>       if (v < 0.0 || branch == 1)
>         res = -res;
>     }
>
> result is always set to PI for very small values of y/x (including y==0. due
> to aforementioned behaviour of frexp) and to -PI/2 for very high numbers
> ("very" means that preceding exponent test sets branch to nonzero value),
> since signs of u and v are effectively ignored. Seems like checks for branch
> here are erroneous. I don't have the textbook mention in s_atangent.c
> available, so can't check whether this snippet comes from its text and whether
> there is any rationale.
>
> Can anybody elaborate on these issues ?

Another bug.  Yes, it appears the check for branch is extraneous and 
causing the problem.  Please try the accompanying patch.

-- Jeff J.



? mathfp.patch
Index: s_atangent.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/mathfp/s_atangent.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_atangent.c
--- s_atangent.c	20 Oct 2003 18:46:38 -0000	1.2
+++ s_atangent.c	9 Jan 2006 22:29:18 -0000
@@ -197,9 +197,9 @@ _DEFUN (atangent, (double, double, doubl
 
   if (arctan2)
     {
-      if (u < 0.0 || branch == 2)
+      if (u < 0.0)
         res = __PI - res;
-      if (v < 0.0 || branch == 1)
+      if (v < 0.0)
         res = -res;
     }
   else if (x < 0.0)
Index: s_frexp.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/mathfp/s_frexp.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_frexp.c
--- s_frexp.c	20 Oct 2003 18:46:38 -0000	1.2
+++ s_frexp.c	9 Jan 2006 22:29:18 -0000
@@ -82,6 +82,17 @@ double frexp (double d, int *exp)
   double f;
   __uint32_t hd, ld, hf, lf;
 
+  /* Check for special values. */
+  switch (numtest (d))
+    {
+      case NAN:
+      case INF:
+        errno = EDOM;
+      case 0:
+        *exp = 0;
+        return (d);
+    }
+
   EXTRACT_WORDS (hd, ld, d);
 
   /* Get the exponent. */
@@ -94,16 +105,6 @@ double frexp (double d, int *exp)
 
   INSERT_WORDS (f, hf, lf);
 
-  /* Check for special values. */
-  switch (numtest (f))
-    {
-      case NAN:
-      case INF:
-        errno = EDOM;
-        *exp = 0;
-        return (f);
-    }
-
   return (f);
 }
 
Index: sf_frexp.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/mathfp/sf_frexp.c,v
retrieving revision 1.1.1.1
diff -u -p -r1.1.1.1 sf_frexp.c
--- sf_frexp.c	17 Feb 2000 19:39:52 -0000	1.1.1.1
+++ sf_frexp.c	9 Jan 2006 22:29:18 -0000
@@ -24,6 +24,17 @@ float frexpf (float d, int *exp)
   float f;
   __int32_t wf, wd;
 
+  /* Check for special values. */
+  switch (numtestf (d))
+    {
+      case NAN:
+      case INF:
+        errno = EDOM;
+      case 0:
+        *exp = 0;
+        return (d);
+    }
+
   GET_FLOAT_WORD (wd, d);
 
   /* Get the exponent. */
@@ -35,16 +46,6 @@ float frexpf (float d, int *exp)
 
   SET_FLOAT_WORD (f, wf);
 
-  /* Check for special values. */
-  switch (numtestf (f))
-    {
-      case NAN:
-      case INF:
-        errno = EDOM;
-        *exp = 0;
-        return (f);
-    }
-
   return (f);
 }