home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #18 / NN_1992_18.iso / spool / comp / lang / fortran / 3108 < prev    next >
Encoding:
Internet Message Format  |  1992-08-21  |  19.2 KB
  1. Path: sparky!uunet!dtix!darwin.sura.net!zaphod.mps.ohio-state.edu!uwm.edu!psuvax1!psuvm!hdk
  2. From: HDK@psuvm.psu.edu (H. D. Knoble)
  3. Newsgroups: comp.lang.fortran
  4. Subject: Re: Ignoring Least Significant Bit in a Real Compare
  5. Message-ID: <92234.181300HDK@psuvm.psu.edu>
  6. Date: 21 Aug 92 22:13:00 GMT
  7. References: <1992Aug21.203055.29668@aero.org>
  8.  <1992Aug21.212225.2961@news.eng.convex.com>
  9. Organization: Penn State University
  10. Lines: 440
  11.  
  12.  
  13.                 Floating-Point System Environment Parameters
  14.                        and Tolerant (Fuzzy) Comparisons
  15.                        H. D. Knoble, HDK@PSUVM.PSU.EDU
  16.                   Penn State Center for Academic Computing
  17.  
  18. Fortran algorithms are given that compute tolerant (fuzzy) comparisons of
  19. real numbers and tolerant (fuzzy) Floor/Ceiling Functions.  The driver
  20. FUZZY FORTRAN illustrates use of these tolerant functions; it also shows
  21. to compute 13 floating-point parameters associated with the
  22. floating-point system being utilized by the Fortran host computer.
  23. (Definitions are given in comments of Cody & Waite's Subroutine MACHAR).
  24.  
  25. As a programming practice, one should not compare real numbers for
  26. equality, as quantities computed with different but algebraically
  27. equivalent formulae will hardly ever be computed equal; yet they may be
  28. equal within a floating-point representation (i.e., plus or minus one bit
  29. in the last place of the mantissa).  Tolerant comparisons enable Fortran
  30. algorithms to compare real numbers within this "tolerance" which is
  31. sometimes called the "machine precision" (distance between floating-point
  32. numbers on the unit interval).  Similarly when computing the greatest
  33. integer less than or equal to an argument X, one would not want for
  34. example, 1 for FLOOR(X=1.999999999999999).  Tolerant Floor/Ceiling
  35. Functions resolve this problem in a most sophisticated way (equivalent to
  36. the algorithms used in APL).
  37.  
  38. Algorithms included in FUZZY FORTRAN include:
  39.  
  40. Driver Program    Illustrates use of following subprograms as well as
  41.                   cases where quantities are not equal with respect to
  42.                   .EQ. but are with respect to TEQ; similarly for
  43.                   Floor and TFLOOR. This may be IBM Fortran dependent
  44.                   in that Z-Format (Hexadecimal) is used to display
  45.                   some results.
  46.  
  47. TEQ(X,Y), TNE(X,Y)  Logical Arithmetic Statement Functions to tolerantly
  48.                     compare Double Precision quantities X and Y.
  49.  
  50. TFLOOR(X,CT)  Returns Double Precision whole number that is tolerantly
  51.               the greatest integer algebraically less than X; CT is the
  52.               Comparison Tolerance. This function, small as it is too
  53.               several years to evolve in the APL community.
  54.  
  55. TCEIL(X,CT)   Returns Double Precision whole number that is tolerantly
  56.               the smallest integer algebraically greater than X; CT is
  57.               the Comparison Tolerance.
  58.  
  59. ROUND(X,CT)   Returns Double Precision tolerant proper rounding of X;
  60.               CT is the Comparison Tolerance.
  61.  
  62. MACHAR( parms )   Returns 13 Characteristics of the Machine's
  63.                   floating-point system. Cody and Waite's work.
  64.  
  65. C**********************************************************************
  66. C  ROUTINE:   FUZZY FORTRAN
  67. C  PURPOSE:   Illustrate Hindmarsh's computation of EPS, and APL        t
  68. C             tolerant comparisons, tolerant CEIL/FLOOR, and Tolerant
  69. C             ROUND functions - implemented in Fortran.
  70. C  CALLS:     none
  71. C  PACKAGE:   SAMPLES
  72. C  PLATFORM:  Test system IBM VS Fortran on IBM 370 and up.
  73. C  AUTHOR:    H. D. Knoble <HDK@PSUVM> 09/22/78
  74. C  ORGANIZATION: Penn State University Center for Academic Computing
  75. C  REVISIONS:
  76. C**********************************************************************
  77. C
  78.       DOUBLE PRECISION EPS,EPS3, X,Y,Z, TFLOOR,TCEIL, EPSNEG,XMIN,XMAX
  79.       LOGICAL TEQ,TNE,TGT,TGE,TLT,TLE
  80. C---Following are Fuzzy Comparison (arithmetic statement) Functions.
  81. C
  82.       TEQ(X,Y)=DABS(X-Y).LE.DMAX1(DABS(X),DABS(Y))*EPS3
  83.       TNE(X,Y)=.NOT.TEQ(X,Y)
  84.       TGT(X,Y)=(X-Y).GT.DMAX1(DABS(X),DABS(Y))*EPS3
  85.       TLE(X,Y)=.NOT.TGT(X,Y)
  86.       TLT(X,Y)=TLE(X,Y).AND.TNE(X,Y)
  87.       TGE(X,Y)=TGT(X,Y).OR.TEQ(X,Y)
  88. C
  89. C---Compute EPS for this computer.  EPS is the smallest real number on
  90. C   this architecture such that 1+EPS>1 and 1-EPS<1.
  91.       CALL MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  92.      1             MAXEXP,EPS,EPSNEG,XMIN,XMAX)
  93.       WRITE(6,*) 'Floating-point environment parameters follow...'
  94.       WRITE(6,*) 'IBETA,IT,IRND=',IBETA,IT,IRND
  95.       WRITE(6,*) 'MACHEP,NEGEP,IEXP=',MACHEP,NEGEP,IEXP
  96.       WRITE(6,*) 'MINEXP,MAXEXP=',MINEXP,MAXEXP
  97.       WRITE(6,*) 'EPS,EPSNEG=',EPS,EPSNEG
  98.       WRITE(6,*) 'XMIN,XMAX=',XMIN,XMAX
  99. C---Accept an exact representation, or one bit on either side
  100. C   (heuristic).
  101.       EPS3=3.D0*EPS
  102.       WRITE(6,1) EPS,EPS, EPS3,EPS3
  103. 1     FORMAT(' EPS=',D16.8,2X,Z16, ', EPS3=',D16.8,2X,Z16)
  104.       WRITE(6,*) ' '
  105. C---Illustrate Tolerant Comparisons using EPS3. Any other magnitudes will
  106. C   behave similarly.
  107.       Z=1.D0
  108.       I=49
  109.         X=1.D0/I
  110.         Y=X*I
  111.         WRITE(6,*) 'X=1.D0/',I,', Y=X*',I,', Z=1.D0'
  112.         IF(Y.EQ.Z) WRITE(6,*) 'Fuzzy Comparisons: Y=Z'
  113.         IF(Y.NE.Z) WRITE(6,*) 'Fuzzy Comparisons: Y<>Z'
  114.         IF(TEQ(Y,Z)) THEN
  115.           WRITE(6,*) 'but TEQ(Y,Z) is .TRUE.'
  116.           WRITE(6,*) 'That is, Y is computationally equal to Z.'
  117.         ENDIF
  118.         IF(TNE(Y,Z)) WRITE(6,*) 'and TNE(Y,Z) is .TRUE.'
  119.       WRITE(6,*) ' '
  120. C---Evaluate Fuzzy FLOOR and CEILing Function values using a Comparison
  121. C   Tolerance, CT, of EPS3=3*EPS.
  122.       X=0.1D0
  123.       Y=X*11.D0-X
  124.       YFLOOR=TFLOOR(Y,EPS3)
  125.       YCEIL=TCEIL(Y,EPS3)
  126. 55    Z=1.D0
  127.       WRITE(6,*) 'X=0.1D0, Y=X*11.D0-X, Z=1.D0'
  128.       IF(Y.EQ.Z) WRITE(6,*) 'Fuzzy FLOOR/CEIL: Y=Z'
  129.       IF(Y.NE.Z) WRITE(6,*) 'Fuzzy FLOOR/CEIL: Y<>Z'
  130.       IF(TFLOOR(Y,EPS3).EQ.TCEIL(Y,EPS3).AND.TFLOOR(Y,EPS3).EQ.Z) THEN
  131.         WRITE(6,*) 'but TFLOOR(Y,EPS3)=TCEIL(Y,EPS3)=Z.'
  132.         WRITE(6,*) 'That is, TFLOOR/TCEIL return exact whole numbers.'
  133.       ENDIF
  134.       STOP
  135.       END
  136.       DOUBLE PRECISION FUNCTION TFLOOR(X,CT)
  137. C===========Tolerant FLOOR Function.
  138. C
  139. C    X  -  is given as a double precision argument to be operated on.
  140. C          it is assumed that X is represented with m mantissa bits.
  141. C    CT -  is   given   as   a   Comparison   Tolerance   such   that
  142. C          0.lt.CT.le.3-Sqrt(5)/2. If the relative difference bewteen
  143. C          X and a whole number is  less  than  CT,  then  TFLOOR  is
  144. C          returned   as   this   whole   number.   By  treating  the
  145. C          floating-point numbers as a finite ordered set  note  that
  146. C          the  heuristic  eps=2.**(-(m-1))   and   CT=3*eps   causes
  147. C          arguments  of  TFLOOR/TCEIL to be treated as whole numbers
  148. C          if they are  exactly  whole  numbers  or  are  immediately
  149. C          adjacent to whole number representations.  Since EPS,  the
  150. C          "distance"  between  floating-point  numbers  on  the unit
  151. C          interval, and m, the number of bits in X's mantissa, exist
  152. C          on  every  floating-point   computer,   TFLOOR/TCEIL   are
  153. C          consistently definable on every floating-point computer.
  154. C
  155. C          For more information see the following references:
  156. C    {1} P. E. Hagerty, "More on Fuzzy Floor and Ceiling," APL  QUOTE
  157. C        QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5.
  158. C    {2} L. M. Breed, "Definitions for Fuzzy Floor and Ceiling",  APL
  159. C        QUOTE QUAD 8(3):16-23, March 1978.
  160. C
  161. C   H. D. KNOBLE, Penn State University.
  162. C=====================================================================
  163.       DOUBLE PRECISION X,Q,RMAX,EPS5,CT,FLOOR,DINT
  164. C---------FLOOR(X) is the largest integer algegraically less than
  165. C         or equal to X; that is, the unfuzzy Floor Function.
  166.       DINT(X)=X-DMOD(X,1.D0)
  167.       FLOOR(X)=DINT(X)-DMOD(2.D0+DSIGN(1.D0,X),3.D0)
  168. C---------Hagerty's FL5 Function follows...
  169.       Q=1.D0
  170.       IF(X.LT.0)Q=1.D0-CT
  171.       RMAX=Q/(2.D0-CT)
  172.       EPS5=CT/Q
  173.       TFLOOR=FLOOR(X+DMAX1(CT,DMIN1(RMAX,EPS5*DABS(1.D0+FLOOR(X)))))
  174.       IF(X.LE.0 .OR. (TFLOOR-X).LT.RMAX)RETURN
  175.       TFLOOR=TFLOOR-1.D0
  176.       RETURN
  177.       END
  178.       DOUBLE PRECISION FUNCTION TCEIL(X,CT)
  179. C==========Tolerant Ceiling Function.
  180. C    See TFLOOR.
  181.       DOUBLE PRECISION X,CT,TFLOOR
  182.       TCEIL= -TFLOOR(-X,CT)
  183.       RETURN
  184.       END
  185.       DOUBLE PRECISION FUNCTION ROUND(X,CT)
  186. C=========Tolerant Round Function
  187. C  See Knuth, Art of Computer Programming, Vol. 1, Problem 1.2.4-5.
  188.       DOUBLE PRECISION TFLOOR,X,CT
  189.       ROUND=TFLOOR(X+0.5D0,CT)
  190.       RETURN
  191.       END
  192.        SUBROUTINE MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  193.      1                    MAXEXP,EPS,EPSNEG,XMIN,XMAX)
  194.        IMPLICIT REAL*8 (A-H,O-Z)
  195. C===See "SOFTWARE MANUAL FOR THE ELEMENTARY FUNCTIONS" by William J.
  196. C   Cody, Jr. and William Waite, Prentice Hall, ISBN 0-13-822064,
  197. C   1980, pages 258-264.
  198. C
  199. C   The subroutine MACHAR is intended to dynamically determine
  200. C   thirteen parameters relating to the floating-point
  201. C   arithmetic system.  The first three--the radix B, the number
  202. C   of base-B digits in the significand, and the question of
  203. C   whether floating-point addition rounds or truncates--are
  204. C   determined with an algorithm developed by M.  Malcolm
  205. C   [1972].  Several of the remaining parameters are determined
  206. C   with "folk algorithms," but most of them are determined
  207. C   using algorithms we believe to be new.  The methods are
  208. C   sometimes ad hoc, and the program is not guaranteed to work
  209. C   on all machines.  Not all of the conditions under which
  210. C   MACHAR will fail are known, but it is known to fail on
  211. C   Honeywell machines where the active arithmetic registers are
  212. C   wider than storage registers.  This particular failure can
  213. C   be corrected by forcing the storage of intermediate results
  214. C   at critical points (see Gentleman and Marovich [1974]).  We
  215. C   have deliberately not followed that procedure in MACHAR
  216. C   because the program is already very complicated and
  217. C   difficult to understand.
  218. C
  219. C   Aside from the Honeywell machines the program has run
  220. C   properly on all other machines we have tried.  These
  221. C   include, in alphabetical order, BESM-6 (Russian), Burroughs
  222. C   6700, CDC 6000/7000 series, CRAY-1, IBM 360/370 series, IBM
  223. C   3033, PDP-10, PDP-11, Univac 1105, Varian V76, and Z-80
  224. C   (with Microsoft fortran).  MACHAR has been translated into
  225. C   Basic and run on the Datapoint 2200, IBM 5100 and PDP-10.
  226. C   It has even been translated and run on the HP-67 and TI-59
  227. C   programmable calculators.  Despite these successes we
  228. C   strongly urge that parameter values returned by this program
  229. C   be examined critically the first time it is run on any
  230. C   machine.
  231. C
  232.        INTEGER I,IBETA,IEXP,IRND,IT,IZ,J,K,MACHEP,MAXEXP,MINEXP,
  233.      1         MX,NEGEP,NGRD
  234. CE     REAL A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN,Y,Z,ZERO
  235.        DOUBLE PRECISION A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,
  236.      1                  XMIN,Y,Z,ZERO
  237. C
  238. C      This subroutine is intended to determine the characteristics
  239. C      of the floating-point arithmetic system that are specified
  240. C      below.  The first three are determined according to an
  241. C      algorithm due to M. Malcolm, CACM 15 (1972), pp. 949-951,
  242. C      incorporating some, but not all, of the improvements
  243. C      suggested by M. Gentleman and S. Marovich, CACM 17 (1974),
  244. C      pp. 276-277.  The version given here is for double precision.
  245. C      Lines containing  CE  in columns 1 and 2 can be used to
  246. C      convert the subroutine to single precision by replacing
  247. C      existing lines in the obvious manner.
  248. C
  249. C
  250. C      IBETA   - the radix of the floating-point representation
  251. C      IT      - the number of base ibeta digits in the floating-point
  252. C                significand
  253. C      IRND    - 0 if floating-point addition chops,
  254. C                1  if floating-point addition rounds
  255. C      NGRD    - the number of guard digits for multiplication.  It is
  256. C                0 if  IRND=1, or if  IRND=0  and only  it  base IBETA
  257. C                  digits participate in the post normalization shift
  258. C                  of the floating-point significand in multiplication
  259. C                1 if  IRND=0  and more than  it  base  IBETA  digits
  260. C                  participate in the post normalization shift of the
  261. C                  floating-point significand in multiplication.
  262. C      MACHEP  - the largest negative integer such that
  263. C                1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
  264. C                MACHEP is bounded below by  -(IT+3)
  265. C      NEGEPS  - the largest negative integer such that
  266. C                1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
  267. C                NEGEPS is bounded below by  -(IT+3)
  268. C      IEXP    - the number of bits (decimal places if IBETA = 10)
  269. C                reserved for the representation of the exponent
  270. C                (including the bias or sign) of a floating-point
  271. C                number
  272. C      MINEXP  - the largest in magnitude negative integer such that
  273. C                FLOAT(IBETA)**MINEXP is a positive floating-point
  274. C                number
  275. C      MAXEXP  - the largest positive integer exponent for a finite
  276. C                floating-point number
  277. C      EPS     - the smallest positive floating-point number such
  278. C                that  1.0+EPS .ne.  1.0. in particular, if either
  279. C                IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
  280. C                otherwise,  EPS =  (FLOAT(IBETA)**MACHEP)/2
  281. C      EPSNEG  - a small positive floating-point number such that
  282. C                1.0-EPSNEG .ne. 1.0. in particular, if IBETA = 2
  283. C                or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
  284. C                otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
  285. C                NEGEPS is bounded below by -(IT+3), EPSNEG may not
  286. C                be the smallest number which can alter 1.0 by
  287. C                subtraction.
  288. C      XMIN    - the smallest non-vanishing floating-point power of
  289. C                the radix. in particular, XMIN = FLOAT(IBETA)**MINEXP
  290. C      XMAX    - the largest finite floating-point number.  In
  291. C                particular   XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
  292. C                Note - on some machines  XMAX  will be only the
  293. C                second, or perhaps third, largest number, being
  294. C                too small by 1 or 2 units in the last digit of
  295. C                the significand.
  296. C
  297. C      Latest revision - October 22, 1979
  298. C
  299. C      Author - W. J. Cody
  300. C               Argonne National Laboratory
  301. C
  302. C-----------------------------------------------------------------
  303. CE     ONE = FLOAT(1)
  304.        ONE = DBLE(FLOAT(1))
  305. CE     ZERO = 0.0E0
  306.        ZERO = 0.0D0
  307. C-----------------------------------------------------------------
  308. C      Determine IBETA,BETA a la Malcolm
  309. C-----------------------------------------------------------------
  310.        A = ONE
  311.    10  A = A + A
  312.          IF (((A+ONE)-A)-ONE .EQ. ZERO) GO TO 10
  313.        B = ONE
  314.    20  B = B + B
  315.          IF ((A+B)-A .EQ. ZERO) GO TO 20
  316. CE    IBETA = INT((A+B)-A)
  317.       IBETA = INT(SNGL((A + B) - A))
  318. CE    BETA = FLOAT(IBETA)
  319.       BETA = DBLE(FLOAT(IBETA))
  320.  
  321. C     Determine IT, IRND
  322. C-----------------------------------------------------------------
  323.       IT = 0
  324.       B  = ONE
  325.   100 IT = IT + 1
  326.          B = B * BETA
  327.          IF (((B+ONE)-B)-ONE .EQ. ZERO) GO TO 100
  328.       IRND = 0
  329.       BETAM1 = BETA - ONE
  330.       IF ((A+BETAM1)-A .NE. ZERO) IRND = 1
  331. C-----------------------------------------------------------------
  332. C     Determine NEGEP, EPSNEG
  333. C-----------------------------------------------------------------
  334.       NEGEP = IT + 3
  335.       BETAIN = ONE / BETA
  336.       A  = ONE
  337. C
  338.       DO 200 I = 1, NEGEP
  339.          A = A * BETAIN
  340.   200 CONTINUE
  341. C
  342.       B = A
  343.   210 IF ((ONE-A)-ONE .NE. ZERO) GO TO 220
  344.          A = A * BETA
  345.          NEGEP = NEGEP - 1
  346.       GO TO 210
  347.   220 NEGEP = -NEGEP
  348.       EPSNEG = A
  349.       IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 300
  350.       A = (A*(ONE+A)) / (ONE+ONE)
  351.       IF ((ONE-A)-ONE .NE. ZERO) EPSNEG = A
  352.  
  353. C     Determine MACHEP, EPS
  354. C-----------------------------------------------------------------
  355.   300 MACHEP = -IT - 3
  356.       A = B
  357.   310 IF((ONE+A)-ONE .NE. ZERO) GO TO 320
  358.           A = A * BETA
  359.           MACHEP = MACHEP + 1
  360.        GO TO 310
  361.   320  EPS = A
  362.        IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 350
  363.        A = (A*(ONE+A)) / (ONE+ONE)
  364.        IF ((ONE+A)-ONE .NE. ZERO) EPS = A
  365. C-----------------------------------------------------------------
  366. C      Determine NGRD
  367. C-----------------------------------------------------------------
  368.   350  NGRD = 0
  369.        IF ((IRND .EQ. 0) .AND. ((ONE+EPS)*ONE-ONE) .NE. ZERO) NGRD = 1
  370. C-----------------------------------------------------------------
  371. C      Determine IEXP, MINEXP, XMIN
  372. C
  373. C      Loop to determine largest I and K = 2**I such that
  374. C         (1/BETA) ** (2**(1))
  375. C      does not underflow.
  376. C      Exit from loop is signaled by an underflow.
  377. C-----------------------------------------------------------------
  378. C---VS Fortran: allow two underflows with no messages.
  379. C      CALL ERRSET(208,3,-1,+1)
  380.        I = 0
  381.        K = 1
  382.        Z = BETAIN
  383.   400  Y = Z
  384.           Z = Y * Y
  385. C-----------------------------------------------------------------
  386. C         Check for Underflow here
  387. C-----------------------------------------------------------------
  388.           A = Z * ONE
  389. CE        IF ((A+A .EQ. ZERO) .OR.  (ABS(Z) .GE. Y)) GO TO 410
  390.           IF ((A+A .EQ. ZERO) .OR.  (DABS(Z) .GE. Y)) GO TO 410
  391.           I = I + 1
  392.           K = K + K
  393.        GO TO 400
  394.   410  IF (IBETA .EQ.  10) GO TO 420
  395.        IEXP = I + 1
  396.        MX = K + K
  397.        GO TO 450
  398. C-----------------------------------------------------------------
  399. C      For decimal machines only
  400. C-----------------------------------------------------------------
  401.   420  IEXP = 2
  402.        IZ = IBETA
  403.   430  IF (K .LT. IZ) GO TO 440
  404.           IZ = IZ * IBETA
  405.           IEXP = IEXP + 1
  406.        GO TO 430
  407.   440  MX = IZ + IZ - 1
  408. C-----------------------------------------------------------------
  409. C      Loop to determine MINEXP, XMIN
  410. C      Exit from loop is signaled by an underflow.
  411. C-----------------------------------------------------------------
  412.   450  XMIN = Y
  413.           Y = Y * BETAIN
  414. C-----------------------------------------------------------------
  415. C         Check for underflow here
  416. C-----------------------------------------------------------------
  417.           A = Y * ONE
  418. CE        IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460
  419.           IF (((A+A) .EQ. ZERO) .OR. (DABS(Y)  .GE. XMIN)) GO TO 460
  420.           K = K + 1
  421.        GO TO 450
  422.   460 MINEXP = -K
  423. C-----------------------------------------------------------------
  424. C      Determine MAXEXP, XMAX
  425. C-----------------------------------------------------------------
  426.        IF ((MX .GT. K+K-3) .OR.  (IBETA .EQ. 10)) GO TO 500
  427.        MX = MX + MX
  428.        IEXP = IEXP + 1
  429.   500  MAXEXP = MX + MINEXP
  430. C-----------------------------------------------------------------
  431. C      Adjust for machines with implicit leading
  432. C      bit in binary significant and machines with
  433. C      radix point at extreme right of significand
  434. C-----------------------------------------------------------------
  435.        I = MAXEXP + MINEXP
  436.        IF ((IBETA .EQ. 2) .AND.  (I .EQ. 0)) MAXEXP = MAXEXP - 1
  437.        IF (I .GT. 20) MAXEXP = MAXEXP - 1
  438.        IF (A .NE. Y) MAXEXP = MAXEXP - 2
  439.        XMAX = ONE - EPSNEG
  440.        IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
  441.        XMAX = XMAX / (BETA * BETA * BETA * XMIN)
  442.        I = MAXEXP + MINEXP + 3
  443.        IF (I .LE. 0) GO TO 520
  444. C
  445.        DO 510 J = 1, I
  446.            IF (IBETA .NE. 2) XMAX = XMAX * BETA
  447.   510  CONTINUE
  448. C
  449.   520  RETURN
  450. C     ---------- Last line of MACHAR ----------
  451.       END
  452.