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