/* Formerly src/appl/machar.c:
 * void machar()  -- computes ALL `machine constants' at once.
 * -------------  -- compare with ../nmath/i1mach.c & ../nmath/d1mach.c
 *		     which use the C  <float.h> constants !
 *      algorithm 665, collected algorithms from acm.
 *      this work published in transactions on mathematical software,
 *      vol. 14, no. 4, pp. 303-311.
 *
 *  this fortran 77 subroutine is intended to determine the parameters
 *   of the floating-point arithmetic system specified below.  the
 *   determination of the first three uses an extension of an algorithm
 *   due to m. malcolm, cacm 15 (1972), pp. 949-951, incorporating some,
 *   but not all, of the improvements suggested by m. gentleman and s.
 *   marovich, cacm 17 (1974), pp. 276-277.  an earlier version of this
 *   program was published in the book software manual for the
 *   elementary functions by w. j. cody and w. waite, prentice-hall,
 *   englewood cliffs, nj, 1980.
 *
 *  the program as given here must be modified before compiling.  if
 *   a single (double) precision version is desired, change all
 *   occurrences of cs (  ) in columns 1 and 2 to blanks.
 *
 *  parameter values reported are as follows:
 *
 *       ibeta   - the radix for the floating-point representation
 *       it      - the number of base ibeta digits in the floating-point
 *                 significand
 *       irnd    - 0 if floating-point addition chops
 *                 1 if floating-point addition rounds, but not in the
 *                   ieee style
 *                 2 if floating-point addition rounds in the ieee style
 *                 3 if floating-point addition chops, and there is
 *                   partial underflow
 *                 4 if floating-point addition rounds, but not in the
 *                   ieee style, and there is partial underflow
 *                 5 if floating-point addition rounds in the ieee style,
 *                   and there is partial underflow
 *       ngrd    - the number of guard digits for multiplication with
 *                 truncating arithmetic.  it is
 *                 0 if floating-point arithmetic rounds, or if it
 *                   truncates and only  it  base  ibeta digits
 *                   participate in the post-normalization shift of the
 *                   floating-point significand in multiplication;
 *                 1 if floating-point arithmetic truncates and more
 *                   than  it  base  ibeta  digits participate in the
 *                   post-normalization shift of the floating-point
 *                   significand in multiplication.
 *       machep  - the largest negative integer such that
 *                 1.0+float(ibeta)**machep .ne. 1.0, except that
 *                 machep is bounded below by  -(it+3)
 *       negeps  - the largest negative integer such that
 *                 1.0-float(ibeta)**negeps .ne. 1.0, except that
 *                 negeps is bounded below by  -(it+3)
 *       iexp    - the number of bits (decimal places if ibeta = 10)
 *                 reserved for the representation of the exponent
 *                 (including the bias or sign) of a floating-point
 *                 number
 *       minexp  - the largest in magnitude negative integer such that
 *                 float(ibeta)**minexp is positive and normalized
 *       maxexp  - the smallest positive power of  beta  that overflows
 *       eps     - the smallest positive floating-point number such
 *                 that  1.0+eps .ne. 1.0. in particular, if either
 *                 ibeta = 2  or  irnd = 0, eps = float(ibeta)**machep.
 *                 otherwise,  eps = (float(ibeta)**machep)/2
 *       epsneg  - a small positive floating-point number such that
 *                 1.0-epsneg .ne. 1.0. in particular, if ibeta = 2
 *                 or  irnd = 0, epsneg = float(ibeta)**negeps.
 *                 otherwise,  epsneg = (ibeta**negeps)/2.  because
 *                 negeps is bounded below by -(it+3), epsneg may not
 *                 be the smallest number that can alter 1.0 by
 *                 subtraction.
 *       xmin    - the smallest non-vanishing normalized floating-point
 *                 power of the radix, i.e.,  xmin = float(ibeta)**minexp
 *       xmax    - the largest finite floating-point number.  in
 *                 particular  xmax = (1.0-epsneg)*float(ibeta)**maxexp
 *                 note - on some machines  xmax  will be only the
 *                 second, or perhaps third, largest number, being
 *                 too small by 1 or 2 units in the last digit of
 *                 the significand.
 *
 *     latest revision - april 20, 1987
 *
 *     author - w. j. cody
 *              argonne national laboratory
 *
 */

/* #include  from ./platform.c , now as
 a "template" to be used with    DTYPE in {double, long double, ..}
                             MACH_NAME in {machar, machar_LD,   ..}
			         ABS() in {fabs(), fabsl()  ,   ..}
*/
static void
MACH_NAME(int *ibeta, int *it, int *irnd, int *ngrd, int *machep, int *negep,
	  int *iexp, int *minexp, int *maxexp,
	  DTYPE *eps, DTYPE *epsneg, DTYPE *xmin, DTYPE *xmax)
{
	volatile DTYPE a, b, beta, betain, betah, one,
		t, temp, tempa, temp1, two, y, z, zero;
	int i, iz, j, k, mx, nxres;

	one = 1;
	two = one+one;
	zero = one-one;

		/* determine ibeta, beta ala malcolm. */
	a = one; // a = <large> = 9.0072e+15 for 'double' is used later
	do {
		a = a + a;
		temp = a + one;
		temp1 = temp - a;
	}
	while(temp1 - one == zero);
#ifdef _no_longer___did_overflow_ // on IBM PowerPPC ('Power 8')
	int itemp;
	b = one;
	do {
		b = b + b;
		temp = a + b;
		itemp = (int)(temp - a);
	}
	while (itemp == 0);
	*ibeta = itemp;
#else
	*ibeta = (int) FLT_RADIX;
#endif
	beta = *ibeta;

		/* determine it, irnd */

	*it = 0;
	b = one;
	do {
		*it = *it + 1;
		b = b * beta;
		temp = b + one;
		temp1 = temp - b;
	}
	while(temp1 - one == zero);
	*irnd = 0;
	betah = beta / two;
	temp = a + betah;
	if (temp - a != zero)
		*irnd = 1;
	tempa = a + beta;
	temp = tempa + betah;
	if (*irnd == 0 && temp - tempa != zero)
		*irnd = 2;

		/* determine negep, epsneg */

	*negep = *it + 3;
	betain = one / beta;
	a = one;
	for(i=1 ; i<=*negep ; i++)
		a = a * betain;
	b = a;
	for(;;) {
		temp = one - a;
		if (temp - one != zero)
			break;
		a = a * beta;
		*negep = *negep - 1;
	}
	*negep = -*negep;
	*epsneg = a;
	if (*ibeta != 2 && *irnd != 0) {
		a = (a * (one + a)) / two;
		temp = one - a;
		if (temp - one != zero)
			*epsneg = a;
	}

		/* determine machep, eps */

	*machep = -*it - 3;
	a = b;
	for(;;) {
		temp = one + a;
		if (temp - one != zero)
			break;
		a = a * beta;
		*machep = *machep + 1;
	}
	*eps = a;
	temp = tempa + beta * (one + *eps);
	if (*ibeta != 2 && *irnd != 0) {
		a = (a * (one + a)) / two;
		temp = one + a;
		if (temp - one != zero)
			*eps = a;
	}

		/* determine ngrd */

	*ngrd = 0;
	temp = one + *eps;
	if (*irnd == 0 && temp * one - one != zero)
		*ngrd = 1;

	/* determine iexp, minexp, xmin */

	/* loop to determine largest i and k = 2**i such that */
	/*        (1/beta) ** (2**(i)) */
	/* does not underflow. */
	/* exit from loop is signaled by an underflow. */

	i = 0;
	k = 1;
	z = betain;
	t = one + *eps;
	nxres = 0;
	for(;;) {
		y = z;
		z = y * y;

		/* check for underflow here */

		a = z * one;
		temp = z * t;
		if (a+a == zero || ABS(z) >= y)
			break;
		temp1 = temp * betain;
		if (temp1 * beta == z)
			break;
		i = i+1;
		k = k+k;
	}
	if (*ibeta != 10) {
		*iexp = i + 1;
		mx = k + k;
	}
	else {
		/* this segment is for decimal machines only */

		*iexp = 2;
		iz = *ibeta;
		while (k >= iz) {
			iz = iz * *ibeta;
			iexp = iexp + 1;
		}
		mx = iz + iz - 1;
	}
	do {
		/* loop to determine minexp, xmin */
		/* exit from loop is signaled by an underflow */

		*xmin = y;
		y = y * betain;

		/* check for underflow here */

		a = y * one;
		temp = y * t;
		if (a+a == zero || ABS(y) >= *xmin)
			goto L10;
		k = k + 1;
		temp1 = temp * betain;
	}
	while(temp1 * beta != y);
	nxres = 3;
	*xmin = y;
L10:	*minexp = -k;

	/* determine maxexp, xmax */

	if (mx <= k + k - 3 && *ibeta != 10) {
		mx = mx + mx;
		*iexp = *iexp + 1;
	}
	*maxexp = mx + *minexp;

	/* adjust irnd to reflect partial underflow */

	*irnd = *irnd + nxres;

	/* adjust for ieee-style machines */

	if (*irnd == 2 || *irnd == 5)
		*maxexp = *maxexp - 2;

	/* adjust for non-ieee machines with partial underflow */

	if (*irnd == 3 || *irnd == 4)
		*maxexp = *maxexp - *it;

	/* adjust for machines with implicit leading bit in binary */
	/* significand, and machines with radix point at extreme */
	/* right of significand. */

	i = *maxexp + *minexp;
	if (*ibeta == 2 && i == 0)
		*maxexp = *maxexp - 1;
	if (i > 20)
		*maxexp = *maxexp - 1;
	if (a != y)
		*maxexp = *maxexp - 2;
	*xmax = one - *epsneg;
	if (*xmax * one != *xmax)
		*xmax = one - beta * *epsneg;
	*xmax = *xmax / (beta * beta * beta * *xmin);
	i = *maxexp + *minexp + 3;
	if (i>0)
		for(j=1 ; j<=i ; j++) {
			if (*ibeta == 2)
				*xmax = *xmax + *xmax;
			if (*ibeta != 2)
				*xmax = *xmax * beta;
		}
}