[Rd] pcauchy precision (PR#6756)
maechler at stat.math.ethz.ch
maechler at stat.math.ethz.ch
Wed Apr 14 15:15:40 CEST 2004
>>>>> "Morten" == Morten Welinder <terra at gnome.org>
>>>>> on Tue, 13 Apr 2004 16:15:15 +0200 (CEST) writes:
Morten> The current code has xbig such that the x^5 term is
Morten> beyond the precision limit (although I would use 2.5
Morten> instead of 5 in the xbig formula). You get xbig
Morten> ~5478, I get xbig ~179513609. That is about log_2
Morten> (5478) = 12.5 bits lost to cancelation in your case
Morten> and 27.5 bits in my case. If you want those numbers
Morten> down and still use a series expansion, you need to
Morten> have a lot more terms so you can get xbig down.
Morten> By the way, your xbig point is valid only for
Morten> log_p==FALSE. For log_p==TRUE it needs to be higher.
ok, that's an important point.
Morten> I don't see why you would be worried over atan's
Morten> performance in the [-1;+1] interval given that the
Morten> main branch already relies on atan in that interval.
I wasn't worried really; just I assume log() and log1p() to be
very accurate (only loose 1 or 2 bits) where (bad) system's atan()'s
may loose more.
Morten> Numbers... There is nothing like numbers.
Indeed,
that was the evidence I've been asking for.
Morten> I hacked up a comparison in IEEE double space,
Morten> i.e., 53-bit mantissa. I calculated a reference
Morten> value using 113-bit mantissa arithmetic (using the
Morten> 1/x method
(which has an obvious "design bias" though -- but I agree it's not
relevant here)
Morten> and Sun's high- performance libraries)
Morten> with the result cast back to 53-bit.
Morten> Fix loc=0 and scale=1 and let eR be the current R
Morten> method's relative error and let eN be the proposed
Morten> method's relative error. Then we have the table
Morten> below.
Morten> Conclusions...
Morten> 1. The series expansion loses 3-4 digits just below xbig.
Morten> 2. The series expansion loses 2-3 digits just above xbig for log_p==TRUE.
Morten> 3. The series version isn't even used for x>xbig and
Morten> lower==TRUE. That makes R's log_p==TRUE performance awful.
"awful" is strong; but otherwise I agree.
Together with the elegance of your proposition, I now agree with
you.
Thank you, Morten!
Martin
More information about the R-devel
mailing list