diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md
index 8d7fb516b..9f7f9d6ce 100644
--- a/doc/specs/stdlib_math.md
+++ b/doc/specs/stdlib_math.md
@@ -91,6 +91,51 @@ program demo_clip_real
 end program demo_clip_real
 ```
 
+### `gcd` function
+
+#### Description
+
+Returns the greatest common divisor of two integers.
+
+#### Syntax
+
+`res = [[stdlib_math(module):gcd(interface)]] (a, b)`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Elemental function.
+
+#### Argument(s)
+
+`a`: One integer with `intent(in)` to get the divisor for.
+`b`: Another integer with `intent(in)` to get the divisor for.
+
+Note: All arguments must be integers of the same `kind`.
+
+#### Output value or Result value
+
+Returns an integer of the same `kind` as that of the arguments.
+
+#### Examples
+
+##### Example 1:
+
+```fortran
+program demo_gcd
+  use stdlib_math, only: gcd
+  implicit none
+  integer :: a, b, c
+  
+  a = 48
+  b = 18
+  c = gcd(a, b) ! returns 6
+end program demo_gcd
+```
+
 ### `linspace` - Create a linearly spaced rank one array
 
 #### Description
@@ -342,4 +387,4 @@ program demo_math_arange
     print *, arange(0.0,2.0,0.0)       !! [0.0,1.0,2.0]. Not recommended: `step` argument is zero!
 
 end program demo_math_arange
-```
\ No newline at end of file
+```
diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp
index 2ca0f543d..3ac37ca78 100644
--- a/src/stdlib_math.fypp
+++ b/src/stdlib_math.fypp
@@ -8,7 +8,7 @@ module stdlib_math
 
     implicit none
     private
-    public :: clip, linspace, logspace
+    public :: clip, gcd, linspace, logspace
     public :: EULERS_NUMBER_SP, EULERS_NUMBER_DP, EULERS_NUMBER_QP
     public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH
     public :: arange
@@ -28,6 +28,16 @@ module stdlib_math
     #:endfor
     end interface clip
 
+    !> Returns the greatest common divisor of two integers
+    !> ([Specification](../page/specs/stdlib_math.html#gcd))
+    !>
+    !> Version: experimental
+    interface gcd
+    #:for k1, t1 in INT_KINDS_TYPES
+        module procedure gcd_${k1}$
+    #:endfor
+    end interface gcd
+
     interface linspace
     !! Version: Experimental
     !!
@@ -292,4 +302,25 @@ contains
     end function clip_${k1}$
 
     #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+    !> Returns the greatest common divisor of two integers of kind ${k1}$
+    !> using the Euclidean algorithm.
+    elemental function gcd_${k1}$(a, b) result(res)
+        ${t1}$, intent(in) :: a
+        ${t1}$, intent(in) :: b
+        ${t1}$ :: res
+
+        ${t1}$ :: rem, tmp
+
+        rem = min(abs(a), abs(b))
+        res = max(abs(a), abs(b))
+        do while (rem /= 0_${k1}$)
+          tmp = rem
+          rem = mod(res, rem)
+          res = tmp
+        end do
+    end function gcd_${k1}$
+
+    #:endfor
 end module stdlib_math
diff --git a/src/tests/math/test_stdlib_math.f90 b/src/tests/math/test_stdlib_math.f90
index 7fafc6bfe..00447d6dd 100644
--- a/src/tests/math/test_stdlib_math.f90
+++ b/src/tests/math/test_stdlib_math.f90
@@ -1,7 +1,7 @@
 ! SPDX-Identifier: MIT
 
 program test_stdlib_math
-    use stdlib_math, only: clip
+    use stdlib_math, only: clip, gcd
     use stdlib_error, only: check
     use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, qp
     implicit none
@@ -95,4 +95,23 @@ program test_stdlib_math
     call check(clip(-55891546.2_qp, -689712245.23_qp, -8958133457.23_qp) == -689712245.23_qp, &
                'clip_qp failed for invalid case', warn=.true.)
 
+
+    ! gcd function
+    ! testing format: check(gcd(a, b) == correct answer)
+    call check(gcd(0, 0) == 0, 'gcd(0, 0) failed.', warn=.true.)
+    call check(gcd(2, 0) == 2, 'gcd(2, 0) failed.', warn=.true.)
+    call check(gcd(0, -2) == 2, 'gcd(0, -2) failed.', warn=.true.)
+    call check(gcd(3, 3) == 3, 'gcd(3, 3) failed.', warn=.true.)
+    call check(gcd(9, 6) == 3, 'gcd(9, 6) failed.', warn=.true.)
+    call check(gcd(6, 9) == 3, 'gcd(6, 9) failed.', warn=.true.)
+    call check(gcd(-9, 6) == 3, 'gcd(-9, 6) failed.', warn=.true.)
+    call check(gcd(9, -6) == 3, 'gcd(9, -6) failed.', warn=.true.)
+    call check(gcd(-9, -6) == 3, 'gcd(-9, -6) failed.', warn=.true.)
+    call check(gcd(97, 91) == 1, 'gcd(97, 91) failed.', warn=.true.)
+
+    call check(gcd(48_int8, 18_int8) == 6_int8, 'gcd(48, 18) failed for int8.', warn=.true.)
+    call check(gcd(48_int16, 18_int16) == 6_int16, 'gcd(48, 18) failed for int16', warn=.true.)
+    call check(gcd(48_int32, 18_int32) == 6_int32, 'gcd(48, 18) failed for int32', warn=.true.)
+    call check(gcd(48_int64, 18_int64) == 6_int64, 'gcd(48, 18) failed for int64', warn=.true.)
+
 end program test_stdlib_math