Use of C99 extra long double math functions after r236148

[Bruce Evans]

On Sun, 22 Jul 2012, Stephen Montgomery-Smith wrote:

> On 07/22/2012 11:12 PM, Bruce Evans wrote:
>> On Sun, 22 Jul 2012, Stephen Montgomery-Smith wrote:
>> 
>>> But I will say that your latest version of clog doesn't do as well as
>>> mine with this input:
>>> 
>>>    x = unur_sample_cont(gen);
>>>    y = unur_sample_cont(gen);
>>>    h = hypot(x,y);
>>>    x = x/h;
>>>    y = y/h;
>>> 
>>> I was able to get ULPs less than 2.  Your program gets ULPs more like
>>> up to 4000.
>> ...
>> What are the actual x and y?  I'm not set up to use mpfr.
>
> The code segment didn't I showed didn't use mpfr.  It used unuran. Basically 
> I am generating random numbers uniformly distributed on the disk |z|=1.  You 
> could also do it using
> x = cos(t)
> y = sin(t)
> where t is a random real number in the interval [0,2 pi].

I see.  You are generating |z| near 1 as a side effect of the inaccuracies
in sin^2 + cos^2 == 1 :-).

I was thinking of generating such numbers using sqrt().  sqrt() is
perfectly rounded, so this gives more control.  Start with x in
[sqrt(2), 1).  (I would first try stepping (uniformly) through all of
float space).  Let y = sqrt(1 - x^2).  Then x^2 + y^2 should be 1, but
due to inaccuracies it won't be.  Since sqrt(2) is irrational and
sqrt() is perfectly rounded, in round-to-nearest mode y should always
be too large or too small by less than half an ulp.  x^2 is imprecise
too, so the error for y may be larger.  Check a few values on both
sides to find the correctly rounded value y.  Then this y or the
`nextafter' value on one side if it makes x^2 + y^2 (in infinite
precision) as close as possible to 1 as possible for this x (by
monotonicity).

I expect the differences to have a distribution that is somewhat unifom
in the mantissa bits, so that the bad cases turn up fast with almost
any distribution of x, but worst cases never turn up with blind
nonexhaustive testing (testing twice as many cases tends to double the
worst case error measured in ulps).

There is no obvious reason why the worst cases for the size of
|x^2 + y^2 - 1| should only be near x = 1.  However, as x moves away
from 1 towards sqrt(2), it becomes easier to calculate x^2 + y^2
accurately enough using only doubled precision.

> I have been putting a lot of work into casinh, casin, cacosh and cacos, 
> getting the branches correct.  That has exhausted me.

I got exhausted with just clog :-).

Bruce