!***********************************************************************************************************************************
!
!                                                         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.
!
!***********************************************************************************************************************************

      PROGRAM DEC2FRAC

      IMPLICIT NONE

      DOUBLE PRECISION, PARAMETER :: TOL = 1.0D-10                                  ! default value of tolerance
      DOUBLE PRECISION :: X, NU, R, T, EPS, M
      INTEGER :: N1, N2, D1, D2
      LOGICAL :: SGN

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

!
!     Get user input.
!

      WRITE (UNIT=*, FMT='(A)', ADVANCE='NO') ' Enter decimal number:  '
      READ (UNIT=*, FMT=*) X

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

      NU = X                                                                        ! make a local copy of X
      SGN = NU .LT. 0.0D0                                                           ! save sign
      NU = ABS(NU)                                                                  ! remove sign from X

!
!     Compute the rational equivalent of X.
!

      D1 = 1
      D2 = 1
      N1 = INT(NU)
      N2 = N1 + 1
      GO TO 300
  100 IF (R .GT. 1.0D0) GO TO 200
      R = 1.0D0/R
  200 N2 = N2 + N1*INT(R)
      D2 = D2 + D1*INT(R)
      N1 = N1 + N2
      D1 = D1 + D2
  300 R = 0.0D0
      IF (NU*D1 .EQ. DBLE(N1)) GO TO 400
      R = (N2-NU*D2)/(NU*D1-N1)
      IF (R .GT. 1.0D0) GO TO 400
      T = N2
      N2 = N1
      N1 = T
      T = D2
      D2 = D1
      D1 = T
  400 IF (SGN) THEN
         WRITE (UNIT=*, FMT='(1X,1H-,I0,1H/,I0)', ADVANCE='NO') N1, D1
      ELSE
         WRITE (UNIT=*, FMT='(1X,I0,1H/,I0)', ADVANCE='NO') N1, D1
      END IF
      WRITE (UNIT=*, FMT='(5X,1H=,ES24.15)', ADVANCE='NO') DBLE(N1)/DBLE(D1)
      EPS = ABS(1.0D0 - (N1/(NU*D1)))
      IF (EPS .LE. TOL) GO TO 600
      M = 1.0D0
  500 M = 10*M
      IF (M*EPS .LT. 1.0D0) GO TO 500
      EPS = (1.0D0/M)*INT(0.5D0+M*EPS)
  600 WRITE (UNIT=*, FMT='(5X,A,ES10.3)') 'error = ', EPS
      IF (EPS .LE. TOL) STOP
      IF (R .NE. 0.0D0) GO TO 100

      STOP

      END PROGRAM DEC2FRAC