/*

This file combines the single and double precision versions of machar,
selected by cc -DSP or cc -DDP.  This feature provided by D. G. Hough,
August 3, 1988.

*/

#ifdef SP
#define REAL float
#define ZERO 0.0
#define ONE 1.0
#define PREC "Single "
#define REALSIZE 1
#endif
 
#ifdef DP
#define REAL double
#define ZERO 0.0e0
#define ONE 1.0e0
#define PREC "Double "
#define REALSIZE 2
#endif
 
#include <math.h>
#include <stdio.h>

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

void
rmachar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,
        maxexp,eps,epsneg,xmin,xmax)

      int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd;
      REAL *eps,*epsneg,*xmax,*xmin;

/*

   This 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 present program is a
    translation of the Fortran 77 program in W. J. Cody, "MACHAR:
    A subroutine to dynamically determine machine parameters".
    TOMS (14), 1988.
 
   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 - August 4, 1988
 
      Author - W. J. Cody
               Argonne National Laboratory
 
*/

{
      int i,iz,j,k;
      int mx,itmp,nxres;
      REAL a,b,beta,betain,one,y,z,zero;
      REAL betah,t,tmp,tmpa,tmp1,two;

      (*irnd) = 1;
      one = (REAL)(*irnd);
      two = one + one;
      a = two;
      b = a;
      zero = 0.0e0;

/*
  determine ibeta,beta ala malcolm
*/

      tmp = ((a+one)-a)-one;

      while (tmp == zero) {
         a = a+a;
         tmp = a+one;
         tmp1 = tmp-a;
         tmp = tmp1-one;
      }

      tmp = a+b;
      itmp = (int)(tmp-a);
      while (itmp == 0) {
         b = b+b;
         tmp = a+b;
         itmp = (int)(tmp-a);
      }

      *ibeta = itmp;
      beta = (REAL)(*ibeta);

/*
  determine irnd, it
*/

      (*it) = 0;
      b = one;
      tmp = ((b+one)-b)-one;

      while (tmp == zero) {
         *it = *it+1;
         b = b*beta;
         tmp = b+one;
         tmp1 = tmp-b;
         tmp = tmp1-one;
      }

      *irnd = 0;
      betah = beta/two;
      tmp = a+betah;
      tmp1 = tmp-a;
      if (tmp1 != zero) *irnd = 1;
      tmpa = a+beta;
      tmp = tmpa+betah;
      if ((*irnd == 0) && (tmp-tmpa != 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;
      tmp = (one-a);
      tmp = tmp-one;

      while (tmp == zero) {
         a = a*beta;
         *negep = *negep-1;
         tmp1 = one-a;
         tmp = tmp1-one;
      }

      (*negep) = -(*negep);
      (*epsneg) = a;

/*
  determine machep, eps
*/

      (*machep) = -(*it) - 3;
      a = b;
      tmp = one+a;

      while (tmp-one == zero) {
         a = a*beta;
         *machep = *machep+1;
         tmp = one+a;
      }

      *eps = a;
      
/*
  determine ngrd
*/

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

/*
  determine iexp, minexp, xmin

  loop to determine largest 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
*/

         a = z * one;
         tmp = z*t;
         if ((a+a == zero) || (ABS(z) > y)) break;
         tmp1 = tmp*betain;
         if (tmp1*beta == z) break;
         i = i + 1;
         k = k+k;
      }

/*
  determine k such that (1/beta)**k does not underflow
    first set  k = 2 ** i
*/

      (*iexp) = i + 1;
      mx = k + k;
      if (*ibeta == 10) {

/*
  for decimal machines only
*/

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

      for (;;) {
         (*xmin) = y;
         y = y * betain;
         a = y * one;
         tmp = y*t;
         tmp1 = a+a;
         if ((tmp1 == zero) || (ABS(y) >= (*xmin))) break;
         k = k + 1;
         tmp1 = tmp*betain;
         tmp1 = tmp1*beta;

         if ((tmp1 == y) && (tmp != y)) {
            nxres = 3;
            *xmin = y;
            break;
         }

      }

      (*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) (*maxexp) = (*maxexp)-2;

/*
  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);
      tmp = (*xmax)*one;
      if (tmp != (*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;
         }

      }
 
    return;

}

typedef union
{
  double d;
  int i[2];
} equiv;

int
main (void)
{
  /* Works for 32 bit machines with 32 bit ints and 64 bit doubles */

  int ibeta, iexp, irnd, it, machep, maxexp, minexp, negep, ngrd;
  REAL eps, epsneg, xmax, xmin;
  int i;
  equiv flt_params[6];

  rmachar (&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
	   &maxexp, &eps, &epsneg, &xmin, &xmax);

  flt_params[1].d = xmin;
  flt_params[2].d = xmax;
  flt_params[3].d = epsneg;
  flt_params[4].d = eps;
  flt_params[5].d = log10 ((double) ibeta);

  printf ("* d1mach.f  Do not edit.  Generated automatically by gen-d1mach.c\n\
      double precision function d1mach(i)\n\
      integer i\n\
      integer i1var (2)\n\
      integer i2var (2)\n\
      integer i3var (2)\n\
      integer i4var (2)\n\
      integer i5var (2)\n\
      double precision dmach(5)\n\
      equivalence (dmach(1), i1var(1))\n\
      equivalence (dmach(2), i2var(1))\n\
      equivalence (dmach(3), i3var(1))\n\
      equivalence (dmach(4), i4var(1))\n\
      equivalence (dmach(5), i5var(1))\n");

  for (i = 1; i < 6; i++)
    printf ("      data i%dvar(1), i%dvar(2) / %ld , %ld /\n",
	    i, i, flt_params[i].i[0], flt_params[i].i[1]);

  printf ("      if (i .lt. 1  .or.  i .gt. 5) goto 999\n\
      d1mach = dmach(i)\n\
      return\n\
  999 write(*,1999) i\n\
 1999 format(' d1mach - i out of bounds', i10)\n\
      call xstopx (' ')\n\
      end\n");

  return 0;
}