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