>> 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: [EMAIL PROTECTED] -
- Department of Mathematics, 110 LCB [EMAIL PROTECTED] [EMAIL PROTECTED] -
- 155 S 1400 E RM 233 [EMAIL PROTECTED] -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe -
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
