C11 conformance of casinl-like functions.

Montgomery-Smith, Stephen 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
>>>> finished,
>>>> 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.
>> http://faculty.missouri.edu/~stephen/software/#catrig
> 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.
> Bruce

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)
as (approximately)
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 mailing list