From terra at gnome.org  Tue Apr 13 16:15:15 2004
From: terra at gnome.org (terra@gnome.org)
Date: Tue Apr 13 16:15:18 2004
Subject: [Rd] pcauchy precision (PR#6756)
Message-ID: <20040413141515.37A1FFE9B@slim.kubism.ku.dk>


The current code has xbig such that the x^5 term is beyond the precision
limit (although I would use 2.5 instead of  5 in the xbig formula).  You
get xbig ~5478, I get xbig ~179513609.  That is about log_2 (5478) = 12.5
bits lost to cancelation in your case and 27.5 bits in my case.  If you
want those numbers down and still use a series expansion, you need to have
a lot more terms so you can get xbig down.

By the way, your xbig point is valid only for log_p==FALSE.  For log_p==TRUE
it needs to be higher.

I don't see why you would be worried over atan's performance in the [-1;+1]
interval given that the main branch already relies on atan in that interval.

Numbers...   There is nothing like numbers.  I hacked up a comparison in
IEEE double space, i.e., 53-bit mantissa.  I calculated a reference value
using 113-bit mantissa arithmetic (using the 1/x method and Sun's high-
performance libraries) with the result cast back to 53-bit.

Fix loc=0 and scale=1 and let eR be the current R method's relative error
and let eN be the proposed method's relative error.  Then we have the table
below.

Conclusions...

1. The series expansion loses 3-4 digits just below xbig.
2. The series expansion loses 2-3 digits just above xbig for log_p==TRUE.
3. The series version isn't even used for x>xbig and lower==TRUE.  That
   makes R's log_p==TRUE performance awful.

Morten


troll:~> for x in 1 10 100 1000 5478 5459 6000 7000 8000 9000 10000 100000 1000000 10000000 100000000; do ./a.out $x 0; done
x=         1  y=       -0.2876820725  lt=1  log=1  eR=                  -0  eN=                  -0
x=         1  y=                0.75  lt=1  log=0  eR=                   0  eN=                   0
x=         1  y=        -1.386294361  lt=0  log=1  eR=                  -0  eN=                  -0
x=         1  y=                0.25  lt=0  log=0  eR=                   0  eN=                   0
x=        10  y=      -0.03223967553  lt=1  log=1  eR=    -1.721827231e-15  eN=                  -0
x=        10  y=        0.9682744826  lt=1  log=0  eR=                   0  eN=                   0
x=        10  y=        -3.450633956  lt=0  log=1  eR=     3.860935836e-16  eN=                  -0
x=        10  y=       0.03172551743  lt=0  log=0  eR=    -1.531015449e-15  eN=     2.187164928e-16
x=       100  y=     -0.003188069262  lt=1  log=1  eR=    -2.122106203e-14  eN=                  -0
x=       100  y=        0.9968170072  lt=1  log=0  eR=     1.113768141e-16  eN=                   0
x=       100  y=        -5.749933404  lt=0  log=1  eR=     3.707222428e-15  eN=                  -0
x=       100  y=      0.003182992765  lt=0  log=0  eR=    -2.111865771e-14  eN=                   0
x=      1000  y=    -0.0003183604514  lt=1  log=1  eR=    -1.323068058e-13  eN=     1.702790293e-16
x=      1000  y=        0.9996816902  lt=1  log=0  eR=                   0  eN=                   0
x=      1000  y=        -8.052485498  lt=0  log=1  eR=     1.632420278e-14  eN=                  -0
x=      1000  y=     0.0003183097801  lt=0  log=0  eR=    -1.323278675e-13  eN=     1.703061358e-16
x=      5478  y=     -5.81086402e-05  lt=1  log=1  eR=    -1.604954363e-12  eN=     1.166137007e-16
x=      5478  y=         0.999941893  lt=1  log=0  eR=      1.11028754e-16  eN=                   0
x=      5478  y=        -9.753225247  lt=0  log=1  eR=     1.644635682e-13  eN=                  -0
x=      5478  y=     5.810695193e-05  lt=0  log=0  eR=    -1.605000994e-12  eN=     1.166170889e-16
x=      5459  y=    -5.831089269e-05  lt=1  log=1  eR=      -3.7128847e-13  eN=                  -0
x=      5459  y=        0.9999416908  lt=1  log=0  eR=                   0  eN=                   0
x=      5459  y=        -9.749750799  lt=0  log=1  eR=     3.807877628e-14  eN=                  -0
x=      5459  y=     5.830919264e-05  lt=0  log=0  eR=    -3.711830826e-13  eN=      1.16212612e-16
x=      6000  y=    -5.305305449e-05  lt=1  log=1  eR=     1.294887935e-12  eN=                  -0
x=      6000  y=        0.9999469484  lt=1  log=0  eR=    -1.110281927e-16  eN=                   0
x=      6000  y=        -9.844244643  lt=0  log=1  eR=                  -0  eN=                  -0
x=      6000  y=     5.305164721e-05  lt=0  log=0  eR=                   0  eN=                   0
x=      7000  y=     -4.54738745e-05  lt=1  log=1  eR=     9.137564969e-13  eN=      1.49014432e-16
x=      7000  y=        0.9999545272  lt=1  log=0  eR=                   0  eN=                   0
x=      7000  y=        -9.998395321  lt=0  log=1  eR=     1.776641933e-16  eN=                  -0
x=      7000  y=     4.547284057e-05  lt=0  log=0  eR=                   0  eN=     1.490178201e-16
x=      8000  y=    -3.978952716e-05  lt=1  log=1  eR=    -2.308452986e-12  eN=                  -0
x=      8000  y=        0.9999602113  lt=1  log=0  eR=     1.110267201e-16  eN=                   0
x=      8000  y=        -10.13192671  lt=0  log=1  eR=                  -0  eN=                  -0
x=      8000  y=     3.978873557e-05  lt=0  log=0  eR=                   0  eN=                   0
x=      9000  y=    -3.536839044e-05  lt=1  log=1  eR=     1.170620714e-12  eN=                  -0
x=      9000  y=        0.9999646322  lt=1  log=0  eR=                   0  eN=                   0
x=      9000  y=        -10.24970975  lt=0  log=1  eR=    -1.733080139e-16  eN=                  -0
x=      9000  y=     3.536776499e-05  lt=0  log=0  eR=                   0  eN=     1.915943397e-16
x=     10000  y=    -3.183149513e-05  lt=1  log=1  eR=     1.955295556e-12  eN=     2.128792113e-16
x=     10000  y=         0.999968169  lt=1  log=0  eR=    -1.110258365e-16  eN=                   0
x=     10000  y=        -10.35507026  lt=0  log=1  eR=                  -0  eN=                  -0
x=     10000  y=     3.183098851e-05  lt=0  log=0  eR=                   0  eN=     2.128825995e-16
x=    100000  y=    -3.183103928e-06  lt=1  log=1  eR=    -5.496886055e-12  eN=     1.330514125e-16
x=    100000  y=        0.9999968169  lt=1  log=0  eR=                   0  eN=                   0
x=    100000  y=        -12.65765535  lt=0  log=1  eR=    -1.403385374e-16  eN=    -1.403385374e-16
x=    100000  y=     3.183098862e-06  lt=0  log=0  eR=     1.330516242e-16  eN=     1.330516242e-16
x=   1000000  y=    -3.183099368e-07  lt=1  log=1  eR=     1.358965789e-10  eN=    -1.663145038e-16
x=   1000000  y=        0.9999996817  lt=1  log=0  eR=                   0  eN=                   0
x=   1000000  y=        -14.96024044  lt=0  log=1  eR=                  -0  eN=                  -0
x=   1000000  y=     3.183098862e-07  lt=0  log=0  eR=     1.663145303e-16  eN=                   0
x=  10000000  y=    -3.183098912e-08  lt=1  log=1  eR=    -3.352301937e-09  eN=                  -0
x=  10000000  y=        0.9999999682  lt=1  log=0  eR=      1.11022306e-16  eN=                   0
x=  10000000  y=        -17.26282554  lt=0  log=1  eR=                  -0  eN=                  -0
x=  10000000  y=     3.183098862e-08  lt=0  log=0  eR=                   0  eN=                   0
x= 100000000  y=    -3.183098867e-09  lt=1  log=1  eR=     1.059916874e-08  eN=                  -0
x= 100000000  y=        0.9999999968  lt=1  log=0  eR=                   0  eN=                   0
x= 100000000  y=        -19.56541063  lt=0  log=1  eR=                  -0  eN=                  -0
x= 100000000  y=     3.183098862e-09  lt=0  log=0  eR=     1.299332268e-16  eN=                   0
x=1000000000  y=    -3.183098862e-10  lt=1  log=1  eR=    -1.986729412e-07  eN=                  -0
x=1000000000  y=        0.9999999997  lt=1  log=0  eR=     1.110223025e-16  eN=                   0
x=1000000000  y=        -21.86799572  lt=0  log=1  eR=                  -0  eN=                  -0
x=1000000000  y=     3.183098862e-10  lt=0  log=0  eR=     1.624165335e-16  eN=     1.624165335e-16
x=     1e+10  y=    -3.183098862e-11  lt=1  log=1  eR=    -1.986729413e-07  eN=                  -0
x=     1e+10  y=                   1  lt=1  log=0  eR=                   0  eN=                   0
x=     1e+10  y=        -24.17058082  lt=0  log=1  eR=                  -0  eN=                  -0
x=     1e+10  y=     3.183098862e-11  lt=0  log=0  eR=     2.030206668e-16  eN=     2.030206668e-16
x=     1e+11  y=    -3.183098862e-12  lt=1  log=1  eR=     6.777064055e-06  eN=                  -0
x=     1e+11  y=                   1  lt=1  log=0  eR=                   0  eN=                   0
x=     1e+11  y=        -26.47316591  lt=0  log=1  eR=                  -0  eN=                  -0
x=     1e+11  y=     3.183098862e-12  lt=0  log=0  eR=                   0  eN=                   0


[Rd] pcauchy precision (PR#6756)
terra@gnome.org