Source code in the fortran programming language

recursive function gcd_rec(u, v) result(gcd)
    integer             :: gcd
    integer, intent(in) :: u, v
    
    if (mod(u, v) /= 0) then
        gcd = gcd_rec(v, mod(u, v))
    else
        gcd = v
    end if
end function gcd_rec


subroutine gcd_iter(value, u, v)
  integer, intent(out) :: value
  integer, intent(inout) :: u, v
  integer :: t

  do while( v /= 0 )
     t = u
     u = v
     v = mod(t, v)
  enddo
  value = abs(u)
end subroutine gcd_iter


function gcd(v, t)
  integer :: gcd
  integer, intent(in) :: v, t
  integer :: c, b, a

  b = t
  a = v
  do
     c = mod(a, b)
     if ( c == 0) exit
     a = b
     b = c
  end do
  gcd = b ! abs(b)
end function gcd


subroutine gcd_bin(value, u, v)
  integer, intent(out) :: value
  integer, intent(inout) :: u, v
  integer :: k, t

  u = abs(u)
  v = abs(v)
  if( u < v ) then
     t = u
     u = v
     v = t
  endif
  if( v == 0 ) then
     value = u
     return
  endif
  k = 1
  do while( (mod(u, 2) == 0).and.(mod(v, 2) == 0) )
     u = u / 2
     v = v / 2
     k = k * 2
  enddo
  if( (mod(u, 2) == 0) ) then
     t = u
  else
     t = -v
  endif
  do while( t /= 0 )
     do while( (mod(t, 2) == 0) )
        t = t / 2
     enddo
     if( t > 0 ) then
        u = t
     else
        v = -t
     endif
     t = u - v
  enddo
  value = u * k
end subroutine gcd_bin


! Stein’s algorithm implemented in Fortran 2008.
! Translated from my implementation for ATS/Postiats.

elemental function gcd (u, v) result (d)
  implicit none
  integer, intent(in) :: u, v
  integer :: d

  integer :: x, y

  ! gcd(x,y) = gcd(u,v), but x and y are non-negative and x <= y.
  x = min (abs (u), abs (v))
  y = max (abs (u), abs (v))

  if (x == 0) then
     d = y
  else
     d = gcd_pos_pos (x, y)
  end if

contains

  elemental function gcd_pos_pos (u, v) result (d)
    integer, intent(in) :: u, v
    integer :: d

    integer :: n
    integer :: x, y
    integer :: p, q

    ! n = the number of common factors of two in u and v.
    n = trailz (ior (u, v))

    ! Remove the common twos from u and v, giving x and y.
    x = ishft (u, -n)
    y = ishft (v, -n)

    ! Make both numbers odd. One of the numbers already was odd.
    ! There is no effect on the value of their gcd.
    x = ishft (x, -trailz (x))
    y = ishft (y, -trailz (y))

    do while (x /= y)
       ! If x > y then swap x and y, renaming them p
       ! and q. Thus p <= q, and gcd(p,q) = gcd(x,y).
       p = min (x, y)
       q = max (x, y)

       x = p                    ! x remains odd.
       q = q - p
       y = ishft (q, -trailz (q)) ! Make y odd again.
    end do

    ! Put the common twos back in.
    d = ishft (x, n)
  end function gcd_pos_pos

end function gcd

program test_gcd
  implicit none

  interface
     elemental function gcd (u, v) result (d)
       integer, intent(in) :: u, v
       integer :: d
     end function gcd
  end interface

  write (*, '("gcd (0, 0) = ", I0)') gcd (0, 0)
  write (*, '("gcd (0, 10) = ", I0)') gcd (0, 10)
  write (*, '("gcd (-6, -9) = ", I0)') gcd (-6, -9)
  write (*, '("gcd (64 * 5, -16 * 3) = ", I0)') gcd (64 * 5, -16 * 3)
  write (*, '("gcd (40902, 24140) = ", I0)') gcd (40902, 24140)
  write (*, '("gcd (-40902, 24140) = ", I0)') gcd (-40902, 24140)
  write (*, '("gcd (40902, -24140) = ", I0)') gcd (40902, -24140)
  write (*, '("gcd (-40902, -24140) = ", I0)') gcd (-40902, -24140)
  write (*, '("gcd (24140, 40902) = ", I0)') gcd (24140, 40902)

end program test_gcd