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

[beebe at math.utah.edu]

Brian Ripley writes today:

>> 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.

I need the same kind of test in my own software, so I made some
experiments and came up with the code below for configure.in.  I've
tested it so far only on NetBSD 1.6 and Sun Solaris 9; it produced the
expected settings

	#define HAVE_BROKEN_LOG1P 1
	#define HAVE_LOG1P 1

on NetBSD and

	/* #undef HAVE_BROKEN_LOG1P */
	#define HAVE_LOG1P 1

on Solaris.  Notice that the test code checks arguments over a range
of powers of two suitable for IEEE 754 with support for subnormals,
but will bail out early if necessary.  It should therefore work as
intended on VAX OpenVMS and IBM S/390; I only rarely have access to
such systems.

AH_TEMPLATE(HAVE_BROKEN_LOG1P,		[Define if log1p() is broken.])
AH_TEMPLATE(HAVE_LOG1P,			[Define if log1p() is available.])
...
AC_SEARCH_LIBS(log1p,          [m], AC_DEFINE(HAVE_LOG1P))
...

dnl Check for broken log1p() implementations
if test "$ac_cv_search_log1p" = "none required"
then
	AC_MSG_CHECKING(for apparently-accurate log1p())
	AC_TRY_RUN(
	[
#include <stdio.h>
#include <math.h>

int
main(void)
{
    int k;
    double d;
    double x;

    /* log(1+x) = x - (1/2)x^2 + (1/3)x^3 - (1/4)x^4 ... */
    /*          = x for x sufficiently small */
    for (k = -54; k > -1074; --k)
    {	
	x = pow(2.0, (double)(k));
	if (x == 0.0)
	    return (0);			/* OK: reached underflow limit */
	d = log1p(x);
	if (d == 0.0)
	    return (1);			/* ERROR: inaccurate log1p() */
	if (d != x)
	    return (1);			/* ERROR: inaccurate log1p() */
    }	
    
    return (0);
}
	],
	[
	    AC_MSG_RESULT(yes)
	],
	[
	    AC_MSG_RESULT(no)
	    AC_DEFINE(HAVE_BROKEN_LOG1P)
	]
	)
fi

-------------------------------------------------------------------------------
- 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  -