In libm, sin(qNaN) doesn't expect FE_INVALID ?

In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Ruinland ChuanTzu Tsai]

Hi all,
sorry for the bothering.

Recently, as I'm testing some modification on libm, I happen to realize
the fact that glibc's testsutie doesn't expect sin(+-qNaN) to trigger
FE_INVALID, which is designed in `math/libm-test-sin.inc` :

TEST_f_f (sin, qnan_value, qnan_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED)
TEST_f_f (sin, -qnan_value, qnan_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED)

Yet I think that if the program issued sin(qNaN), it should be okay to
be given an invalid operation excpetion. 

If that's the case, then appending `INVALID_EXCEPTION_OK` to the
expected exception list of sin(qNaN) and sin(-qNaN) should be benign.

Though I'm neither an expert on libm nor IEEE standards, I'm wondering
will there be any concern of such behavior ( raising FE_INVALID on
sin(+-qNaN) ) ?

Sincerely,
Ruinland

Re: In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Adhemerval Zanella]

On 03/09/2020 09:34, Ruinland ChuanTzu Tsai wrote:
> Hi all,
> sorry for the bothering.
> 
> Recently, as I'm testing some modification on libm, I happen to realize
> the fact that glibc's testsutie doesn't expect sin(+-qNaN) to trigger
> FE_INVALID, which is designed in `math/libm-test-sin.inc` :

AFAIK this is the expected behavior for *quiet* NaN for most symbols, it should
not signal invalid exception where applicable.

> 
> TEST_f_f (sin, qnan_value, qnan_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED)
> TEST_f_f (sin, -qnan_value, qnan_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED)
> 
> Yet I think that if the program issued sin(qNaN), it should be okay to
> be given an invalid operation excpetion. 
> 
> If that's the case, then appending `INVALID_EXCEPTION_OK` to the
> expected exception list of sin(qNaN) and sin(-qNaN) should be benign.

I think we should honor the quiet nan specification where applicable. Why
exactly are you trying to do that is requiring return a signaling NaN
for quiet NaN input?

> 
> Though I'm neither an expert on libm nor IEEE standards, I'm wondering
> will there be any concern of such behavior ( raising FE_INVALID on
> sin(+-qNaN) ) ?
> 
> Sincerely,
> Ruinland
> 

Re: In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Joseph Myers]

On Thu, 3 Sep 2020, Ruinland ChuanTzu Tsai wrote:

> Hi all,
> sorry for the bothering.
> 
> Recently, as I'm testing some modification on libm, I happen to realize
> the fact that glibc's testsutie doesn't expect sin(+-qNaN) to trigger
> FE_INVALID, which is designed in `math/libm-test-sin.inc` :

See subclause 6.2 of IEEE 754: "Every general-computational and 
quiet-computational operation involving one or more input NaNs, none of 
them signaling, shall signal no exception, except fusedMultiplyAdd might 
signal the invalid operation exception (see 7.2).".

The special rule about fma(0, Inf, qNaN) (where the standard leaves 
unspecified whether the NaN operand takes precedence over the 
multiplication 0 * Inf that would raise "invalid" if not fused with an 
addition) is handled through INVALID_EXCEPTION_OK in libm-test-fma.inc.  
The cases of INVALID_EXCEPTION_OK for various complex functions are where 
Annex G in the C standard mentions an optional "invalid" exception, and 
those for the significand function are empirical, reflecting that that 
function is not part of any standard.  No function from any standard 
should be using INVALID_EXCEPTION_OK in the testsuite without explicit 
standard wording about such an exception being permitted but optional.

Note that conversions from floating-point formats to *integer* formats 
signal "invalid" for quiet NaN inputs, since the NaN can't be represented 
in the output format.  Likewise, some but not all comparison operations 
signal "invalid" for quiet NaN arguments (see subclause 5.8 of IEEE 754).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Ruinland ChuanTzu Tsai]

Hi Dr. Myers, Zanella and all,
thanks very, very, very much for your precious comment and explanation.
It's a wonderful piece for I to read and learn from.
I was caught up in the middle of something and I'm so sorry for this
late reply.

Currrently I'm trying to improve several trigonal functions inside libm
and conducting some expriments to check against testsuite.

If I may, I have another question upon sin,atan(±0x4p-1076), the test-
suite expects FE_UNDERFLOW being raised.

In `math/libm-test-{sin,atan}.inc` , I cannot see the corresponding
expectation being set.  
And I discovered that in __sin() [sysdeps/ieee754/dbl-64/s_sin.c], it
will deliberately trigger the exception - - 

```
if (k < 0x3e500000)       /* if x->0 =>sin(x)=x */
  {
    math_check_force_underflow (x);
    retval = x;
  }
```

I've read IEEE Standard Section 7.5 which regulates that if a tiny non-
zero number is detected, an exception shall be raised. However I'm not
so sure about the reason why the magic number is set to `0x3e500000`.

Could you kindly point me to some directions about this magic number ?

Really appreciate the work and the comments from the community,
Ruinland ChuanTzu Tsai

Re: In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Joseph Myers]

On Tue, 8 Sep 2020, Ruinland ChuanTzu Tsai wrote:

> If I may, I have another question upon sin,atan(±0x4p-1076), the test-
> suite expects FE_UNDERFLOW being raised.
> 
> In `math/libm-test-{sin,atan}.inc` , I cannot see the corresponding
> expectation being set.  

That's because tests with finite inputs and finite mathematical results 
(possibly overflowing the floating-point type) generally go in 
auto-libm-test-in to have the expected results generated automatically by 
gen-auto-libm-tests; libm-test-*.inc are mainly for cases where the inputs 
or results involve exact infinities or NaNs.

> if (k < 0x3e500000)       /* if x->0 =>sin(x)=x */
>   {
>     math_check_force_underflow (x);
>     retval = x;
>   }
> ```
> 
> I've read IEEE Standard Section 7.5 which regulates that if a tiny non-
> zero number is detected, an exception shall be raised. However I'm not
> so sure about the reason why the magic number is set to `0x3e500000`.

That is a check for whether the argument is small enough that the argument 
itself is an accurate return value (then math_check_force_underflow does 
the separate check for whether it's small enough that the underflow 
exception must be raised).  As sin(x) = x - x^3/6 + ..., if x^2/6 is 
around 2^-54 or smaller, the error from returning x as the result is 
small.  That reasoning leads to the comparison of the leading bits of the 
representation against 0x3e500000 (i.e. comparing |x| against 0x1p-26).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Ruinland ChuanTzu Tsai]

Hi Dr. Myers,
thanks again for your detailed letter :-)

If you could bear with a little more, I have a final quetsion and would
like to have someone to discuss with - -

I'm a little bit confused by the implementation of ULPDIFF() inside
`math/libm-test-support.c` which is :

```
#define ULPDIFF(given, expected) \                                                                                                                            (FUNC(fabs) ((given) - (expected)) / ulp (expected)
```

and it looks _not_ really the same as the formula inside glibc's docu-
mentation [1] :
` |d.d...d - (z / 2^e)| / 2^(p - 1) `

( For a number z with the representation d.d…d·2^e and p is the number
of bits in the mantissa of the floating-point number representation. )

The denominator part of these two seems to have different meaning ?

Besides this issue, I would like to know that is there any written
policy for loosening or tightening the ULPs for mathematic functions ?
For instance, the libm-test-ulps for i386 has varied several times.
Yet I cannot find concrete discussions about making the precission
check either more strict or tolerant.

And if someone is introducing a new platform to glibc, are there any
rules to regulate ? e.g. "ccosh" mustn't have a ulp more than ......

Really appreciate your kindness,
Ruinland ChuanTzu Tsai

[1] https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html

Re: In libm, sin(qNaN) doesn't expect FE_INVALID ?

[Joseph Myers]

On Thu, 10 Sep 2020, Ruinland ChuanTzu Tsai wrote:

> I'm a little bit confused by the implementation of ULPDIFF() inside
> `math/libm-test-support.c` which is :
> 
> ```
> #define ULPDIFF(given, expected) \                                                                                                                            (FUNC(fabs) ((given) - (expected)) / ulp (expected)
> ```
> 
> and it looks _not_ really the same as the formula inside glibc's docu-
> mentation [1] :
> ` |d.d...d - (z / 2^e)| / 2^(p - 1) `
> 
> ( For a number z with the representation d.d…d·2^e and p is the number
> of bits in the mantissa of the floating-point number representation. )
> 
> The denominator part of these two seems to have different meaning ?

It looks like that formula from the manual should actually be multiplying 
by 2^(p-1), not dividing, to get an actual figure in ulps.

Note that there are at least two different measures of errors in ulps.  
The one used in glibc is that we take an ideal correctly rounded result, 
take the absolute value of the difference between that and the result 
returned by the function, and divide that by a unit in the last place of 
the correctly rounded result.  This gives an error that is almost always 
an integer number of ulps (it can be a non-integer if the result returned 
has a lower exponent than the correctly rounded result).  A correctly 
rounded result has a 0 ulps error by this definition (but that's not 
sufficient for being correctly rounded; correct rounding also requires the 
correct sign of 0 and correct exceptions).

Another version sometimes seen in the literature defines ulps not for a 
correctly rounded result but for the infinite-precision mathematical 
result.  When that's given as "the absolute value of the difference 
between the two floating-point numbers closest to x, one of which may 
equal x", note that if x is a power of 2 (and exceeds the magnitude of the 
least normal value), or rounds away from zero to a power of 2, then this 
gives a definition of ulp that's half the one used by glibc (and thus an 
error that's twice that of the glibc definition).  Then the error in a 
function is determined by comparing the rounded value to the 
infinite-precision value, in terms of ulps of the infinite-precision 
value.  With this definition, a correctly rounded result has error at most 
0.5 ulps in round-to-nearest mode and less than 1 ulp in other modes (but 
again, that's not sufficient for being correctly rounded).

> Besides this issue, I would like to know that is there any written
> policy for loosening or tightening the ULPs for mathematic functions ?

Only the functions bound to IEEE operations (sqrt, fma, etc.) are expected 
to be correctly rounded.  For others, people have typically found 
performance can be improved without introducing large errors.

My guess is that most functions could be made to achieve 1ulp errors in 
round-to-nearest and 2ulp in other modes (whichever definition is used) 
without making performance worse.

> And if someone is introducing a new platform to glibc, are there any
> rules to regulate ? e.g. "ccosh" mustn't have a ulp more than ......

The general rule for new platforms is to avoid having 
architecture-specific function implementations that aren't actually 
needed, and to improve performance by improving the generic C 
implementations instead; see <https://sourceware.org/glibc/wiki/NewPorts>.  
Architecture-specific versions of functions such as fma that are fully 
bound to IEEE operations may make sense, where there are relevant hardware 
instructions.  Architecture-specific versions of transcendental functions 
are almost surely a bad idea.  Once you're using the 
architecture-independent implementations, you should have the same ulps as 
for most other platforms (modulo minor differences arising from compiler 
choices in whether to contract operations, if one of the platforms has 
fused multiply-add instructions).

The only architecture-specific implementation of ccosh is for m68k (the 
alpha version is just dealing with compatibility for past ABI changes).  
The m68k version really ought to go away because of its use of fsincos 
(see bug 13742 regarding use of fsincos on m68k, and note that emulators 
may well not accurately reflect hardware inaccuracy there).

-- 
Joseph S. Myers
joseph@codesourcery.com