module math_module
! This file:
! http://ftp.aset.psu.edu/pub/ger/fortran/hdk/MathModuleDemo.f95
!
! F-compiler compatible; example adapted from Walt Brainerd,
! F-Compiler distribution, example, math_module.f95.
! Date: 5 December 2003

   integer, parameter, private :: double = &
      selected_real_kind(precision(0.0)+1)
!     selected_real_kind(15)
! integer, parameter, public :: int_kind = kind(0)
  integer, parameter, public :: int_kind = &
      selected_int_kind(9)

   public :: gcd

   real, parameter, public :: pi =  &
      3.1415926535897932384626433832795028841972

   real, parameter, public :: e =  &
      2.7182818284590452353602874713526624977572

   real, parameter, public :: gamma =  &
      0.5772156649015328606065120900824024310422

   real, parameter, public :: phi =  &
      1.6180339887498948482045868343656381177203

   real(kind=double), parameter, public :: pi_double =  &
      3.1415926535897932384626433832795028841972_double

   real(kind=double), parameter, public :: e_double =  &
      2.7182818284590452353602874713526624977572_double

   real(kind=double), parameter, public :: gamma_double =  &
      0.5772156649015328606065120900824024310422_double

   real(kind=double), parameter, public :: phi_double =  &
      1.6180339887498948482045868343656381177203_double

contains

  elemental function gcd( a, b) result(gcd_result)

    integer( kind= int_kind) :: gcd_result
    integer( kind= int_kind), intent( in) :: a, b
    integer( kind= int_kind) :: a_gcd, b_gcd, rnp1, rn, rnm1

! if a or b zero, abs( other) is gcd

    if( a == 0_int_kind )then
      gcd_result = abs( b)                     ! note this might return zero
      return                                   ! done
    endif

    if( b == 0_int_kind )then
      gcd_result = abs( a)
      return                                   ! done
    endif

! set |a| >= |b| ( > 0)
! r1 = a .mod. b
! r0 = b

    a_gcd = max( abs( a), abs( b))              ! a is the greater
    b_gcd = min( abs( a), abs( b))              ! b is the lesser
    rn = modulo(a_gcd, b_gcd)                   ! n = 1
    rnm1 = b_gcd                                ! n-1 = 0

! while rn /= 0
!    compute rn+1 = rn .mod. rn-1
! gcd() = rnm1

    do
      if( rn == 0_int_kind) then
        exit
      endif
      rnp1 = modulo(rnm1, rn)                    ! rn+1
      rnm1 = rn                                  ! n --> n-1
      rn = rnp1                                  ! n+1 --> n

    enddo                                        ! til zero remainder

    gcd_result = rnm1                            ! done

  end function gcd

end module math_module

program testmath
  use math_module
  print *,"DP pi=",pi_double,", SP e=",e
  print *,"GCD(15,21)=",gcd(15,21)
  print *,"DP Huge(e)=",huge(e_double)
end program testmath