/*
 * Copyright 1/5/95 Zhishun Alex Liu.  All Rights Reserved.
 *
 ***************************
 * WORK IN PROGRESS
 ***************************
 *
 * Written by Zhishun Alex Liu under direction of Professor William Kahan.
 * The author's current electronic mail address as of December 1994 is
 * Alex.Liu@Eng.Sun.COM.
 *
 * Use of this software is granted with the understanding that all recipients
 * should regard themselves as participants in an ongoing research project and
 * hence should feel obligated to report their experiences (good or bad) with
 * this elementary functions test suite to the author.
 */
#ifndef lint
static char sccsid[] = "@(#)er_l1p.c	1.1 (BeEF) 1/5/95";
#endif
/*
 * er_l1p(.) computes the error committed by log1p(x) to with 0.0625 ULPs.
 * (Bound proved: 0.052 ULPs).
 * extern GENERIC getulp(.), cf_atn(.);
 * extern void double_check(.);
 */
#include "sys.h"
extern GENERIC off_log;				/* computed in bgn_log(.) */
extern GENERIC y_ln_l[], y_ln_h[];		/* unpacked in bgn_log(.) */

static FLOAT x_l1p[] = {
    -2.75390625e-01,				/*  0 */
    -2.30468750e-01,				/*  1 */
    -2.07031250e-01,				/*  2 */
    -1.83593750e-01,				/*  3 */
    -1.60156250e-01,				/*  4 */
    -1.36718750e-01,				/*  5 */
     0.00000000e+00,				/*  6 */
     1.32812500e-01,				/*  7 */
     1.58203125e-01,				/*  8 */
     1.93359375e-01,				/*  9 */
     2.28515625e-01,				/*  0 */
     2.65625000e-01,				/* 11 */
     3.24218750e-01,				/* 12 */
     4.08203125e-01,				/* 13 */
};
#define CTR_L1P (sizeof(x_l1p) / sizeof(x_l1p[0]))

static struct rational xpq_l1p[ ] = {
	{ -282, 1024},		/*  0: -2.7539062500e-01 */
	{ -236, 1024},		/*  1: -2.3046875000e-01 */
	{ -212, 1024},		/*  2: -2.0703125000e-01 */
	{ -188, 1024},		/*  3: -1.8359375000e-01 */
	{ -164, 1024},		/*  4: -1.6015625000e-01 */
	{ -140, 1024},		/*  5: -1.3671875000e-01 */
	{    0, 1024},		/*  6:  0.0000000000e+00 */
	{  136, 1024},		/*  7:  1.3281250000e-01 */
	{  162, 1024},		/*  8:  1.5820312500e-01 */
	{  198, 1024},		/*  9:  1.9335937500e-01 */
	{  234, 1024},		/* 10:  2.2851562500e-01 */
	{  272, 1024},		/* 11:  2.6562500000e-01 */
	{  332, 1024},		/* 12:  3.2421875000e-01 */
	{  418, 1024},		/* 13:  4.0820312500e-01 */
};

void
#ifdef __STDC__
dbc_l1p(void)
#else
dbc_l1p( )
#endif
{
	int i;

	for (i = 0; i < CTR_L1P; i++)
		double_check(i, "x_l1p", &xpq_l1p[i], &x_l1p[i]);
}

#if TEST_L1P
/* NOTE: In order to run this routine, A_atn[], B_atn[] & trm_atn
 *      must be initialized (currently they are initialized in bgn_atn( ).
 *      This routine also calls cf_atn( ).
 */
GENERIC
#ifdef __STDC__
er_l1p(int i, GENERIC *_x, GENERIC *_y)
#else
er_l1p(i, _x, _y)
GENERIC *_x, *_y;
#endif
{
	GENERIC x = *_x, y = *_y, hxsq, third, rho, scale, error, ulp;

#if GUARD_AGAINST_OVERFLOW
	if (y == ZERO)
		return ZERO;
#endif
	x = (x - (GENERIC) x_l1p[i]) / (GENERIC) (ONE + x_l1p[i]);
	hxsq = x * x * HALF;
	third = hxsq * x; third /= x + TWO;
	rho = x - hxsq; rho += third;
	error = (y - y_ln_h[i]) - x; error += hxsq; error -= third;
	if (rho > off_log || rho < -off_log) {
		ulp = (GENERIC) (-12) / (rho * rho);
		error += rho / cf_atn(&ulp);
	}
	error -= y_ln_l[i];

	ulp = y - error; ulp = getulp(&ulp, &scale);
	error *= scale;
#if GUARD_AGAINST_OVERFLOW
	if (GABS(error) / GHUGE > ulp)
		return ZERO;
#endif
	return error / ulp;
}
#endif	/* TEST_L1P */
