Use of C99 extra long double math functions after r236148
Stephen Montgomery-Smith
stephen at missouri.edu
Sun Aug 12 23:10:47 UTC 2012
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.
>
> I may have broken the double version when working mostly on the float
> version recently.
>
> 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].
> Since the float version gets errors of 4096 ulps (12 bits wrong), the
> double version is sure to get errors of [much more than] 12 + (53-24)
> = 41 bits wrong. That is 2 tera ulps. Not noticing such enormous
> errors indicates that the problematic cases haven't been tested. I
> think you are right that it needs more like tripled double precision
> -- with merely doubled double precision, it can probably get all 53
> mantissa bits and the sign bit wrong too (sign of (|z|^2 - 1)). That
> is total loss of precision (TLOSS), and should be handled by returning
> NaN. Sign errors are especially interesting with complex functions
> and even for real log() applied to a real function, since they may
> change the branch. I got TLOSS including sign errors in the loghypotf()
> result in intermediate version due to bugs in the doubling of float
> precision. Before the attempted doubling, TLOSS might have been the
> usual case for z near 1!
I will have another go at this code, maybe tomorrow, maybe later.
I have been putting a lot of work into casinh, casin, cacosh and cacos,
getting the branches correct. That has exhausted me.
>> I have to say that I consider a ULP of 4000 under these very extreme
>> circumstances to be acceptable. Definitely acceptable if the code
>> goes a whole lot faster than code that has a ULP of less than 2.
>
> "An ULP of 4000" is unusual terminology. An ulp is a unit, not a count.
I am so different than you guys. I have no problem with lack of
consistency in notation, no problem with language being abused, no
problem with people using different programming styles - certainly no
problems with ULP's on the large side (oops, errors with large numbers
of ULP's) , and ...
> I haven't figured out how to cut down the amount of mail generated by
> this thread. Sorry to add to it :-).
... no problems deleting emails that I don't want (but surely with all
the email volume we have, this does suggest we need a freebsd-numerics
mailing list, doesn't it?), and ...
>
>> On 07/22/2012 02:29 PM, Bruce Evans wrote:
>>> Replying again to this...
>
> Top posting is one way :-).
... most definitely no problem with top posting! I try not to top post
on the FreeBSD mailing list because I know so many people there are
bothered by it. But most of the people I communicate with (professors
in academic institutions, friends and relatives) top post all the time.
But I should have remembered not to top post with you guys.
I understand that you guys and I think differently. I don't think that
I am right and you guys are wrong. But conversely, I hope that you
don't think I am wrong and you are right! And I hope that you can
appreciate the qualities I bring to the discussion, just as I appreciate
the qualities that you bring.
Stephen
More information about the freebsd-numerics
mailing list