/* NUMmachar.c */
/* David Weenink, 20010719 */

#include "NUMmachar.h"
#include "melder.h"

#ifndef _NUM_h_
	#include "NUM.h"
#endif
#ifdef _WIN32
	#pragma global_optimizer off
	#pragma optimization_level 1
#endif

#define SEPARATE  (void) Melder_hasError ();
static struct machar_Table machar_table;
machar_Table NUMfpp = NULL;

static void NUMmachar_d (void)
{
	volatile int i, itemp, iz, j, k, mx, nxres;
	volatile int base, emin, emax, iexp, it, irnd, machep, ngrd, negep;
	volatile double a, b, beta, betah, betain, eps, epsneg, one, rmin, rmax;
	volatile double t, temp, temp1, temp2, tempa, two, y, z, zero;
	
	one = 1;
	two = one + one;
	zero = one - one;
	a = one;
	
	/*
		Determine base and beta by the method of M. Malcolm.
	*/
	
	do
	{
		a += a;
		SEPARATE
		temp = a + one;
		SEPARATE
		temp1 = temp - a;
		SEPARATE
		temp2 = temp1 - one;
		SEPARATE
	} while (temp2 == zero);
	
	b = one;
	do
	{
		b += b;
		SEPARATE
		temp = a + b;
		SEPARATE
		itemp = (int) (temp - a);
		SEPARATE
	} while (itemp == 0);
	base = itemp;
	
	/*
		Determine it and irnd.
	*/
	
	beta = base;
	it = 0;
	b = one;
	do
	{
		++it;
		SEPARATE
		b *= beta;
		SEPARATE
		temp = b + one;
		SEPARATE
		temp1 = temp - b;
		SEPARATE
		temp2 = temp1 - one;
		SEPARATE
	} while (temp2 == zero);
	irnd = 0;
	betah = beta / two;
	temp = a + betah;
	SEPARATE
	temp2 = temp - a;
	SEPARATE
	if (temp2 != zero) irnd = 1;
	tempa = a + beta;
	SEPARATE
	temp = tempa + betah;
	SEPARATE
	temp2 = temp - tempa;
	SEPARATE
	if (irnd == 0 && temp2 != zero) irnd = 2;
	
	/*
		Determine negep and epsneg.
	*/
	
	negep = it + 3;
	betain = one / beta;
	a = one;
	for (i=1; i <= negep; i++) a *= betain;
	b = a;
	for (;;)
	{
		temp = one - a;
		SEPARATE
		temp2 = temp - one;
		SEPARATE
		if (temp2 != zero) break;
		a *= beta;
		--negep;
	}
	negep = -negep;
	epsneg = a;
	
	/*
		Determine machep and eps.
	*/
	
	machep = -it - 3;
	a = b;
	for (;;)
	{
		temp = one + a;
		SEPARATE
		temp2 = temp - one;
		SEPARATE
		if (temp2 != zero) break;
		a *= beta;
		++machep;
	}
	eps = a;
	
	/*
		Determine ngrd.
	*/
	
	ngrd = 0;
	temp = one + eps;
	SEPARATE
	temp1 = temp * one;
	SEPARATE
	temp2 = temp1 - one;
	SEPARATE
	if (irnd == 0 && temp2 != zero) ngrd = 1;
	
	/*
		Determine iexp.
	*/
	
	i = 0;
	k = 1;
	z = betain;
	t = one + eps;
	nxres = 0;
	for (;;)
	{
		/*
			Loop until an underflow occurs, then exit.
		*/
		y = z;
		SEPARATE
		z = y * y;
		SEPARATE
		a = z * one;
		temp = z * t;
		SEPARATE
		temp1 = a + a;
		SEPARATE
		temp2 = fabs (z);
		SEPARATE
		if (temp1 == zero || temp2 >= y) break;
		temp1 = temp * betain;
		SEPARATE
		temp2 = temp1 * beta;
		SEPARATE
		if (temp2 == z) break;
		++i;
		k += k;
	}
	
	if (base != 10)
	{
		iexp = i + 1;
		mx = k + k;
	}
	else /* For decimal machines only. */
	{
		iexp = 2;
		iz = (base);
		while (k >= iz)
		{
			iz *= base;
			++iexp;
		}
		mx = iz + iz - 1;
	}
	
	/*
		Determine emin and xmin.
		Loop until an underflow occurs.
	*/
	
	for (;;)
	{
		rmin = y;
		y *= betain;
		SEPARATE
		a = y * one;
		temp = y * t;
		SEPARATE
		temp1 = a + a;
		temp2 = fabs (y);
		SEPARATE
		if (temp1 != zero && temp2 < rmin)
		{
			++k;
			temp1 = temp * betain;
			SEPARATE
			temp2 = temp1 * beta;
			SEPARATE
			if (temp2 == y && temp != y)
			{
				nxres = 3;
				rmin = y;
				break;
			}
		}
		else break;
	}
	
	/*
		Determine emax, rmax.
	*/
	
	emin = -k;
	if (mx <= k + k - 3 && base != 10)
	{
		mx += mx;
		++iexp;
	}
	emax = mx + emin;
	
	/*
		Adjust irnd to reflect partial underflow.
	*/
	
	irnd += nxres;
	
	/*
		Adjust for IEEE-style machines.
	*/
	
	if (irnd >= 2) emax -= 2;
	i = emax + emin;
	
	/*
		Adjust for machines with implicit leading bit in mantissa, and
		machines with radixpoint at extreme right of mantissa.
	*/
	
	if (base == 2 && !i) --emax;
	if (i > 20) --emax;
	if (a != y) emax -= 2;
	rmax = one - epsneg;
	SEPARATE
	temp1 = rmax * one;
	SEPARATE
	if (temp1 != rmax) {
		temp1 = beta * epsneg;
		SEPARATE
		rmax = one - temp1;
	}
	temp1 = rmin * beta;
	SEPARATE
	temp1 *= beta;
	SEPARATE
	temp1 *= beta;
	SEPARATE
	rmax /= temp1;
	i = emax + emin + 3;
	for (j=1; j <= i;j++)
	{
		if (base == 2) rmax += rmax;
		else rmax *= beta;
	}
	
	NUMfpp = & machar_table;
	
	NUMfpp -> base_d = base; 
	NUMfpp -> it_d = it; 
	NUMfpp -> irnd_d = irnd;
	NUMfpp -> ngrd_d = ngrd; 
	NUMfpp -> machep_d = machep; 
	NUMfpp -> negep_d = negep;
	NUMfpp -> iexp_d = iexp; 
	NUMfpp -> emin_d = emin; 
	NUMfpp -> emax_d = emax;
	NUMfpp -> eps_d = eps; 
	NUMfpp -> epsneg_d = epsneg; 
	NUMfpp -> rmin_d = rmin;
	NUMfpp -> rmax_d = rmax;
}

static void NUMmachar_f (void)
{
	volatile int i, itemp, iz, j, k, mx, nxres;
	volatile int base, emin, emax, iexp, it, irnd, machep, ngrd, negep;
	volatile float a, b, beta, betah, betain, eps, epsneg, one, rmin, rmax;
	volatile float t, temp, temp1, temp2, tempa, two, y, z, zero;

	one = 1;
	two = one + one;
	zero = one - one;
	a = one;
	
	/*
		Determine base and beta by the method of M. Malcolm.
	*/
	
	do
	{
		a += a;
		SEPARATE
		temp = a + one;
		SEPARATE
		temp1 = temp - a;
		SEPARATE
		temp2 = temp1 - one;
		SEPARATE
	} while (temp2 == zero);
	
	b = one;
	do
	{
		b += b;
		SEPARATE
		temp = a + b;
		SEPARATE
		itemp = (int) (temp - a);
		SEPARATE
	} while (itemp == 0);
	base = itemp;
	
	/*
		Determine it and irnd.
	*/
	
	beta = base;
	it = 0;
	b = one;
	do
	{
		++it;
		SEPARATE
		b *= beta;
		SEPARATE
		temp = b + one;
		SEPARATE
		temp1 = temp - b;
		SEPARATE
		temp2 = temp1 - one;
		SEPARATE
	} while (temp2 == zero);
	irnd = 0;
	betah = beta / two;
	temp = a + betah;
	SEPARATE
	temp2 = temp - a;
	SEPARATE
	if (temp2 != zero) irnd = 1;
	tempa = a + beta;
	SEPARATE
	temp = tempa + betah;
	SEPARATE
	temp2 = temp - tempa;
	SEPARATE
	if (irnd == 0 && temp2 != zero) irnd = 2;
	
	/*
		Determine negep and epsneg.
	*/
	
	negep = it + 3;
	betain = one / beta;
	a = one;
	for (i=1; i <= negep; i++) a *= betain;
	b = a;
	for (;;)
	{
		temp = one - a;
		SEPARATE
		temp2 = temp - one;
		SEPARATE
		if (temp2 != zero) break;
		a *= beta;
		--negep;
	}
	negep = -negep;
	epsneg = a;
	
	/*
		Determine machep and eps.
	*/
	
	machep = -it - 3;
	a = b;
	for (;;)
	{
		temp = one + a;
		SEPARATE
		temp2 = temp - one;
		SEPARATE
		if (temp2 != zero) break;
		a *= beta;
		++machep;
	}
	eps = a;
	
	/*
		Determine ngrd.
	*/
	
	ngrd = 0;
	temp = one + eps;
	SEPARATE
	temp1 = temp * one;
	SEPARATE
	temp2 = temp1 - one;
	SEPARATE
	if (irnd == 0 && temp2 != zero) ngrd = 1;
	
	/*
		Determine iexp.
	*/
	
	i = 0;
	k = 1;
	z = betain;
	t = one + eps;
	nxres = 0;
	for (;;)
	{
		/*
			Loop until an underflow occurs, then exit.
		*/
		y = z;
		SEPARATE
		z = y * y;
		SEPARATE
		a = z * one;
		temp = z * t;
		SEPARATE
		temp1 = a + a;
		SEPARATE
		temp2 = fabs (z);
		SEPARATE
		if (temp1 == zero || temp2 >= y) break;
		temp1 = temp * betain;
		SEPARATE
		temp2 = temp1 * beta;
		SEPARATE
		if (temp2 == z) break;
		++i;
		k += k;
	}
	
	if (base != 10)
	{
		iexp = i + 1;
		mx = k + k;
	}
	else /* For decimal machines only. */
	{
		iexp = 2;
		iz = (base);
		while (k >= iz)
		{
			iz *= base;
			++iexp;
		}
		mx = iz + iz - 1;
	}
	
	/*
		Determine emin and xmin.
		Loop until an underflow occurs.
	*/
	
	for (;;)
	{
		rmin = y;
		y *= betain;
		SEPARATE
		a = y * one;
		temp = y * t;
		SEPARATE
		temp1 = a + a;
		temp2 = fabs (y);
		SEPARATE
		if (temp1 != zero && temp2 < rmin)
		{
			++k;
			temp1 = temp * betain;
			SEPARATE
			temp2 = temp1 * beta;
			SEPARATE
			if (temp2 == y && temp != y)
			{
				nxres = 3;
				rmin = y;
				break;
			}
		}
		else break;
	}
	
	/*
		Determine emax, rmax.
	*/
	
	emin = -k;
	if (mx <= k + k - 3 && base != 10)
	{
		mx += mx;
		++iexp;
	}
	emax = mx + emin;
	
	/*
		Adjust irnd to reflect partial underflow.
	*/
	
	irnd += nxres;
	
	/*
		Adjust for IEEE-style machines.
	*/
	
	if (irnd >= 2) emax -= 2;
	i = emax + emin;
	
	/*
		Adjust for machines with implicit leading bit in mantissa, and
		machines with radixpoint at extreme right of mantissa.
	*/
	
	if (base == 2 && !i) --emax;
	if (i > 20) --emax;
	if (a != y) emax -= 2;
	rmax = one - epsneg;
	SEPARATE
	temp1 = rmax * one;
	SEPARATE
	if (temp1 != rmax) {
		temp1 = beta * epsneg;
		SEPARATE
		rmax = one - temp1;
	}
	temp1 = rmin * beta;
	SEPARATE
	temp1 *= beta;
	SEPARATE
	temp1 *= beta;
	SEPARATE
	rmax /= temp1;
	i = emax + emin + 3;
	for (j=1; j <= i;j++)
	{
		if (base == 2) rmax += rmax;
		else rmax *= beta;
	}
	
	NUMfpp = & machar_table;
	
	NUMfpp -> base_f = base; NUMfpp -> it_f = it;
	NUMfpp -> irnd_f = irnd; NUMfpp -> ngrd_f = ngrd; 
	NUMfpp -> machep_f = machep; NUMfpp -> negep_f = negep;
	NUMfpp -> iexp_f = iexp; NUMfpp -> emin_f = emin;
	NUMfpp -> emax_f = emax; NUMfpp -> eps_f = eps;
	NUMfpp -> epsneg_f = epsneg; NUMfpp -> rmin_f = rmin;
	NUMfpp -> rmax_f = rmax;
}

void NUMmachar (void)
{
	if (NUMfpp) return;
	NUMmachar_d ();
	NUMmachar_f ();
}

/* End of file NUMmachar.c */
