#include <math.h>
#include <stdio.h>

#define ZERO 0.0
#define ABS(xxx) ((xxx>=ZERO)?(xxx):(-xxx))

/* do the hexadecimal output according to whether float or double used */
#ifdef SINGLE
	union intorreal{
		int n;
		float x;
	} outbits;
#define HEXOUT(xxx) outbits.x = xxx; printf("%8x\n",outbits.n);
#endif
#ifdef DOUBLE
	union intorreal{
		int n[2];
		double x;
	} outbits;
#define HEXOUT(xxx) outbits.x = xxx; printf("%8x %8x\n",outbits.n[0], outbits.n[1]);
#endif

main()
{
	/*
	c     **********
	c
	c     this program prints hardware-determined
	c     machine constants obtained by smchar, a subroutine due to
	c     w. j. cody.
	c
	c     descriptions of the machine constants are
	c     given in the prologue comments of smchar.
	c
	c     subprograms called
	c
	c       smchar
	c
	c     minpack. version of january 1979.
	c     burton s. garbow, kenneth e. hillstrom, jorge j. more
	c
	c     c version Fri Oct 25 22:54:26 BST 1985
	c       Tim Hopkins
	c
	c     **********
*/
	int ibeta,iexp,irnd,it,machep,maxexp,minexp,negep,ngrd;
#ifdef SINGLE
	float
#endif
#ifdef DOUBLE
	double
#endif
	    dwarf,eps,epsmch,epsneg,giant,xmax,xmin;
	machar(&ibeta,&it,&irnd,&ngrd,&machep,&negep,
	&iexp,&minexp,&maxexp,&eps,&epsneg,&xmin,&xmax);
	printf("smchar constants\n");
	printf("ibeta = %d\n",ibeta);
	printf("it    = %d\n",it);
	printf("irnd  = %d\n",irnd);
	printf("ngrd  = %d\n",ngrd);
	printf("machep = %d\n",machep);
	printf("negep = %d\n",negep);
	printf("iexp = %d\n",iexp);
	printf("minexp = %d\n",minexp);
	printf("maxexp = %d\n",maxexp);
	printf("eps = "); HEXOUT(eps);
	printf("epsneg = "); HEXOUT(epsneg);
	printf("xmin = "); HEXOUT(xmin);
	printf("xmax = "); HEXOUT(xmax);
}

machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,eps,epsneg,xmin,xmax)
int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd;
#ifdef SINGLE
	float
#endif
#ifdef DOUBLE
	double
#endif
*eps,*epsneg,*xmax,*xmin;
{
	int i, iz, j, k, mx;
#ifdef SINGLE
	float
#endif
#ifdef DOUBLE
	double
#endif
	    a, b, beta, betain, betam1, one, temp,  y, z, zero;
	/*
	      this subroutine is intended to determine the characteristics
	      of the floating-point arithmetic system that are specified
	      below.  the first three are determined according to 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.  the version given here is for single precision.
	      cards containing  cd  in columns 1 and 2 can be used to
	      convert the subroutine to double precision by replacing
	      existing cards in the obvious manner.
	 
	 
	        *ibeta   - the radix of 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
	        *ngrd    - the number of guard digits for multiplication.  it is
	                  0 if  irnd=1, or if  irnd=0  and only  it  base  ibeta
	                    digits participate in the post normalization shift
	                    of the floating-point significand in multiplication
	                  1 if  irnd=0  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+real(ibeta)***machep .ne. 1.0, except that
	                  *machep is bounded below by  -(it+3)
	        *negeps  - the largest negative integer such that
	                  1.0-real(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
	                  real(ibeta) ** *minexp is a positive floating-point
	                  number
	        *maxexp  - the largest positive integer exponent for a finite
	                  floating-point number
	        *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 = real(ibeta)***machep.
	                  otherwise,  eps = (real(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 = real(ibeta) ** *negeps.
	                  otherwise,  *epsneg = (ibeta ** *negeps)/2.  because
	                  *negeps is bounded below by -(it+3), *epsneg may not
	                  be the smallest number which can alter 1.0 by
	                  subtraction.
	        *xmin    - the smallest non-vanishing floating-point power of the
	                  radix.  in particular,  *xmin = real(ibeta) ** *minexp
	        *xmax    - the largest finite floating-point number.  in
	                  particular   *xmax = (1.0-*epsneg)*real(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 - october 22, 1979
	 
	      author - w. j. cody
	               argonne national laboratory
	 
	 -----------------------------------------------------------------
	 */

#ifdef SINGLE
	one = (float)1;
	zero = (float)0;
#endif
#ifdef DOUBLE
	one = (double)1;
	zero = (double)0;
#endif
	/*
	-----------------------------------------------------------------
	     determine ibeta,beta ala malcolm
	-----------------------------------------------------------------
	*/
	a = one;
	do {
		a = a+a;
		temp = a + one;
		temp = temp -a;
		temp = temp - one;
	} 
	while(temp == zero);
	b = one;
	do{
		b = b+b;
		temp = a+b;
		temp = temp - a;
	} 
	while (temp == zero);
	*ibeta = (int)((a+b)-a);
#ifdef SINGLE
	beta = (float)(*ibeta);
#endif
#ifdef DOUBLE
	beta = (double)(*ibeta);
#endif
	/*
	-----------------------------------------------------------------
	     determine it, irnd
	-----------------------------------------------------------------
	*/
	*it = 0;
	b = one;
	do {
		(*it)++;
		b = b*beta;
		temp = b + one;
		temp = temp - b;
		temp = temp - one;
	} 
	while(temp == zero);
	*irnd = 0;
	betam1 = beta-one;
	temp = a + betam1;
	temp = temp - a;
	if (temp != zero)
		*irnd = 1;
	/*
	-----------------------------------------------------------------
	     determine *negep, *epsneg
	-----------------------------------------------------------------
	*/
	*negep = *it+3;
	betain = one/beta;
	a = one;

	for (i=1; i<=*negep; i++)
		a = a*betain;

	b = a;
	temp = one - a;
	temp = temp - one;
	while (temp == zero) {
		a = a*beta;
		*negep = *negep-1;
		temp = one - a;
		temp = temp - one;
	}
	*negep = -*negep;
	*epsneg = a;
	if (*ibeta!=2&&*irnd!=0) {
		a = (a*(one+a))/(one+one);
		temp = one - a;
		temp = temp - one;
		if (temp!=zero)
			*epsneg = a;
	}
	/*
	-----------------------------------------------------------------
	     determine *machep, eps
	-----------------------------------------------------------------
	*/
	*machep = -*it-3;
	a = b;
	temp = one + a;
	temp = temp - one;
	while (temp==zero) {
		a *= beta;
		(*machep)++;
		temp = one + a;
		temp = temp - one;
	}
	*eps = a;
	if (*ibeta!=2&&*irnd!=0) {
		a = (a*(one+a))/(one+one);
		temp = one + a;
		temp = temp - one;
		if (temp!=zero)
			*eps = a;
	}
	/*
	-----------------------------------------------------------------
	     determine *ngrd
	-----------------------------------------------------------------
	*/
	*ngrd = 0;
	temp = one + *eps;
	temp = temp * one;
	temp = temp - one;
	if (*irnd==0&& temp !=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;
	for (; ;  ) {
		y = z;
		z = y*y;
		/*
		-----------------------------------------------------------------
		        check for underflow here
		-----------------------------------------------------------------
		*/
		a = z*one;
		temp = a + a;
		if (temp==zero||ABS(z)>=y)
			break ;
		i = i+1;
		k = k+k;
	}
	if (*ibeta!=10) {
		*iexp = i+1;
		mx = k+k;
	}
	else {
		/*
		-----------------------------------------------------------------
		     for decimal machines only
		-----------------------------------------------------------------
		*/
		*iexp = 2;
		iz = *ibeta;
		while (k>=iz) {
			iz = *ibeta*iz;
			(*iexp)++ ;
		}
		mx = iz+iz-1;
	}
	for ( ; ;) {
		/*
		-----------------------------------------------------------------
		     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 = a + a;
		if (temp==zero||ABS(y)>=*xmin)
			break;
		k = k+1;
	}
	*minexp = -k;
	/*
	-----------------------------------------------------------------
	     determine *maxexp, *xmax
	-----------------------------------------------------------------
	*/
	if (mx<=k+k-3&&*ibeta!=10) {
		mx = mx+mx;
		(*iexp)++;
	}
	*maxexp = mx+*minexp;
	/*
	-----------------------------------------------------------------
	     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;
	temp = *xmax * one;
	if (temp!=*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;
		}
}