#include <ncbi.h>
#include <gishlib.h>

/*
	Return values accurate to 16 digits (when such hardware precision
	is available) for the quantity log(x+1) for all valid values for x,
	large and small.
*/
double _cdecl
fct_log1p(x)
	double	x;
{
	int	i;
	double	sum, y;

	if (ABS(x) >= 0.2)
		return log(x+1.);

	for (i=0, sum=0., y=x; ; ) {
		sum += y/++i;
		if (y < 1.e-16)
			break;
		y *= x;
		sum -= y/++i;
		if (y < 1.e-16)
			break;
		y *= x;
	}
	return sum;
}