/*
 * 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_em1.c	1.1 (BeEF) 1/5/95";
#endif
/*
 * er_em1(.) computes the error committed by expm1(x) to within 0.0625 ULPs.
 * (Bound proved: 0.052 ULPs).
 *
 * er_exp(.) computes the error committed by exp(x) to within 0.03125 ULPs.
 * (Bound proved: 0.028 ULPs).
 *
 * extern GENERIC ld_2_n(.), getulp(.);
 * extern void unpack(.), double_check(.);
 */
#ifndef vms
#include <stdio.h>
#include <assert.h>
#else	/* vms */
#define assert(A)
#endif	/* vms */
#include "sys.h"
extern int sig_bits,
	mid_exp,				/* declared in driver.c */
	PCS_LOG;				/* declared in er_log.c */
extern LONG ln2_p[];				/* declared in er_log.c */
extern FLOAT brk_em1[];

static GENERIC log2_h[2],			/* 0: for n < 0 (rndd down) */
	log2_l[2];				/* 1: for n > 0 (rndd up) */

#define	TRM_EM1	16
#define MID_EM1 7

static GENERIC off_em1;
static int trm_em1;				/* # of terms needed */
static GENERIC A_em1[TRM_EM1],			/* MUST be of type GENERIC */
	      B_em1[TRM_EM1];

static FLOAT x_em1[] = {			/* centers */
    -9.8730468750e-01,		/*  0:-1011/1024 */
    -8.8574218750e-01,		/*  1: -907/1024 */
    -7.7539062500e-01,		/*  2: -794/1024 */
    -6.5332031250e-01,		/*  3: -669/1024 */
    -5.1074218750e-01,		/*  4: -523/1024 */
    -3.1835937500e-01,		/*  5: -326/1024 */
    -1.7382812500e-01,		/*  6: -178/1024 */
     0.0000000000e+00,		/*  7:=  MID_EM1 */
     2.0996093750e-01,		/*  8:  215/1024 */
     3.9746093750e-01,		/*  9:  407/1024 */
     5.9765625000e-01,		/* 10:  612/1024 */
     7.3144531250e-01,		/* 11:  749/1024 */
     9.2773437500e-01,		/* 12:  950/1024 */
};
#define CTR_EM1	(sizeof(x_em1) / sizeof(x_em1[0]))

static struct rational xpq_em1[ ] = {
	{-1011, 1024},		/*  0: -9.8730468750e-01 */
	{ -907, 1024},		/*  1: -8.8574218750e-01 */
	{ -794, 1024},		/*  2: -7.7539062500e-01 */
	{ -669, 1024},		/*  3: -6.5332031250e-01 */
	{ -523, 1024},		/*  4: -5.1074218750e-01 */
	{ -326, 1024},		/*  5: -3.1835937500e-01 */
	{ -178, 1024},		/*  6: -1.7382812500e-01 */
	{    0, 1024},		/*  7:  0.0000000000e+00 */
	{  215, 1024},		/*  8:  2.0996093750e-01 */
	{  407, 1024},		/*  9:  3.9746093750e-01 */
	{  612, 1024},		/* 10:  5.9765625000e-01 */
	{  749, 1024},		/* 11:  7.3144531250e-01 */
	{  950, 1024},		/* 12:  9.2773437500e-01 */
};

static GENERIC y_em1l[CTR_EM1],
	       y_em1h[CTR_EM1];

static LONG y_em1p[] = {
    -1,  -2, 0x1413d4, 0x0950b, 0x7b4c3, 0xe34ef, 0xc7da7,
              0xdbf87, 0x2c994, 0x512b8, 0x0d1f7, 0x59f82,
    -1,  -2, 0x12cd8d, 0xca033, 0x0ac5a, 0xedce7, 0x35895,
              0xa4579, 0x310f8, 0x43e63, 0x8a503, 0xb3a69,
    -1,  -2, 0x114363, 0x7aa69, 0xd7147, 0x64c22, 0xca981,
              0x88c8e, 0xb1c00, 0x86744, 0x6b561, 0x58f54,
    -1,  -4, 0x1eb327, 0x78e8d, 0x2b51e, 0x7d5a3, 0xcc724,
              0x75cda, 0xefc47, 0xe6661, 0x0d794, 0xa1115,
    -1,  -4, 0x1998c7, 0x9df3c, 0xf97fa, 0x29bf1, 0xdf7d0,
              0x481c1, 0xbc1b3, 0x0343a, 0xed5ad, 0x368ce,
    -1,  -4, 0x11733d, 0x40ce8, 0x48436, 0x7029c, 0x16f1e,
              0x4828b, 0xd1eb9, 0x41b6f, 0xf3565, 0x8b6f9,
    -1,  -8, 0x146c6b, 0x159f3, 0x46316, 0x85128, 0x1f1a1,
              0x37717, 0xbad10, 0x8825f, 0xbe89c, 0xdf672,
     1,   1, 0x000000, 0x00000, 0x00000, 0x00000, 0x00000,
              0x00000, 0x00000, 0x00000, 0x00000, 0x00000,
     1,  -8, 0x1de795, 0x66421, 0xdf785, 0x06eae, 0x67666,
              0x00e05, 0x0f97d, 0xd027f, 0xfc973, 0x0e5b3,
     1,  -4, 0x1f3c13, 0x1cdb9, 0x90e6b, 0xcaff6, 0x520c2,
              0xab2cf, 0x74628, 0x646ea, 0x2d158, 0x0a401,
     1,  -2, 0x1a2bda, 0x7ecfc, 0xf7660, 0x768cb, 0xc27b0,
              0x815ef, 0x73298, 0xf4dac, 0xe2ea6, 0x557a3,
     1,   1, 0x113fd2, 0xd2ba8, 0x5bdd9, 0x7f284, 0xf8128,
              0xbd785, 0xabff5, 0x49626, 0x4c8b1, 0x8ff71,
     1,   1, 0x1875db, 0x20de2, 0x3988d, 0x218a0, 0x96305,
              0x544ca, 0xb9efa, 0xee1f6, 0x6611c, 0xbd1e5,
};
#define PCS_EM1 ((sizeof(y_em1p) / sizeof(y_em1p[0])) / CTR_EM1)

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

	for (i = 0; i < CTR_EM1; i++)
		double_check(i, "x_em1", &xpq_em1[i], &x_em1[i]);
}

static GENERIC
#ifdef __STDC__
_A_em1(int m)
#else
_A_em1(m)
int m;
#endif
{
	LONG t = 4 * m - 3;

	return -SIX/ (GENERIC) (t * (t + 4));
}

static GENERIC
#ifdef __STDC__
_B_em1(int m)
#else
_B_em1(m)
int m;
#endif
{
	LONG t = 4 * m + 1;

	return -NINE/ (GENERIC) ((t - 2) * t * t * (t + 2));
}

void
#ifdef __STDC__
bgn_em1(void)
#else
bgn_em1( )
#endif
{
	int i, pieces;

	trm_em1 = (sig_bits + 12) / 16 + 1;		/* (m + 4 + 13) / 26 + 1 will do */
	off_em1 = ld_2_n(-sig_bits / 2 - 2);

	A_em1[0] = -ONE; A_em1[0] /= FIVE;	/* for sin / cos */
	for (i = 1; i <= trm_em1; i++) {
		A_em1[i] = _A_em1(i);
		B_em1[i] = _B_em1(i);
	}
	pieces = (sig_bits + 15) / 20 + 1;

	unpack(pieces, (int) PCS_EM1, (int) (CTR_EM1 - 1), y_em1p, y_em1l, y_em1h);
}

void
#ifdef __STDC__
bgn_exp(void)
#else
bgn_exp( )
#endif
{
	int pieces;
	GENERIC gtmp, chopped, ulp, scale;

	pieces = sig_bits / 10 + 1;
	unpack(pieces, (int) (PCS_LOG + 5), (int) 0, ln2_p, log2_l, log2_h);

	ulp = getulp(&log2_h[0], &scale);
	assert(scale == ONE);
	if (log2_l[0] < ZERO) {			/* rndd up (to t bits)? */
		log2_h[0] -= ulp;		/* chop it */
		log2_l[0] += ulp;
	}
	assert(log2_l[0] >= ZERO);

	ulp *= (GENERIC) BITS_2M;		/* ulp (to t - m bits) */

	gtmp = log2_h[0]* (GENERIC) BITS_2M;
	chopped = gtmp + log2_h[0];
	chopped -= gtmp;
	if (chopped > log2_h[0])		/* rndd up (to t - m bits)? */
		chopped -= ulp;				/* chop it */
	assert(chopped <= log2_h[0]);

	log2_l[0] += log2_h[0] - chopped;
	assert(log2_l[0] > ZERO);
	log2_h[0] = chopped;			/* rndd down (to t - m bits) */

	log2_l[1] = log2_l[0] - ulp;
	assert(log2_l[1] < ZERO);
	log2_h[1] = chopped + ulp;		/* rndd up (to t - m bits) */
}

#if TEST_EXP || TEST_EM1
static GENERIC
#ifdef __STDC__
cf_tan(GENERIC *_u)
#else
cf_tan(_u)
GENERIC *_u;			/* input argument */
#endif
{
	int m;
	GENERIC u = *_u, y = u, gtmp;

	for (m = trm_em1; m >= 1; m--) {
		gtmp = A_em1[m] + B_em1[m] / y;
		y = u + gtmp;
	}
	return y;
}

static GENERIC
#ifdef __STDC__
R_em1(GENERIC *_x)
#else
R_em1(_x)			/* returns expm1(x) - x - x^2 / 2 */
GENERIC *_x;			/* input argument */
#endif
{
	GENERIC x = *_x, hx, cf, gtmp;

	if (x <= off_em1 && x >= -off_em1)
		return ZERO;
	hx = x * HALF;
	gtmp = -THREE / (hx * hx); cf = hx / cf_tan(&gtmp);
	gtmp = cf + cf;
	gtmp += hx * hx * x;
	return (gtmp + (x + x * hx) * cf) / (ONE - (hx + cf));
}

static GENERIC
#ifdef __STDC__
raw_em1(int rid, GENERIC *_x, GENERIC *_y)
#else
raw_em1(rid, _x, _y)
int rid;
GENERIC *_x, *_y;		/* input arguments */
#endif
{
	GENERIC x = *_x, y = *_y, xi, hi, lo, error;

	if (rid == MID_EM1)
		error = ((y - x) - HALF * x * x) - R_em1(&x);
	else {
		xi = x- (GENERIC) x_em1[rid];
		hi = HALF * xi * xi;
		lo = R_em1(&xi);
		error = (((((((((y - y_em1h[rid]) - xi) - y_em1h[rid] * xi) -
			hi) - y_em1h[rid] * hi) - lo) - y_em1h[rid] * lo) -
			y_em1l[rid]) - y_em1l[rid] * xi) - y_em1l[rid] * hi;
	}
	return error;
}

#if TEST_EM1
GENERIC
#ifdef __STDC__
er_em1(int rid, GENERIC *_x, GENERIC *_y)
#else
er_em1(rid, _x, _y)
int rid;
GENERIC *_x, *_y;		/* input arguments */
#endif
{
	GENERIC x = *_x, y = *_y, error, ulp, scale;

	if (y == ZERO) return ZERO;
	error = raw_em1(rid, &x, &y);

	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

#if TEST_EXP
static GENERIC
#ifdef __STDC__
red_exp(int *rid, GENERIC *x, GENERIC *y)
#else
red_exp(rid, x, y)				/* returns expm1(-n * log2_l) */
int *rid;
GENERIC *x, *y;
#endif
{
	static int nold = 0;			/* cache n in nold */
	static GENERIC p2_mn,			/* cache 2^-n */
		nlog2h,				/* cache n*log2_h[i] */
		z;				/* cache E(-n*log2_l[i]) */
	int n, i;
	GENERIC gtmp;

	if (*x < ZERO)
		n = *x / log2_h[0] - HALF;		/* assume truncated */
	else
		n = *x / log2_h[1] + HALF;		/* assume truncated */
	assert(n != 0 && n > -BITS_2M && n < BITS_2M);
	if (n != nold) {
		nold = n;
		i = n < 0 ? 0 : 1;
		p2_mn = ld_2_n(-n);
		nlog2h = log2_h[i] * (GENERIC) n;
		gtmp = -log2_l[i] * (GENERIC) n;
		z = R_em1(&gtmp);
		z += gtmp * gtmp * HALF;
		z += gtmp;
	}
	*y *= p2_mn;
	*x -= nlog2h;

	*rid = mid_exp;
	if (*x < (GENERIC) brk_em1[*rid]) {
		while (*x < (GENERIC) brk_em1[--(*rid)])
			;
		assert(*rid >= mid_exp - 2);
	}
	else {
		while (*x >= (GENERIC) brk_em1[++(*rid)])
			;
		(*rid)--;
		assert(*rid <= mid_exp + 2);
	}
	return z;
}

GENERIC
#ifdef __STDC__
er_exp(int rid, GENERIC *_x, GENERIC *_y)
#else
er_exp(rid, _x, _y)
int rid;
GENERIC *_x, *_y;		/* input arguments */
#endif
{
	GENERIC x = *_x, y = *_y, z = ZERO, error, ulp, scale, gtmp;

	if (rid < 0) {
		z = red_exp(&rid, &x, &y);
		assert(z >= ZERO);
	}
	gtmp = y - ONE;
	error = raw_em1(rid, &x, &gtmp);
	if (z != ZERO) {
		gtmp -= error;			/* approximated E(xi) */
		error -= z;
		error -= gtmp * z;
	}

	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_EXP */
#endif	/* TEST_EXP || TEST_EM1 */

#if TEST_COS || TEST_SIN
void
#ifdef __STDC__
cf_tan2(GENERIC *_u, GENERIC *U, GENERIC *V)
#else
cf_tan2(_u, U, V)					/* USED BY er_sin(.) */
GENERIC *_u, *U, *V;		/* _u is input argument */
#endif
{
	int m;
	GENERIC u = *_u, y = u, gtmp;

	for (m = trm_em1; m >= 2; m--) {
		gtmp = A_em1[m] + B_em1[m] / y;
		y = u + gtmp;
	}
	y = B_em1[1] / y;
	*U = A_em1[0] + y;
	*V = A_em1[1] + y;
}
#endif	/* TEST_COS || TEST_SIN */
