C11 conformance of casinl-like functions.
stephen at missouri.edu
Wed Feb 8 17:28:59 UTC 2017
On 02/08/17 09:10, Bruce Evans wrote:
> On Wed, 8 Feb 2017, Montgomery-Smith, Stephen wrote:
>> On 02/08/17 06:47, mokhi wrote:
>>>> I think you mean acosl, asinl, ...
>>> yeah :D
>>>> These were implemented quite well in 2012-2013, but not quite
>>>> and not committed. Only the float and double version were committed.
>>>> The raw versions are still available in
>>>> https://people.freebsd.org/~stephen/catrig*.c. These have rotted and
>>>> require some editing. Compare with the committed parts to see most
>>>> of the necessary editing.
>>> Okay, I'll do.
>>> Would you like that I add you on Phabricator for reviewing when my
>>> editing is done?
>>> Thanks and best wishes, Mokhi.
>> I would love to see these functions implemented. I spent a lot of time
>> on them, and they went through a lot of testing to make sure they were
>> accurate, even in very extreme edge cases.
> It was you that implemented them :-). Except for many onerous details
> like testing on all supported arches and committing them.
>> Here are more details about the implementation.
> It says that the the results are good to 4 ulps except in float precision.
> I dout that that is correct. Higher precisions are just less likely to
> get near the corner cases where the errors are larger.
> Corner cases are probably only near zeros and poles, but I test the
> real and imaginary parts separately and this gives a dubious size for
> an ulp (for the error relative to each part instead of relative to
> the absolute value of the result). This gives 1-dimensional sets of
> corner cases. The functions try to get a small error by this measure
> and mostly succeed.
I did some very extensive testing. In fact, I was able to find an edge
case in which the paper by Hull, Fairgrieve and Tang had a mistake if
one measured the real and imaginary parts separately for cacos/cacosh.
I corrected this in the Boost library:
The algorithm in the paper by Hull, Fairgrieve and Tang would compute
cacos(1.00000002785941 + I*5.72464869028403e-200)
0 - I*0.000236048349018331
whereas it should be (approximately)
2.42520172707401e-196 - I*0.000236048349018331.
(This is not only close to a zero of acos, I think it is also the end of
a branch point.)
More information about the freebsd-numerics