Deploying to stdlib-fpm from @ 65b74f8 🚀

Author: milancurcic

src/stdlib_math.f90

@@ -5,7 +5,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
@@ -29,6 +29,17 @@ module stdlib_math
         module procedure clip_qp
     end interface clip
 
+    !> Returns the greatest common divisor of two integers
+    !> ([Specification](../page/specs/stdlib_math.html#gcd))
+    !>
+    !> Version: experimental
+    interface gcd
+        module procedure gcd_int8
+        module procedure gcd_int16
+        module procedure gcd_int32
+        module procedure gcd_int64
+    end interface gcd
+
     interface linspace
     !! Version: Experimental
     !!
@@ -648,4 +659,77 @@ elemental function clip_qp(x, xmin, xmax) result(res)
         res = max(min(x, xmax), xmin)
     end function clip_qp
 
+
+    !> Returns the greatest common divisor of two integers of kind int8
+    !> using the Euclidean algorithm.
+    elemental function gcd_int8(a, b) result(res)
+        integer(int8), intent(in) :: a
+        integer(int8), intent(in) :: b
+        integer(int8) :: res
+
+        integer(int8) :: rem, tmp
+
+        rem = min(abs(a), abs(b))
+        res = max(abs(a), abs(b))
+        do while (rem /= 0_int8)
+          tmp = rem
+          rem = mod(res, rem)
+          res = tmp
+        end do
+    end function gcd_int8
+
+    !> Returns the greatest common divisor of two integers of kind int16
+    !> using the Euclidean algorithm.
+    elemental function gcd_int16(a, b) result(res)
+        integer(int16), intent(in) :: a
+        integer(int16), intent(in) :: b
+        integer(int16) :: res
+
+        integer(int16) :: rem, tmp
+
+        rem = min(abs(a), abs(b))
+        res = max(abs(a), abs(b))
+        do while (rem /= 0_int16)
+          tmp = rem
+          rem = mod(res, rem)
+          res = tmp
+        end do
+    end function gcd_int16
+
+    !> Returns the greatest common divisor of two integers of kind int32
+    !> using the Euclidean algorithm.
+    elemental function gcd_int32(a, b) result(res)
+        integer(int32), intent(in) :: a
+        integer(int32), intent(in) :: b
+        integer(int32) :: res
+
+        integer(int32) :: rem, tmp
+
+        rem = min(abs(a), abs(b))
+        res = max(abs(a), abs(b))
+        do while (rem /= 0_int32)
+          tmp = rem
+          rem = mod(res, rem)
+          res = tmp
+        end do
+    end function gcd_int32
+
+    !> Returns the greatest common divisor of two integers of kind int64
+    !> using the Euclidean algorithm.
+    elemental function gcd_int64(a, b) result(res)
+        integer(int64), intent(in) :: a
+        integer(int64), intent(in) :: b
+        integer(int64) :: res
+
+        integer(int64) :: rem, tmp
+
+        rem = min(abs(a), abs(b))
+        res = max(abs(a), abs(b))
+        do while (rem /= 0_int64)
+          tmp = rem
+          rem = mod(res, rem)
+          res = tmp
+        end do
+    end function gcd_int64
+
 end module stdlib_math

test/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