/*
 * 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_sin.c	1.1 (BeEF) 1/5/95";
#endif
/*
 * er_sin(.) and er_cos(.) compute the error committed by sin(x) and cos(x)
 * to within 0.0625 ULPs.  (Bounds proved: 0.0611 ULPs).
 *
 * Requires a C compiler that distinguishes cases among identifiers.
 * extern GENERIC ld_2_n(.), getulp(.);
 * extern void unpack(.), cf_tan2(.), double_check(.);
 */
#ifdef __STDC__
extern void exit(int);
#endif
#include <stdio.h>
#include <assert.h>
#include "sys.h"
extern int sig_bits;

#define BRK_CS	3
#define VAL_CS	(BRK_CS + BRK_CS)

static int ntl_sin;	/* #terms in the Taylor expansion of sin(x) - x + x^3 / 6 */
static int ntm_sin[] = {
	 -4,	  1,	  9,	 17,	 26,	 35,	 45,	 55,
	 65,	 76,	 87,	 98,	110,	121,	133,	145,
	157,	170,	182,	195,	208,	221,	234,	247,
};
#define NTM_SIN	(sizeof(ntm_sin) / sizeof(ntm_sin[0]))

static GENERIC B_sin[NTM_SIN];	/* B[k] := (2k + 4)(2k + 5), 1 <= k <= ntl_sin - 1 */

static GENERIC p2m20;		/* 2^-20, used by mkjpi2() and mknpi2() */
static GENERIC off_sin;
static GENERIC val_l[VAL_CS + 1], val_h[VAL_CS + 1];
static GENERIC *theta_h, *sine_h, *theta_l, *sine_l;
static GENERIC adjust[BRK_CS];
static FLOAT cosine[] = {
	0.875,
	0.75,
	-0.5,
};

static struct rational cospq[ ] = {
	{  896, 1024},		/*  0:  8.7500000000e-01 */
	{  768, 1024},		/*  1:  7.5000000000e-01 */
	{ -512, 1024},		/*  2: -5.0000000000e-01 */
};

static LONG val_p[] = {					   /* mod 17 */
     1,  -2, 0x102be9, 0xce0b8, 0x7cd1e, 0x5d09d, 0xa2e0f, /* acos(7 / 8) */
              0x0423b, 0xfa5d4, 0x8cd07, 0xedf3e, 0x3a37b,
              0x68b74, 0xc4714, 0x437a3, 0x8ccc3, 0x5d087,
     1,  -2, 0x1720a3, 0x92c1d, 0x95485, 0x1badb, 0xd6cd2, /* acos(3 / 4) */
              0xd8e79, 0x2c125, 0x263fd, 0xd2063, 0xde0d7,
              0xa0d4c, 0xf8da6, 0x36f57, 0xabf3d, 0x7dc51,
     1,  -2, 0x10c152, 0x382d7, 0x36584, 0x65bb3, 0x2e0f5, /* %pi / 6 */
              0x67ad1, 0x16e15, 0x8680b, 0x63351, 0x09aad,
              0x64fe3, 0x2f96f, 0x79831, 0x70d60, 0xa212f,
     1,  -4, 0x1efbde, 0xb14f4, 0xed9b1, 0x7ae80, 0x7907e, /* sqrt(15) / 8 */
              0x1e88f, 0xe92c4, 0xda189, 0x48f5b, 0x8c75f,
              0x950f3, 0xed92b, 0xbbd30, 0xb2a35, 0x9506c,
     1,  -2, 0x152a7f, 0xa9d2f, 0x8e9b7, 0x8e753, 0xf30fe, /* sqrt(7) / 4 */
              0x1bd10, 0x6aa53, 0xbff3d, 0xbc20a, 0xefe47,
              0x1a537, 0xf4716, 0xfc952, 0x72da0, 0x464e4,
     1,  -2, 0x1bb67a, 0xe8584, 0xcaa73, 0xb2574, 0x2d707, /* sqrt(3) / 2 */
              0x8b83b, 0x8925d, 0x834cc, 0x53da4, 0x798c7,
              0x20a64, 0x86e45, 0xa6e24, 0x90bcf, 0xd95ef,
};
#define PCS_CS ((sizeof(val_p) / sizeof(val_p[0])) / VAL_CS)

static LONG piov2_p[] = {	/* no leading sign + power-of-2 info */
	0x1921fb, 0x54442, 0xd1846, 0x9898c, 0xc5170,
	0x1b839, 0xa2520, 0x49c11, 0x14cf9, 0x8e804,
	0x177d4, 0xc7627, 0x3644a, 0x29410, 0xf31c6,
	0x809bb, 0xdf2a3, 0x3679a, 0x74863, 0x66056,
	0x14dbe, 0x4be28, 0x6e9fc, 0x26ada, 0xdaa38,
	0x48bc9, 0x0b6ae, 0xcc4bc, 0xfd8de, 0x89885,
	0xd34c6, 0xfdad6,
};
#define	PCS_PI2	(sizeof(piov2_p) / sizeof(piov2_p[0]))

#define	H	0
#define	M	1
#define	L	2
#define	HML	3
#define	MAX_PI2	4

static GENERIC npiov2[PCS_PI2 + HML * (MAX_PI2 + 1) + 1];	/* i * %pi / 2, i = 0..MAX_PI2 */
#define	PI2x(i)	(npiov2 + HML * i)

#define	APRX_SIN(y, x)	FUNC_SIN(y, x)
#define	APRX_COS(y, x)	FUNC_COS(y, x)

#if MACHINE_PI
static LONG mpi2_p[] = {	/* no leading sin + power-of-2 info */
	MPI2_P
};
#define	PCS_MP2	(sizeof(mpi2_p) / sizeof(mpi2_p[0]))

#if BEEF_QA
GENERIC mpi2[HML];		/* used by hXxx_() in BEEF_QA mode */
#endif
GENERIC mpi_mu;			/* (true_pi-machine_PI) / true_pi_rounded */
				/* also used in er_atn.c */
#endif	/* MACHINE_PI */

static int
#ifdef __STDC__
mkjpi2(int j, int chunks, LONG *lp)
#else
mkjpi2(j, chunks, lp)
int j, chunks;
LONG *lp;
#endif
{
	register int i;
	GENERIC p2m20i;

	p2m20i = ONE;
	for (i = 0; i < chunks; i++) {
		PI2x(j)[i] = (GENERIC) ((LONG) j * lp[i]);
		PI2x(j)[i] *= p2m20i;			/* to avoid underflow */
		PI2x(j)[i] *= p2m20;			/* we must "*" twice */
		p2m20i *= p2m20;
	}
	chunks = distd(chunks, PI2x(j));
	for (i = chunks; i <= L; i++)
		PI2x(j)[i] = ZERO;
	return chunks;
}

static void
#ifdef __STDC__
mknpi2(void)
#else
mknpi2()
#endif
{
	register int chunks, i, j;
	LONG *lp;
#if MACHINE_PI
	GENERIC pi2[2 * HML];
#endif

	p2m20 = ld_2_n(-20);
	chunks = (sig_bits * HML - 1 + (20 - 1)) / 20;	/* need HML good pieces */
	i = mkjpi2(0, chunks, piov2_p);
	assert(i == 0);
	i = mkjpi2(1, chunks, piov2_p);
	assert(i == HML || i == HML + 1);
#if MACHINE_PI
	for (i = 0; i <= L; i++)
		pi2[i] = PI2x(1)[i];
	chunks = PCS_MP2;
	lp = mpi2_p;
	j = 1;
#else
	lp = piov2_p;
	j = 2;
#endif
	for (; j <= MAX_PI2; j++) {
		i = mkjpi2(j, chunks, lp);
#if MACHINE_PI
		assert(i >= 1);
#else
		assert(i == HML || i == HML + 1);
#endif
	}
#if MACHINE_PI
	assert(PI2x(1)[H] == pi2[H]);
	for (i = 0; i <= L; i++)
		pi2[HML + i] = -PI2x(1)[i];
	i = distd(2 * HML, pi2);
	assert(i >= 1);
	mpi_mu = pi2[H] / PI2x(1)[H];
#if BEEF_QA
	for (i = 0; i <= L; i++)
		mpi2[i] = PI2x(1)[i];
#endif	/* BEEF_QA */
#endif	/* MACHINE_PI */
}

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

	for (i = 0; i < 3; i++)
		double_check(i, "cosine", &cospq[i], &cosine[i]);
}

/* NOTE that in order to run this routine, A_em1[], B_em1[] & trm_em1
 *      must be initialized (currently they are initialized in bgn_em1(.)).
 */
void
#ifdef __STDC__
bgn_sin(void)
#else
bgn_sin( )
#endif
{
	int k = 2, pieces;

	off_sin = ld_2_n(-sig_bits / 2 - 2);
	pieces = (sig_bits + sig_bits + 10) / 20;
	for (k = 2; k < NTM_SIN; k++)
		if (sig_bits <= ntm_sin[k])
			break;
	if (k == NTM_SIN) {
		(void) printf(
			"Botch in bgn_sin(.): ntl_sin cannot be determined.\n");
		exit(0);
	}
	ntl_sin = k;
	for (k = 1; k < ntl_sin; k++)		/* initialize B_sin */
		B_sin[k] = (GENERIC) ((k + k + 4) * (k + k + 5));

	unpack(pieces, (int) PCS_CS, (int) (VAL_CS - 1), val_p, val_l, val_h);

	theta_h = val_h; sine_h = val_h + BRK_CS;
	theta_l = val_l; sine_l = val_l + BRK_CS;
	for (k = 0; k < BRK_CS; k++)
		adjust[k] = sine_l[k] - (GENERIC) cosine[k] * theta_l[k];

	mknpi2();
}

#if TEST_COS || TEST_SIN
				/* 1: increasing, 0: decreasing */
#define	IORD(rid, c, x)	{\
	switch ((rid + c) % 4) {\
	case 0: *id = 1; break;\
	case 1: *id = x <= PI2x(rid)[H]; break;\
	case 2: *id = 0; break;\
	case 3: *id = x > PI2x(rid)[H]; break;\
	default: exit(1); break;\
	}\
}

static void
#ifdef __STDC__
RcRs(GENERIC *_x, GENERIC *_hx2, GENERIC *Rc, GENERIC *Rs)
#else
RcRs(_x, _hx2, Rc, Rs)	/* cos(x) - 1 + x^2 / 2 in Rc */
			/* sin(x) - x in Rs if Rs non-null */
GENERIC *_x, *_hx2, *Rc, *Rs;	/* _x and _hx2 are input arguments */
#endif
{				/* hx2 := x^2 / 2; hx_2 := (x / 2)^2 */
	GENERIC x = *_x, hx2 = *_hx2, hx_2, U, V, denum, gtmp;

	if (x > off_sin || x < -off_sin) {
		hx_2 = hx2 * HALF;
		gtmp = THREE / hx_2; cf_tan2(&gtmp, &U, &V);
		U *= hx_2; U += THREE;
		V *= hx_2; V += THREE;
		gtmp = U * U;
		denum = V * V; denum /= hx_2; denum += gtmp;
		*Rc = (gtmp - (U + V)) / denum; *Rc *= hx2;
		if (Rs) {
			*Rs = (V - gtmp) / denum; *Rs *= x;
		}
	}
	else {
		*Rc = ZERO;
		if (Rs)
		    *Rs = ZERO;
	}
}

static GENERIC
#ifdef __STDC__
Rs0(GENERIC *_x, GENERIC *R)
#else
Rs0(_x, R)		/* Rs0 := -x^3 / 6; *R := sin(x) - x + x^3 / 6 */
GENERIC *_x, *R;		/* _x is input argument */
#endif
{
	register k;
	GENERIC x = *_x, r = ONE, xsq = x*x;

	for (k = ntl_sin - 1; k >= 1; k--) {
		r = -xsq * r; r /= B_sin[k]; r += ONE;
	}
	x *= xsq;				/* x <- x^3 */
	r *= xsq * x;				/* r <- r * (x^5) */
	*R = r / (GENERIC) 120;			/* r <- r * (x^5 / 120) */
	return x / (GENERIC) (-6);			/* return -x^3 / 6 */
}

static GENERIC
#ifdef __STDC__
er_cs(int id, GENERIC *_x, GENERIC *_y)
#else
er_cs(id, _x, _y)
int id;
GENERIC *_x, *_y;		/* input arguments */
#endif
{
	GENERIC x = *_x, y = *_y, xh, xl, hx2h, hx2l;

	x -= theta_h[id];
	xh = x - theta_l[id];
	xl = (x - xh) - theta_l[id];
	hx2l = xh * xl;
#define	gtmp xl
	gtmp = xh * xh * HALF;
	hx2h = gtmp + hx2l;
	hx2l += gtmp - hx2h;
#undef	gtmp

#define	error y
	error = y - sine_h[id];
	error -= (GENERIC) cosine[id] * x;
	error += sine_h[id] * hx2h;
#define	Rs xl
#define	Rc x
	RcRs(&xh, &hx2h, &Rc, &Rs);
	error -= (GENERIC) cosine[id] * Rs;
#undef	Rs
	error -= sine_h[id] * Rc;
	error -= adjust[id];
	error += sine_h[id] * hx2l;
	error += sine_l[id] * (hx2h - Rc);
#undef	Rc
	return error;
#undef	error
}

static GENERIC
#ifdef __STDC__
ersin1(GENERIC *_x, GENERIC *_y)
#else
ersin1(_x, _y)
GENERIC *_x, *_y;		/* input arguments */
#endif
{
	register int neg_x;
	GENERIC x = *_x, y = *_y, error, R;

	neg_x = x < ZERO;
	if (neg_x) {
		x = -x;
		y = -y;
	}
	if (x < (GENERIC) 0.4375) {	/* from this point on x >= ZERO */
		error = y - x;
		if (x > off_sin) {
			error -= Rs0(&x, &R);
			error -= R;
		}
	}
	else
		error = er_cs(x >= (GENERIC) 0.5625, &x, &y);
	return neg_x ? -error : error;
}

static GENERIC
#ifdef __STDC__
ercos1(GENERIC *_x, GENERIC *_y)
#else
ercos1(_x, _y)
GENERIC *_x, *_y;		/* input arguments */
#endif
{
	GENERIC x = *_x, y = *_y, error, Rc, hx2;

	if (x < ZERO)
		x = -x;
	if (x < (GENERIC) 0.3125) {	/* from this point on x >= ZERO */
		hx2 = x * x * HALF;
		RcRs(&x, &hx2, &Rc, (GENERIC *) 0);
		error = y - ONE; error += hx2; error -= Rc;
	}
	else
		error = er_cs(2, &x, &y);
	return error;
}

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

	IORD(rid, 0, x)
	if (rid == 0) {
#define	gtmp scale
		error = ersin1(&x, &y);
#if MACHINE_PI
		APRX_COS(gtmp, x)
		gtmp *= mpi_mu * x;
		error -= gtmp;
#endif
	}
	else {
		register int i;
		GENERIC xi[HML];

		xi[H] = PI2x(rid)[H] - x;
		xi[M] = PI2x(rid)[M];
		xi[L] = PI2x(rid)[L];
		for (i = distd(HML, xi); i <= M; i++)
			xi[i] = ZERO;
		if ((rid - 1) & 2)
			y = -y;
		if (rid & 1) {
			error = ercos1(&xi[H], &y);
			APRX_SIN(gtmp, xi[H])
			error += gtmp * xi[M];
#if MACHINE_PI
			gtmp *= mpi_mu * xi[H];
			error += gtmp;
#endif
		}
		else {
			error = ersin1(&xi[H], &y);
			APRX_COS(gtmp, xi[H])
			error -= gtmp * xi[M];
#if MACHINE_PI
			gtmp *= mpi_mu * xi[H];
			error -= gtmp;
#endif
		}
		if ((rid - 1) & 2) {
			y = -y;
			error = -error;
		}
#undef	gtmp
	}
	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_SIN */

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

	IORD(rid, 1, x)
	if (rid == 0) {
		error = ercos1(&x, &y);
#if MACHINE_PI
		APRX_SIN(gtmp, x)
		gtmp *= mpi_mu * x;
		error += gtmp;
#endif
	}
	else {
		register int i;
		GENERIC xi[HML];

		xi[H] = PI2x(rid)[H] - x;
		xi[M] = PI2x(rid)[M];
		xi[L] = PI2x(rid)[L];
		for (i = distd(HML, xi); i <= M; i++)
			xi[i] = ZERO;
		if (rid & 2)
			y = -y;
		if (rid & 1) {
			error = ersin1(&xi[H], &y);
			APRX_COS(gtmp, xi[H])
			error -= gtmp * xi[M];
#if MACHINE_PI
			gtmp *= mpi_mu * xi[H];
			error -= gtmp;
#endif
		}
		else {
			error = ercos1(&xi[H], &y);
			APRX_SIN(gtmp, xi[H])
			error += gtmp * xi[M];
#if MACHINE_PI
			gtmp *= mpi_mu * xi[H];
			error += gtmp;
#endif
		}
		if (rid & 2) {
			y = -y;
			error = -error;
		}
	}

	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_COS */
#endif	/* TEST_COS || TEST_SIN */
