From ripley at stats.ox.ac.uk  Fri Aug 29 22:03:28 2003
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Fri Aug 29 22:08:13 2003
Subject: (PR#3979) Re: [Rd] Re: R 1.7.x and inaccurate log1p() on OpenBSD
	3.2 and NetBSD 1.6 (PR#3979)
In-Reply-To: <Pine.LNX.4.44.0308251809580.5741-100000@gannet.stats>
Message-ID: <Pine.LNX.4.44.0308292101500.28090-100000@gannet.stats>

The fix I proposed is now in the R-devel sources, if you could please try 
it.  (It should appear in the next snapshot after now.)

On Mon, 25 Aug 2003, Prof Brian Ripley wrote:

> 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


(PR#3979) Re: [Rd] Re: R 1.7.x and inaccurate log1p() on OpenBSD 3.2 and NetBSD 1.6 (PR#3979)
Prof Brian Ripley