*=======================================================================
*   Euclid's Algorithm for finding the Greatest Common Divisor 
*
*   refer to the Algol procedure by Robert Claussen, 1957
*   General Electric Co., Cincinnati, Ohio
*   published as Algorithm #7 of the ACM
*   This implementation by Andy Allinger, 2011, public domain
*
* 	Each pair of integers must have a greatest common divisor (GCD)
*  Input : R, S                 integers 
*  Output: R                    junk
*          S                    greatest common divisor
*=======================================================================
       SUBROUTINE EUCLID (R, S)
 
        IMPLICIT NONE
        INTEGER Q, R, S

*        Require positive input
       IF ((R < 0) .OR. (S < 0)) 
     &            STOP 'euclid: input must not be negative'
       IF ((R .EQ. 0) .OR. (S .EQ. 0)) 
     &            STOP 'euclid: input must not be zero'

*       Ensure that R > S
       IF (R < S ) THEN     ! swap
         Q = R
         R = S 
         S = Q
       END IF

*        Take the remainder Q = R % S
  50   Q = MOD(R, S)
       IF (Q .EQ. 0) RETURN
       R = S 
       S = Q 
       GO TO 50

      END ! of euclid

      
*=======================================================================
* GCD - the greatest common divisor
*      A            array of integers
*      N            size of the array
*=======================================================================
      FUNCTION GCD (A, N)
      
       IMPLICIT NONE
       INTEGER A, N, GCD 
       DIMENSION A(N)
       
       INTEGER J, K
      
       GCD = A(1)
       DO J = 2, N
         K = A(J)
         CALL EUCLID (K, GCD)
       END DO   
      
      END ! of GCD
      
      
*=======================================================================
* LCD - the lowest common denominator 
*        A           array of integers
*        N           size of the array
*        fault       error code, 1 = overflow
*=======================================================================
      FUNCTION LCD (A, N, fault)
      
       IMPLICIT NONE
       INTEGER A, N, fault, LCD 
       DIMENSION A(N)
       
       INTEGER I, J, K, L, BITSZ
       REAL lmax
       LOGICAL first
       SAVE lmax, first
       DATA first / .TRUE. / 
        
*            comparison value to test for overflow
       IF (first) THEN
         lmax = LOG( (2.0 **(BITSZ(I)-1) -1.) )
         first = .FALSE.
       END IF
       
       fault = 0
       LCD = A(1)
       DO I = 2, N
         J = A(I)
         K = J
         L = LCD
         CALL EUCLID (L, K)
         
*             ask if LCD @ (J / K) may be computed
         IF (LOG( FLOAT(LCD) * FLOAT(J) / FLOAT(K) ) .GE. lmax) THEN
           fault = 1
           RETURN
         END IF
         LCD = LCD * (J / K) 
       END DO
      
      END ! of LCD
      
        

*=======================================================================
*       
*=======================================================================
      PROGRAM gcdlcd
      
       IMPLICIT NONE
       INTEGER N 
       PARAMETER (N=30)
       INTEGER F(N), G(N), fault, GCD, LCD, numer, denom
       
!       DATA F / 80, 64, 16, 40, 32, 8, 16, 8, 16, 4, 8 /
!       DATA G / 3, 3, 1, 3, 3, 1, 3, 3, 9, 3, 9 /
       DATA F / 19664, 15608, 4128, 9832, 2600, 2064, 4912, 3904, 1032,
     &   2456, 1952, 1552, 1232, 976, 776, 616, 488, 128, 364, 80, 64, 
     &   152, 40, 32, 80, 64, 16, 40, 32, 8 /
       DATA G / 3, 3, 1, 3, 1, 1, 3, 3, 1, 3, 3, 3, 3, 3, 2, 3, 3, 1, 
     &   3, 1, 1, 3, 1, 1, 3, 3, 1, 3, 3, 1 /    
       
       numer = GCD(F, N)
       denom = LCD(G, N, fault)
       PRINT *, 'The GCD of the first array is ', numer
       PRINT *, 'The LCD of the second array is ', denom

      END