#include "ucbtest.h"

/*
 * 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[ ] = "@(#)others.c	1.2 (BeEF) 1/5/95";
#endif

#ifdef __STDC__
#include <stdio.h>
#endif
#include "sys.h"

int sig_bits, sig_decimals;
static GENERIC ulphalf, p2_20, gtiny;

/*
 * gsrand(seed): routine to re-seed the pseudo-random number generator
 * calling sequence:
 *	unsigned seed;
 *	(void) gsrand(seed);
 * Generator is re-started with the new seed.  No value is returned.
 *
 * grand(.): routine to return random GENERICs
 * calling sequence:
 *	GENERIC d, grand(.);
 *	d = grand(.);
 * The value returned ranges from 0.0 (inclusive) to 1.0 (exclusive).
 */
void
#ifdef __STDC__
gsrand(unsigned seed)
#else
gsrand(seed)
unsigned seed;
#endif
{
	srandm(seed);
}

GENERIC
#ifdef __STDC__
grand(void)
#else
grand( )
#endif
{
	long i = randm( );

	return (GENERIC) (i == RANDMAX ? randm( ) : i) / (GENERIC) RANDMAX;
}

GENERIC
#ifdef __STDC__
ld_2_n(int n)
#else
ld_2_n(n)
int n;
#endif
{
	GENERIC x, y;

	if (!n)
		return ONE;
	else if (n > 0)
		x = TWO;
	else {
		x = HALF;
		n = -n;
	}
	while (!(n & 1)) {
		x *= x;
		n >>= 1;
	}
	for (y = x; --n > 0; y *= x)
		while (!(n & 1)) {
			x *= x;
			n >>= 1;
		}
	return y;		/* returns 2^n */
}

GENERIC
#ifdef __STDC__
getulp(GENERIC *_y, GENERIC *t)
#else
getulp(_y, t)
GENERIC *_y, *t;			/* _y is input argument */
#endif
{
	GENERIC y = *_y, ulpy, tmpy;

	*t = ONE;			/* initialize scaling factor */
	if (y == ZERO) {		/* if 0 then ulp(y) is undefined */
		ulpy = gtiny;			/* give it a tiny value... */
		*t /= ulphalf;			/* and a huge scaling factor */
	}
	else {
		if (y < ZERO)		/* force a positive y */
			y = -y;
		tmpy = y * ulphalf;
		if (tmpy == ZERO) {	/* if tmpy underflows... */
			tmpy = y;		/* use y as tmpy */
			y /= ulphalf;		/* scale y up */
			*t /= ulphalf;		/* set scaling factor */
		}
		ulpy = y + tmpy;
		ulpy -= y;
		if (ulpy == ZERO) {	/* IEEE-style round-to-nearst-even */
			tmpy += tmpy;
			ulpy = y + tmpy;
			ulpy -= y;
		}
	}
	return ulpy;
}

/*
 * This routine takes whatever stored in xs[1] and returns the 2 nearby
 * numbers in
 * xs[0] := xs[1] - 1ULP_of_xs[1]
 * xs[2] := xs[1] + 1ULP_of_xs[1]
 */
void
#ifdef __STDC__
nearby(GENERIC *xs)
#else
nearby(xs)
GENERIC *xs;
#endif
{
	GENERIC x, tmpx, gtmp, t;

	x = xs[1];
	if (x == ZERO) {
		xs[0] = -gtiny;
		xs[2] = gtiny;
	}
	else {
		t = ulphalf;		/* scaling factor */
		if (x < ZERO)		/* force same sign on t and x */
			t = -t;
		tmpx = x * t;		/* tmpx >= 0 */
		if (tmpx == ZERO) {	/* if tmpx underflows... */
			tmpx = x;		/* use x as tmpx */
			x /= t;			/* scale x up */
		}
		else			/* otherwise... */
			t = ZERO;		/* flag normal */
		xs[0] = x - tmpx;
		if (xs[0] == x) {	/* IEEE-style round-to-nearest-even */
			gtmp = tmpx + tmpx;
			xs[0] = x - gtmp;
		}
		xs[2] = x + tmpx;
		if (xs[2] == x) {	/* IEEE-style round-to-nearest-even */
			gtmp = tmpx + tmpx;
			xs[2] = x + gtmp;
		}
		if (t != ZERO) {	/* scale back down if needed */
			xs[0] *= t;
			xs[2] *= t;
		}
	}
}

void
#ifdef __STDC__
unpack(int pcs_needed, int pcs, int n, LONG *yp, GENERIC *yl, GENERIC *yh)
#else
unpack(pcs_needed, pcs, n, yp, yl, yh)
LONG *yp;
GENERIC *yl, *yh;
#endif
{
	int i, j, k;
	GENERIC t, expnt;

	for (i = 0; i <= n; i++) {
		t = ONE;
		k = pcs * i + 1;
		expnt = (GENERIC) yp[k];
		if ((int) yp[k] < 0)
			expnt = -ONE / expnt;
		for (yh[i] = ZERO, j = k + pcs_needed; j >= k + 1; j--) {
			yh[i] += (GENERIC) yp[j];
			yh[i] /= p2_20;
		}
		for (yl[i] = yh[i], j = k + 1; j <= k + pcs_needed; j++) {
			yl[i] *= p2_20;
			yl[i] -= (GENERIC) yp[j];
			t /= p2_20;
		}
		yl[i] *= -t;
		yl[i] *= expnt;
		yh[i] *= expnt;
		if ((int) yp[k - 1] < 0) {
			yl[i] = -yl[i];
			yh[i] = -yh[i];
		}
	}
}

void
#ifdef __STDC__
update(int n, int m, int is, GENERIC *_x, GENERIC *_error, struct summary *s)
#else
update(n, m, is, _x, _error, s)
int n, m, is;
GENERIC *_x, *_error;		/* input arguments */
struct summary *s;
#endif
{
	if (*_error > s->ep) {
		s->xp = *_x;
		s->np = n;
		s->ep = *_error;
	} else if (*_error < s->en) {
		s->xn = *_x;
		s->nn = n;
		s->en = *_error;
	}
	if (m && !s->mcnt) {
		s->xm = *_x;
		s->nm = n;
	}
	s->mcnt += m;
	s->scnt += is;				/* for sin, cos and atan */
}

void
#ifdef __STDC__
double_check(int i, char *name, struct rational *rat, FLOAT *x)
#else
double_check(i, name, rat, x)
int i;
char *name;
struct rational *rat;
FLOAT *x;
#endif
{
	FLOAT z;

	if (!rat->q)
		return;
	z = (FLOAT) rat->p / (FLOAT) rat->q;
	if (*x == z)
		return;
	(void) printf(
		"-->COMPILER FLAW: %s[%d] (%d/%d) was not converted exactly!\n",
		name, i, rat->p, rat->q);
	*x = z;
}

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

	if (significand_length == 0) significand_length=get_significand();
	sig_bits = significand_length;
	for (i = 2; i <= sig_bits; i++)
		if ((gtiny = ld_2_n(i - BITS_2M)) > ZERO)
			break;
	if (gtiny == ZERO) {
		(void) printf("--> 2**n == 0 for n := %d through %d.\n",
			2 - BITS_2M, sig_bits - BITS_2M);
		(void) printf("--> Is BITS_M (%d) defined correctly?\n",
			BITS_M);
		gtiny = (GENERIC) 2.938736e-39;
	}
	ulphalf = ld_2_n(-sig_bits);
	p2_20 = ld_2_n(20);
	sig_decimals = (FLOAT) sig_bits * 0.30103 + 0.5;
	if (sig_decimals <= 17)
		sig_decimals = 17;
	else if (sig_decimals >= 46)
		sig_decimals = 46;
}
