/*
 * 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_log.c	1.1 (BeEF) 1/5/95";
#endif
/*
 * er_log(.) computes the error committed by log(x) to within 0.0625 ULPs.
 * Covers [2^-16.5, 2^16.5).  (Bound proved: 0.052 ULPs).
 * extern GENERIC ld_2_n(.), cf_atn(.), getulp(.);
 * extern void unpack(.), update(.), double_check(.);
 */
#include "sys.h"
extern int sig_bits;
#if BEEF_QA
extern GENERIC h_diff( );
#endif

#define	MID_LOG	6
#define	MAXPCS	10			/* should be >= 6 */
#define	MASKOFF	5
#define	MAXEXPO	(1 << (MASKOFF - 1))	/* = 2^(MASKOFF - 1) */

GENERIC off_log;			/* ALSO USED in: ac_l1p(.) */

static GENERIC p2_max, ln2_hi, ln2_lo;

LONG ln2_p[] = {			/* ALSO USED in: bgn_exp(.) */
  1,  -2, 0x162e42, 0xfefa3, 0x9ef35, 0x793c7, 0x67300,		/* log(2) */
	   0x7e5ed, 0x5e81e, 0x6864c, 0xe5316, 0xc5b14,
	   0x1a2eb, 0x71755, 0xf457c, 0xf70ec, 0x40dbd,
};

static FLOAT x_log[] = {
     7.246093750e-01,			/*  0 */
     7.695312500e-01,			/*  1 */
     7.929687500e-01,			/*  2 */
     8.164062500e-01,			/*  3 */
     8.398437500e-01,			/*  4 */
     8.632812500e-01,			/*  5 */
     1.000000000e+00,			/*  6 */
     1.132812500e+00,			/*  7 */
     1.158203125e+00,			/*  8 */
     1.193359375e+00,			/*  9 */
     1.228515625e+00,			/* 10 */
     1.265625000e+00,			/* 11 */
     1.324218750e+00,			/* 12 */
     1.408203125e+00,			/* 13 */
};
#define CTR_LOG (sizeof(x_log) / sizeof(x_log[0]))

static struct rational xpq_log[ ] = {
	{  742, 1024},		/*  0:  7.2460937500e-01 */
	{  788, 1024},		/*  1:  7.6953125000e-01 */
	{  812, 1024},		/*  2:  7.9296875000e-01 */
	{  836, 1024},		/*  3:  8.1640625000e-01 */
	{  860, 1024},		/*  4:  8.3984375000e-01 */
	{  884, 1024},		/*  5:  8.6328125000e-01 */
	{ 1024, 1024},		/*  6:  1.0000000000e+00 */
	{ 1160, 1024},		/*  7:  1.1328125000e+00 */
	{ 1186, 1024},		/*  8:  1.1582031250e+00 */
	{ 1222, 1024},		/*  9:  1.1933593750e+00 */
	{ 1258, 1024},		/* 10:  1.2285156250e+00 */
	{ 1296, 1024},		/* 11:  1.2656250000e+00 */
	{ 1356, 1024},		/* 12:  1.3242187500e+00 */
	{ 1442, 1024},		/* 13:  1.4082031250e+00 */
};

GENERIC y_ln_l[CTR_LOG], y_ln_h[CTR_LOG];	/* ALSO USED in: er_l1p(.) */

static LONG y_ln_p[] = {
    -1,  -4, 0x149da7, 0xf3bcc, 0x41ecc, 0xd36bd, 0x2e66a,	/*  0 */
              0x6c718, 0x21b02, 0xec7a5, 0x1b6b8, 0x0735d,
    -1,  -4, 0x10c42d, 0x67616, 0x2e311, 0x62c79, 0xd5d11,	/*  1 */
              0xee41e, 0x3b351, 0xff419, 0x49216, 0xca302,
    -1,  -8, 0x1db13d, 0xb0d48, 0x94035, 0x423a9, 0x3f2d9,	/*  2 */
              0x71062, 0xf5613, 0x9580f, 0xd566f, 0x151cc,
    -1,  -8, 0x19f6c4, 0x07089, 0x66413, 0x5a196, 0x05e67,	/*  3 */
              0xef382, 0xd7c64, 0xd5883, 0x4b04b, 0x5f89b,
    -1,  -8, 0x16574e, 0xbe8c1, 0x339f1, 0x65878, 0x5cef2,	/*  4 */
              0x095f4, 0xf00ef, 0xf4801, 0xcfcd1, 0x34661,
    -1,  -8, 0x12d161, 0x0c868, 0x139d6, 0xccb81, 0xb4a0d,	/*  5 */
              0x41109, 0x0848d, 0x6f582, 0xf0e24, 0x72971,
     1,   1, 0x000000, 0x00000, 0x00000, 0x00000, 0x00000,	/*  6 */
              0x00000, 0x00000, 0x00000, 0x00000, 0x00000,
     1, -16, 0x1fec91, 0x31dbe, 0xabaaa, 0x2e519, 0x9f932,	/*  7 */
              0x4e3bf, 0xe91e2, 0xba812, 0x02ec6, 0x15272,
     1,  -8, 0x12cca0, 0xf5f5f, 0x25087, 0x37280, 0x7703f,	/*  8 */
              0xa7911, 0xbc279, 0x27e20, 0x0e1e1, 0x557ac,
     1,  -8, 0x16a079, 0xd0f7a, 0xad1fc, 0x22468, 0xa7ab0,	/*  9 */
              0x1d0d2, 0x2f4e8, 0x2ae81, 0x8a21a, 0x51d83,
     1,  -8, 0x1a57df, 0x28244, 0xdcce4, 0x650ec, 0xd5db1,	/* 10 */
              0xc724d, 0x16788, 0x2a217, 0x9ef90, 0xa1d3c,
     1,  -8, 0x1e2707, 0x6e2af, 0x2e5e9, 0xea87f, 0xfe1fe,	/* 11 */
              0x9e155, 0xdb94e, 0xbc401, 0x7f6f9, 0x57dd0,
     1,  -4, 0x11f8ff, 0x9e48a, 0x2f28d, 0x80819, 0x7ced3,	/* 12 */
              0xe58cf, 0x23e43, 0x622b0, 0xb6ee3, 0x7d610,
     1,  -4, 0x15e87b, 0x20c29, 0x549f4, 0x63ddc, 0xe3e81,	/* 13 */
              0xd7ac0, 0xf4aba, 0x8be5b, 0x934db, 0xa4ae6,
};

int PCS_LOG =					/* ALSO USED in er_em1.c */
	((sizeof(y_ln_p) / sizeof(y_ln_p[0])) / CTR_LOG);

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

	for (i = 0; i < CTR_LOG; i++)
		double_check(i, "x_log", &xpq_log[i], &x_log[i]);
}

/* NOTE that in order to run this routine, A_atn[], B_atn[] & trm_atn
 *      must be initialized (currently they are initialized in bgn_atn(.)
 *      and this routine calls cf_atn(.) in computing accurate log(.).
 */
void
#ifdef __STDC__
bgn_log(void)
#else
bgn_log( )
#endif
{
	int pieces;
	GENERIC gtmp, chopped;

	off_log = ld_2_n(-sig_bits / 2 - 2);
	p2_max = ld_2_n(MAXEXPO);

	pieces = (sig_bits + 15) / 20 + 1;
	unpack(pieces, (int) PCS_LOG, (int) (CTR_LOG - 1),
		y_ln_p, y_ln_l, y_ln_h);

	pieces = sig_bits / 10 + 1;
	unpack(pieces, (int) (PCS_LOG + 5), (int) 0,
		ln2_p, &ln2_lo, &ln2_hi);

	gtmp = ln2_hi * ld_2_n(MASKOFF + 1);
	chopped = gtmp + ln2_hi;
	chopped = chopped - gtmp;
	ln2_lo += ln2_hi - chopped;
	ln2_hi = chopped;
}

#if TEST_LOG
static int
#ifdef __STDC__
ac_log(int i, GENERIC *_x, GENERIC *y)
#else
ac_log(i, _x, y)
int i;
GENERIC *_x, *y;		/* _x is input argument */
#endif
{
	GENERIC x = *_x, xm1, xm1sq, rho;

	xm1 = (x - (GENERIC) x_log[i]) / (GENERIC) x_log[i];
	x /= (GENERIC) x_log[i];	/* MUST come *after* xm1 is computed */
	xm1sq = xm1 * xm1;
	rho = xm1 - xm1sq / (x + ONE);

	y[0] = y_ln_h[i];
	y[1] = xm1;
	y[2] = -xm1sq * HALF;
	y[3] = -xm1 * y[2]; y[3] /= x + ONE;
	if (rho > off_log || rho < -off_log) {
		xm1 = (GENERIC) (-12.0) / (rho * rho);
		y[4] = -rho / cf_atn(&xm1);
	}
	else
		y[4] = ZERO;
	y[5] = y_ln_l[i];
	return 6;			/* 6 should be less than MAXPCS */
}

void
#ifdef __STDC__
er_log(int fid, int rid, GENERIC *xs_frz, struct summary *s)
#else
er_log(fid, rid, xs_frz, s)
int fid, rid;
GENERIC *xs_frz;
struct summary *s;
#endif
{
	register int n, i, j;
	int terms[3];
	GENERIC y, prevy = ZERO, xs[3], ac_y[MAXPCS * 3];
	GENERIC error, ulp, scale;

#ifdef lint
	fid = fid;
#endif	/* lint */
	for (i = 0; i < 3; i++) {
		terms[i] = ac_log(rid, &xs_frz[i], ac_y + i * MAXPCS);
		xs[i] = xs_frz[i] / p2_max;
	}
	for (n = -MAXEXPO; n <= MAXEXPO; n++) {
		for (i = 0; i < 3; i++) {
			FUNC_LOG(y, xs[i])
			error = y - (GENERIC) n * ln2_hi;
			for (j = i * MAXPCS; j < i * MAXPCS + terms[i]; j++)
				error -= ac_y[j];
			error -= (GENERIC) n * ln2_lo;

			ulp = y - error; ulp = getulp(&ulp, &scale);
			error *= scale;
#if GUARD_AGAINST_OVERFLOW
			if (GABS(error) / GHUGE > ulp)    /* we gave up */
				error = ZERO;
#endif
			error /= ulp;
#if BEEF_QA
			error = h_diff(fid, &xs[i], &y, &error);
#endif
			j = i && y < prevy;
			update(n, j, 0, &xs_frz[i], &error, s);
			prevy = y;
			xs[i] += xs[i];
		}
	}
}
#endif	/* TEST_LOG */
