Convert Decimal to Fraction. Version without gotos.

gistfile1.f90

!***********************************************************************************************************************************
!
!                                                         D E C 2 F R A C
!
!  Module:       DEC2FRAC
!
!  Programmer:   David G. Simpson
!                NASA Goddard Space Flight Center
!                Greenbelt, Maryland  20771
!
!  Date:         December 28, 2005
!
!  Language:     Fortran-95
!
!  Version:      1.00b  (December 28, 2005)
!
!  Description:  This program will convert a decimal number to a fraction (i.e. rational number) with increasing
!                degrees of accuracy.  In other words, in computes a series of rational number approximations to
!                any given decimal number.
!
!  Files:        Source files:
!
!                   dec2frac.f95
!
!  Note:         The algorithm is taken from "An Atlas of Functions" by Spanier and Oldham, Springer-Verlag, 1987, pp. 665-667.
!
!***********************************************************************************************************************************
      SUBROUTINE DEC2FRACINT(X,TOL,N1,D1)
      USE ieee_arithmetic
      IMPLICIT NONE
      REAL*8, INTENT(IN) :: X, TOL
      INTEGER, INTENT(OUT) :: N1, D1


      REAL*8 :: NU, R, T, EPS, M
      INTEGER ::  N2,  D2, IDX
      LOGICAL :: SGN

!     Save the sign of X, and make it positive.

      NU = X
      SGN = NU .LT. 0.0D0
      NU = ABS(NU)

!
!     Compute the rational equivalent of X.
!

      D1 = 1
      D2 = 1
      N1 = INT(NU)
      N2 = N1 + 1
      IDX = 0
      EPS = ieee_value(EPS,ieee_positive_inf)
      DO WHILE(EPS .GT. TOL)
          IDX = IDX + 1
          IF (IDX.GT.1) THEN 
              ! recurrence.
              IF (R .LE. 1.0D0) THEN
                  R = 1.0D0/R
              END IF
              N2 = N2 + N1*INT(R)
              D2 = D2 + D1*INT(R)
              N1 = N1 + N2
              D1 = D1 + D2
          END IF

          ! Error ratio of newstep/oldstep
          R = 0.0D0
          IF (NU*D1 .NE. DBLE(N1)) THEN
              R = (N2-NU*D2)/(NU*D1-N1)
          END IF

          ! Swap N1,N2 and D1,D2 if ratio < 1 
          IF (R .LE. 1.0D0) THEN
              T = N2
              N2 = N1
              N1 = T
              T = D2
              D2 = D1
              D1 = T
          END IF

          ! Check realitive error.
          EPS = ABS(1.0D0 - (N1/(NU*D1)))
          IF (EPS .LE. TOL) EXIT

          ! Find first significant digit of error.
          M = 1.0D0
          DO WHILE(M*EPS .LT. 1.0D0) 
            M = 10*M
          END DO
          
          ! round(M*EPS)/M. Check absolute error.
          EPS = (1.0D0/M)*INT(0.5D0+M*EPS)
          IF (EPS .LE. TOL) EXIT
      END DO


      END SUBROUTINE 

      PROGRAM DEC2FRAC
      IMPLICIT NONE

      REAL*8, PARAMETER :: TOL = 1.0D-4
      REAL*8 :: X
      INTEGER :: N1, D1

!-----------------------------------------------------------------------------------------------------------------------------------

!
!     Get user input.
!

      WRITE (UNIT=*, FMT='(A)', ADVANCE='NO') ' Enter decimal number:  '
      READ (UNIT=*, FMT=*) X
      CALL DEC2FRACINT(X,TOL,N1,D1)


      WRITE(*,*),int(sign(dble(n1),X)),d1

      STOP

      END PROGRAM DEC2FRAC