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