[Rd] pcauchy precision (PR#6756)

[maechler at stat.math.ethz.ch]

>>>>> "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