*comment
	declarations of various floating point type letters
	with corresponding byte sizes for real type declarations
	and format declarations
	all as lists
*endcomment
*mlist(sizes,size_real_list,size_format_list)
'F'	'*4'	'e20.10'
'D'	'*8'	'e20.10'
'G'	'*8'	'e20.10'
'H'	'*16'	'e50.35e6'
*endmlist
*set(list_size=length(sizes))
*comment
	declaration section done-start main program
*endcomment
	PROGRAM TEST_CHAR
	IMPLICIT NONE
*do(index_counter=1,list_size)
	call do_$sizes
*enddo
	stop
	end
*page
*comment
	define macro to call check routine and output results
*endcomment
*macro(execute_size_check)
	*if(size_letter.eq.'G')
	 options/g_floating
	*endif	
	subroutine do_$size_letter
	integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
	real$size_real eps,epsneg,xmin,xmax
      	call size_letter$_MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1                   MAXEXP,EPS,EPSNEG,XMIN,XMAX)
	write(6,100)IBETA,IT,iexp,ngrd
100	format(	' size_letter$_Floating point',/,
	1	' Radix ',t30,I10,/,
	1	' Digits in base in mantisa',t30,I10,/,
	1	' Digits in base in exponent',t30,i10,/,
	1	' Guard digits in base in mantisa',t30,I10)
	if(irnd.eq.0)then
		write(6,101)
101		format(' F.P. Add chops')
	elseif(irnd.eq.1)then
		write(6,102)
102		format(' F.P. Add rounds, not in the IEEE style')
	elseif(irnd.eq.2)then
		write(6,103)
103		format(' F.P. Add rounds in the IEEE style')
	elseif(irnd.eq.3)then
		write(6,104)
104		format(' F.P. Add chops, partial underflow')
	elseif(irnd.eq.4)then
		write(6,105)
105	format(' F.P. Add rounds,not in the IEEE style,partial underflow')
	elseif(irnd.eq.5)then
		write(6,106)
106	format(' F.P. Add rounds in the IEEE style,partial underflow')
	else
		write(6,107)
107		format(' Unknown rounding and underflow')
	endif
	write(6,108)ibeta,MACHEP,eps,
	1	ibeta,NEGEP,epsneg,
	1	ibeta,MINEXP,XMIN,
     1                   ibeta,MAXEXP,XMAX
108	format(
	1 ' Smallest no. such that 1.+N.ne.1.0 =(',I3,'**',I10,')',
	1 size_format,/,
	1 ' Smallest no. such that 1.-N.ne.1.0 =(',I3,'**',I10,')',
	1 size_format,/,
	1 ' min non zero number is ',I3,'**',I10,' = ',size_format,/,
	1 ' Min number that overflows is ',I3,'**',I10,/,
	1 ' Largest rep number ',size_format)
	return
	end
*endmacro
*do(index_counter=1,list_size)
 *execute_size_check(size_letter='$sizes$')(size_real='$size_real_list$')(size_format='$size_format_list$')
*enddo
*page
*comment
	define macro to generate routine to do actuall check
*endcomment
*macro(get_machine_char)
	*if(size_letter.eq.'G')
	 options/g_floating
	*endif	
      SUBROUTINE size_letter$_MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1                   MAXEXP,EPS,EPSNEG,XMIN,XMAX)
C-----------------------------------------------------------------------
C  THIS FORTRAN 77 SUBROUTINE IS INTENDED TO DETERMINE THE PARAMETERS
C   OF THE FLOATING-POINT ARITHMETIC SYSTEM SPECIFIED BELOW.  THE
C   DETERMINATION OF THE FIRST THREE USES AN EXTENSION OF AN ALGORITHM
C   DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, INCORPORATING SOME,
C   BUT NOT ALL, OF THE IMPROVEMENTS SUGGESTED BY M. GENTLEMAN AND S.
C   MAROVICH, CACM 17 (1974), PP. 276-277.  AN EARLIER VERSION OF THIS
C   PROGRAM WAS PUBLISHED IN THE BOOK SOFTWARE MANUAL FOR THE
C   ELEMENTARY FUNCTIONS BY W. J. CODY AND W. WAITE, PRENTICE-HALL,
C   ENGLEWOOD CLIFFS, NJ, 1980.  THE PRESENT VERSION IS DOCUMENTED IN 
C   W. J. CODY, "MACHAR: A SUBROUTINE TO DYNAMICALLY DETERMINE MACHINE
C   PARAMETERS," TOMS 14, DECEMBER, 1988.
C
C  THE PROGRAM AS GIVEN HERE MUST BE MODIFIED BEFORE COMPILING.  IF
C   A SINGLE (DOUBLE) PRECISION VERSION IS DESIRED, CHANGE ALL
C   OCCURRENCES OF CS (CD) IN COLUMNS 1 AND 2 TO BLANKS.
C
C  PARAMETER VALUES REPORTED ARE AS FOLLOWS:
C
C       IBETA   - THE RADIX FOR THE FLOATING-POINT REPRESENTATION
C       IT      - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING-POINT
C                 SIGNIFICAND
C       IRND    - 0 IF FLOATING-POINT ADDITION CHOPS
C                 1 IF FLOATING-POINT ADDITION ROUNDS, BUT NOT IN THE
C                   IEEE STYLE
C                 2 IF FLOATING-POINT ADDITION ROUNDS IN THE IEEE STYLE
C                 3 IF FLOATING-POINT ADDITION CHOPS, AND THERE IS
C                   PARTIAL UNDERFLOW
C                 4 IF FLOATING-POINT ADDITION ROUNDS, BUT NOT IN THE
C                   IEEE STYLE, AND THERE IS PARTIAL UNDERFLOW
C                 5 IF FLOATING-POINT ADDITION ROUNDS IN THE IEEE STYLE,
C                   AND THERE IS PARTIAL UNDERFLOW
C       NGRD    - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION WITH
C                 TRUNCATING ARITHMETIC.  IT IS
C                 0 IF FLOATING-POINT ARITHMETIC ROUNDS, OR IF IT
C                   TRUNCATES AND ONLY  IT  BASE  IBETA DIGITS
C                   PARTICIPATE IN THE POST-NORMALIZATION SHIFT OF THE
C                   FLOATING-POINT SIGNIFICAND IN MULTIPLICATION;
C                 1 IF FLOATING-POINT ARITHMETIC TRUNCATES AND MORE
C                   THAN  IT  BASE  IBETA  DIGITS PARTICIPATE IN THE
C                   POST-NORMALIZATION SHIFT OF THE FLOATING-POINT
C                   SIGNIFICAND IN MULTIPLICATION.
C       MACHEP  - THE LARGEST NEGATIVE INTEGER SUCH THAT
C                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, EXCEPT THAT
C                 MACHEP IS BOUNDED BELOW BY  -(IT+3)
C       NEGEPS  - THE LARGEST NEGATIVE INTEGER SUCH THAT
C                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, EXCEPT THAT
C                 NEGEPS IS BOUNDED BELOW BY  -(IT+3)
C       IEXP    - THE NUMBER OF BITS (DECIMAL PLACES IF IBETA = 10)
C                 RESERVED FOR THE REPRESENTATION OF THE EXPONENT
C                 (INCLUDING THE BIAS OR SIGN) OF A FLOATING-POINT
C                 NUMBER
C       MINEXP  - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT
C                 FLOAT(IBETA)**MINEXP IS POSITIVE AND NORMALIZED
C       MAXEXP  - THE SMALLEST POSITIVE POWER OF  BETA  THAT OVERFLOWS
C       EPS     - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH
C                 THAT  1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER
C                 IBETA = 2  OR  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
C                 OTHERWISE,  EPS = (FLOAT(IBETA)**MACHEP)/2
C       EPSNEG  - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT
C                 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2
C                 OR  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
C                 OTHERWISE,  EPSNEG = (IBETA**NEGEPS)/2.  BECAUSE
C                 NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT
C                 BE THE SMALLEST NUMBER THAT CAN ALTER 1.0 BY
C                 SUBTRACTION.
C       XMIN    - THE SMALLEST NON-VANISHING NORMALIZED FLOATING-POINT
C                 POWER OF THE RADIX, I.E.,  XMIN = FLOAT(IBETA)**MINEXP
C       XMAX    - THE LARGEST FINITE FLOATING-POINT NUMBER.  IN
C                 PARTICULAR  XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
C                 NOTE - ON SOME MACHINES  XMAX  WILL BE ONLY THE
C                 SECOND, OR PERHAPS THIRD, LARGEST NUMBER, BEING
C                 TOO SMALL BY 1 OR 2 UNITS IN THE LAST DIGIT OF
C                 THE SIGNIFICAND.
C
C     LATEST REVISION - DECEMBER 4, 1987
C
C     AUTHOR - W. J. CODY
C              ARGONNE NATIONAL LABORATORY
C
C-----------------------------------------------------------------------
      INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP,
     1        MINEXP,MX,NEGEP,NGRD,NXRES
CD    DOUBLE PRECISION
      REAL$size_real
     1   A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA,
     2   TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
C-----------------------------------------------------------------------
*if(size_letter.eq.'F')
      CONV(I) = REAL(I)
*elseif(size_letter.eq.'D')
      CONV(I) = DBLE(I)
*elseif(size_letter.eq.'G')
      CONV(I) = DBLE(I)
*elseif(size_letter.eq.'H')
	conv(i)=qext(i)
*endif
      ONE = CONV(1)
      TWO = ONE + ONE
      ZERO = ONE - ONE
C-----------------------------------------------------------------------
C  DETERMINE IBETA, BETA ALA MALCOLM.
C-----------------------------------------------------------------------
      A = ONE
   10 A = A + A
         TEMP = A+ONE
         TEMP1 = TEMP-A
         IF (TEMP1-ONE .EQ. ZERO) GO TO 10
      B = ONE
   20 B = B + B
         TEMP = A+B
         ITEMP = INT(TEMP-A)
         IF (ITEMP .EQ. 0) GO TO 20
      IBETA = ITEMP
      BETA = CONV(IBETA)
C-----------------------------------------------------------------------
C  DETERMINE IT, IRND.
C-----------------------------------------------------------------------
      IT = 0
      B = ONE
  100 IT = IT + 1
         B = B * BETA
         TEMP = B+ONE
         TEMP1 = TEMP-B
         IF (TEMP1-ONE .EQ. ZERO) GO TO 100
      IRND = 0
      BETAH = BETA / TWO
      TEMP = A+BETAH
      IF (TEMP-A .NE. ZERO) IRND = 1
      TEMPA = A + BETA
      TEMP = TEMPA+BETAH
      IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2
C-----------------------------------------------------------------------
C  DETERMINE NEGEP, EPSNEG.
C-----------------------------------------------------------------------
      NEGEP = IT + 3
      BETAIN = ONE / BETA
      A = ONE
      DO 200 I = 1, NEGEP
         A = A * BETAIN
  200 CONTINUE
      B = A
  210 TEMP = ONE-A
         IF (TEMP-ONE .NE. ZERO) GO TO 220
         A = A * BETA
         NEGEP = NEGEP - 1
      GO TO 210
  220 NEGEP = -NEGEP
      EPSNEG = A
C-----------------------------------------------------------------------
C  DETERMINE MACHEP, EPS.
C-----------------------------------------------------------------------
      MACHEP = -IT - 3
      A = B
  300 TEMP = ONE+A
         IF (TEMP-ONE .NE. ZERO) GO TO 320
         A = A * BETA
         MACHEP = MACHEP + 1
      GO TO 300
  320 EPS = A
C-----------------------------------------------------------------------
C  DETERMINE NGRD.
C-----------------------------------------------------------------------
      NGRD = 0
      TEMP = ONE+EPS
      IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1
C-----------------------------------------------------------------------
C  DETERMINE IEXP, MINEXP, XMIN.
C
C  LOOP TO DETERMINE LARGEST I AND K = 2**I SUCH THAT
C         (1/BETA) ** (2**(I))
C  DOES NOT UNDERFLOW.
C  EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW.
C-----------------------------------------------------------------------
      I = 0
      K = 1
      Z = BETAIN
      T = ONE + EPS
      NXRES = 0
  400 Y = Z
         Z = Y * Y
C-----------------------------------------------------------------------
C  CHECK FOR UNDERFLOW HERE.
C-----------------------------------------------------------------------
         A = Z * ONE
         TEMP = Z * T
         IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410
         TEMP1 = TEMP * BETAIN
         IF (TEMP1*BETA .EQ. Z) GO TO 410
         I = I + 1
         K = K + K
      GO TO 400
  410 IF (IBETA .EQ. 10) GO TO 420
      IEXP = I + 1
      MX = K + K
      GO TO 450
C-----------------------------------------------------------------------
C  THIS SEGMENT IS FOR DECIMAL MACHINES ONLY.
C-----------------------------------------------------------------------
  420 IEXP = 2
      IZ = IBETA
  430 IF (K .LT. IZ) GO TO 440
         IZ = IZ * IBETA
         IEXP = IEXP + 1
      GO TO 430
  440 MX = IZ + IZ - 1
C-----------------------------------------------------------------------
C  LOOP TO DETERMINE MINEXP, XMIN.
C  EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW.
C-----------------------------------------------------------------------
  450 XMIN = Y
         Y = Y * BETAIN
C-----------------------------------------------------------------------
C  CHECK FOR UNDERFLOW HERE.
C-----------------------------------------------------------------------
         A = Y * ONE
         TEMP = Y * T
         IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460
         K = K + 1
         TEMP1 = TEMP * BETAIN
         IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN
               GO TO 450
            ELSE
               NXRES = 3
               XMIN = Y
         END IF
  460 MINEXP = -K
C-----------------------------------------------------------------------
C  DETERMINE MAXEXP, XMAX.
C-----------------------------------------------------------------------
      IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500
      MX = MX + MX
      IEXP = IEXP + 1
  500 MAXEXP = MX + MINEXP
C-----------------------------------------------------------------
C  ADJUST IRND TO REFLECT PARTIAL UNDERFLOW.
C-----------------------------------------------------------------
      IRND = IRND + NXRES
C-----------------------------------------------------------------
C  ADJUST FOR IEEE-STYLE MACHINES.
C-----------------------------------------------------------------
      IF (IRND .GE. 2) MAXEXP = MAXEXP - 2
C-----------------------------------------------------------------
C  ADJUST FOR MACHINES WITH IMPLICIT LEADING BIT IN BINARY
C  SIGNIFICAND, AND MACHINES WITH RADIX POINT AT EXTREME
C  RIGHT OF SIGNIFICAND.
C-----------------------------------------------------------------
      I = MAXEXP + MINEXP
      IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1
      IF (I .GT. 20) MAXEXP = MAXEXP - 1
      IF (A .NE. Y) MAXEXP = MAXEXP - 2
      XMAX = ONE - EPSNEG
      IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
      XMAX = XMAX / (BETA * BETA * BETA * XMIN)
      I = MAXEXP + MINEXP + 3
      IF (I .LE. 0) GO TO 520
      DO 510 J = 1, I
          IF (IBETA .EQ. 2) XMAX = XMAX + XMAX
          IF (IBETA .NE. 2) XMAX = XMAX * BETA
  510 CONTINUE
  520 RETURN
C---------- LAST CARD OF MACHAR ----------
      END
*endmacro
*do(index_counter=1,list_size)
 *get_machine_char(size_letter='$SIZEs$')(size_real='$size_real_list$')
*enddo
*end