Confused with frexp and atangent (mathfp branch)

Jeff Johnston jjohnstn@redhat.com
Tue Jan 10 16:53:00 GMT 2006


Great.  Patch checked in.

-- Jeff J.

Roman Belenov wrote:
> Thanks, the patch fixes the problem.
> 
> Jeff Johnston <jjohnstn@redhat.com> writes:
> 
> 
>>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);
>> }
>> 
> 
> 




More information about the Newlib mailing list