beebe@math.utah.edu
2003-Aug-25 16:34 UTC
[Rd] Re: R 1.7.x and inaccurate log1p() on OpenBSD 3.2 and NetBSD 1.6 (PR#3979)
>> I have come across your reported log1p error (#2837) on a NetBSD (1.6W) >> system.I've just made further experiments on the deficient log1p() function on OpenBSD 3.2 and NetBSD 1.6 with this test program: % cat bug-log1p.c #include <stdio.h> #include <stdlib.h> #include <math.h> int main(int argc, char* argv[]) { int k; double x; for (k = 0; k <= 100; ++k) { x = pow(2.0,(double)(-k)); printf("%3d\t%.15e\t%.15e\n", k, log1p(x), log(1.0 + x)); } return (EXIT_SUCCESS); } % cc bug-log1p.c -lm && ./a.out 0 6.931471805599453e-01 6.931471805599453e-01 1 4.054651081081644e-01 4.054651081081644e-01 2 2.231435513142098e-01 2.231435513142098e-01 ... 51 4.440892098500625e-16 4.440892098500625e-16 52 2.220446049250313e-16 2.220446049250313e-16 53 0.000000000000000e+00 0.000000000000000e+00 54 0.000000000000000e+00 0.000000000000000e+00 ... 99 0.000000000000000e+00 0.000000000000000e+00 100 0.000000000000000e+00 0.000000000000000e+00 Evidently, on these systems, log1p(x) is carelessly implemented as log(1+x). Correct output from FreeBSD 5.0, Sun Solaris 9, ... looks like this: % cc bug-log1p.c -lm && ./a.out 0 6.931471805599453e-01 6.931471805599453e-01 1 4.054651081081644e-01 4.054651081081644e-01 2 2.231435513142098e-01 2.231435513142098e-01 ... 51 4.440892098500625e-16 4.440892098500625e-16 52 2.220446049250313e-16 2.220446049250313e-16 53 1.110223024625157e-16 0.000000000000000e+00 54 5.551115123125783e-17 0.000000000000000e+00 ... 99 1.577721810442024e-30 0.000000000000000e+00 100 7.888609052210118e-31 0.000000000000000e+00 The whole point of log1p(x) is to return accurate results for |x| << 1, and the OpenBSD/FreeBSD folks failed to understand that. The simple solution for a missing log1p() that I adopted in hoc is this internal function: fp_t Log1p(fp_t x) { #if defined(HAVE_LOG1PF) || defined(HAVE_LOG1P) || defined(HAVE_LOG1PL) return (log1p(x)); #else fp_t u; /* Use log(), corrected to first order for truncation loss */ u = FP(1.0) + x; if (u == FP(1.0)) return (x); else return (log(u) * (x / (u - FP(1.0)) )); #endif } I have yet to put in an accuracy test in hoc's configure.in that will check for a broken log1p(), and use the internal fallback implementation instead. Here is a test comparing accuracy of the two log1p() implementations on Sun Solaris 9, which has a good log1p() implementation: % cat cmp-log1p.c #include <stdio.h> #include <stdlib.h> #include <math.h> double LOG1P(double x) { double u; u = 1.0 + x; if (u == 1.0) return (x); else return (log(u) * (x / (u - 1.0))); } int main(int argc, char* argv[]) { int k; double d; double x; for (k = 0; k <= 100; ++k) { x = pow(2.0,(double)(-k)); printf("%3d\t%.15e\t%.15e\t%.2e\n", k, log1p(x), LOG1P(x), (LOG1P(x) - log1p(x))/LOG1P(x)); } return (EXIT_SUCCESS); } % cc cmp-log1p.c -lm && ./a.out 0 6.931471805599453e-01 6.931471805599453e-01 0.00e+00 1 4.054651081081644e-01 4.054651081081644e-01 0.00e+00 2 2.231435513142098e-01 2.231435513142098e-01 0.00e+00 ... 51 4.440892098500625e-16 4.440892098500625e-16 0.00e+00 52 2.220446049250313e-16 2.220446049250313e-16 0.00e+00 53 1.110223024625157e-16 1.110223024625157e-16 0.00e+00 54 5.551115123125783e-17 5.551115123125783e-17 0.00e+00 ... 98 3.155443620884047e-30 3.155443620884047e-30 0.00e+00 99 1.577721810442024e-30 1.577721810442024e-30 0.00e+00 100 7.888609052210118e-31 7.888609052210118e-31 0.00e+00 At least for test arguments of the form 2^(-k), my LOG1P() is identical to log1p(). A simple change to that test program, inserting x *= (double)rand() / (double)(RAND_MAX); after the assignment to x to pick a random value near a power of k, produces output like this: % cc cmp-log1p-2.c -lm && ./a.out 0 4.146697237286190e-01 4.146697237286190e-01 0.00e+00 1 8.421502722841255e-02 8.421502722841256e-02 1.65e-16 2 7.432648260535767e-02 7.432648260535767e-02 0.00e+00 ... 48 2.771522173451896e-15 2.771522173451896e-15 1.42e-16 49 1.346294923235749e-15 1.346294923235749e-15 0.00e+00 50 8.498507032336806e-16 8.498507032336806e-16 0.00e+00 51 1.246870549827746e-17 1.246870549827746e-17 0.00e+00 52 7.077345664348359e-17 7.077345664348359e-17 0.00e+00 ... 98 2.127061943360297e-30 2.127061943360297e-30 0.00e+00 99 1.276978671673724e-30 1.276978671673724e-30 0.00e+00 100 1.252374165764246e-31 1.252374165764246e-31 0.00e+00 For all random test arguments x < 2^(-49), the relative error of LOG1P() vs log1p() is zero. ------------------------------------------------------------------------------- - Nelson H. F. Beebe Tel: +1 801 581 5254 - - Center for Scientific Computing FAX: +1 801 581 4148 - - University of Utah Internet e-mail: beebe@math.utah.edu - - Department of Mathematics, 110 LCB beebe@acm.org beebe@computer.org - - 155 S 1400 E RM 233 beebe@ieee.org - - Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe -
Prof Brian Ripley
2003-Aug-25 19:19 UTC
(PR#3979) Re: [Rd] Re: R 1.7.x and inaccurate log1p() on OpenBSD 3.2 and NetBSD 1.6 (PR#3979)
There is already a usable log1p implementation in src/nmath/log1p, for platforms without it. All we need to do is to arrange to use it on those systems with broken versions. That's not easy without access to such a platform to test it, though. On Mon, 25 Aug 2003 beebe@math.utah.edu wrote:> >> I have come across your reported log1p error (#2837) on a NetBSD (1.6W) > >> system. > > I've just made further experiments on the deficient log1p() function > on OpenBSD 3.2 and NetBSD 1.6 with this test program: > > % cat bug-log1p.c > #include <stdio.h> > #include <stdlib.h> > #include <math.h> > > int > main(int argc, char* argv[]) > { > int k; > double x; > > for (k = 0; k <= 100; ++k) > { > x = pow(2.0,(double)(-k)); > printf("%3d\t%.15e\t%.15e\n", k, log1p(x), log(1.0 + x)); > } > > return (EXIT_SUCCESS); > } > > % cc bug-log1p.c -lm && ./a.out > 0 6.931471805599453e-01 6.931471805599453e-01 > 1 4.054651081081644e-01 4.054651081081644e-01 > 2 2.231435513142098e-01 2.231435513142098e-01 > ... > 51 4.440892098500625e-16 4.440892098500625e-16 > 52 2.220446049250313e-16 2.220446049250313e-16 > 53 0.000000000000000e+00 0.000000000000000e+00 > 54 0.000000000000000e+00 0.000000000000000e+00 > ... > 99 0.000000000000000e+00 0.000000000000000e+00 > 100 0.000000000000000e+00 0.000000000000000e+00 > > Evidently, on these systems, log1p(x) is carelessly implemented as > log(1+x). Correct output from FreeBSD 5.0, Sun Solaris 9, ... looks > like this: > > % cc bug-log1p.c -lm && ./a.out > 0 6.931471805599453e-01 6.931471805599453e-01 > 1 4.054651081081644e-01 4.054651081081644e-01 > 2 2.231435513142098e-01 2.231435513142098e-01 > ... > 51 4.440892098500625e-16 4.440892098500625e-16 > 52 2.220446049250313e-16 2.220446049250313e-16 > 53 1.110223024625157e-16 0.000000000000000e+00 > 54 5.551115123125783e-17 0.000000000000000e+00 > ... > 99 1.577721810442024e-30 0.000000000000000e+00 > 100 7.888609052210118e-31 0.000000000000000e+00 > > The whole point of log1p(x) is to return accurate results for > |x| << 1, and the OpenBSD/FreeBSD folks failed to understand that. > > The simple solution for a missing log1p() that I adopted in hoc is > this internal function: > > fp_t > Log1p(fp_t x) > { > #if defined(HAVE_LOG1PF) || defined(HAVE_LOG1P) || defined(HAVE_LOG1PL) > return (log1p(x)); > #else > fp_t u; > /* Use log(), corrected to first order for truncation loss */ > u = FP(1.0) + x; > if (u == FP(1.0)) > return (x); > else > return (log(u) * (x / (u - FP(1.0)) )); > #endif > } > > I have yet to put in an accuracy test in hoc's configure.in that will > check for a broken log1p(), and use the internal fallback > implementation instead. > > Here is a test comparing accuracy of the two log1p() implementations > on Sun Solaris 9, which has a good log1p() implementation: > > % cat cmp-log1p.c > #include <stdio.h> > #include <stdlib.h> > #include <math.h> > > double > LOG1P(double x) > { > double u; > > u = 1.0 + x; > if (u == 1.0) > return (x); > else > return (log(u) * (x / (u - 1.0))); > } > > > int > main(int argc, char* argv[]) > { > int k; > double d; > double x; > > for (k = 0; k <= 100; ++k) > { > x = pow(2.0,(double)(-k)); > > printf("%3d\t%.15e\t%.15e\t%.2e\n", > k, log1p(x), LOG1P(x), (LOG1P(x) - log1p(x))/LOG1P(x)); > } > > return (EXIT_SUCCESS); > } > > % cc cmp-log1p.c -lm && ./a.out > 0 6.931471805599453e-01 6.931471805599453e-01 0.00e+00 > 1 4.054651081081644e-01 4.054651081081644e-01 0.00e+00 > 2 2.231435513142098e-01 2.231435513142098e-01 0.00e+00 > ... > 51 4.440892098500625e-16 4.440892098500625e-16 0.00e+00 > 52 2.220446049250313e-16 2.220446049250313e-16 0.00e+00 > 53 1.110223024625157e-16 1.110223024625157e-16 0.00e+00 > 54 5.551115123125783e-17 5.551115123125783e-17 0.00e+00 > ... > 98 3.155443620884047e-30 3.155443620884047e-30 0.00e+00 > 99 1.577721810442024e-30 1.577721810442024e-30 0.00e+00 > 100 7.888609052210118e-31 7.888609052210118e-31 0.00e+00 > > At least for test arguments of the form 2^(-k), my LOG1P() is > identical to log1p(). > > A simple change to that test program, inserting > > x *= (double)rand() / (double)(RAND_MAX); > > after the assignment to x to pick a random value near a power of k, > produces output like this: > > % cc cmp-log1p-2.c -lm && ./a.out > 0 4.146697237286190e-01 4.146697237286190e-01 0.00e+00 > 1 8.421502722841255e-02 8.421502722841256e-02 1.65e-16 > 2 7.432648260535767e-02 7.432648260535767e-02 0.00e+00 > ... > 48 2.771522173451896e-15 2.771522173451896e-15 1.42e-16 > 49 1.346294923235749e-15 1.346294923235749e-15 0.00e+00 > 50 8.498507032336806e-16 8.498507032336806e-16 0.00e+00 > 51 1.246870549827746e-17 1.246870549827746e-17 0.00e+00 > 52 7.077345664348359e-17 7.077345664348359e-17 0.00e+00 > ... > 98 2.127061943360297e-30 2.127061943360297e-30 0.00e+00 > 99 1.276978671673724e-30 1.276978671673724e-30 0.00e+00 > 100 1.252374165764246e-31 1.252374165764246e-31 0.00e+00 > > For all random test arguments x < 2^(-49), the relative error of > LOG1P() vs log1p() is zero. > > ------------------------------------------------------------------------------- > - Nelson H. F. Beebe Tel: +1 801 581 5254 - > - Center for Scientific Computing FAX: +1 801 581 4148 - > - University of Utah Internet e-mail: beebe@math.utah.edu - > - Department of Mathematics, 110 LCB beebe@acm.org beebe@computer.org - > - 155 S 1400 E RM 233 beebe@ieee.org - > - Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe - > > ______________________________________________ > R-devel@stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-devel > >-- Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595