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

[beebe at math.utah.edu]

>> 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 at math.utah.edu  -
- Department of Mathematics, 110 LCB        beebe at acm.org  beebe at computer.org -
- 155 S 1400 E RM 233                       beebe at ieee.org                    -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe  -