diff --git a/.github/workflows/fpm-deployment.yml b/.github/workflows/fpm-deployment.yml
index 2e4a3203c..19cf58d8b 100644
--- a/.github/workflows/fpm-deployment.yml
+++ b/.github/workflows/fpm-deployment.yml
@@ -35,9 +35,12 @@ jobs:
         with:
           fpm-version: 'v0.10.0'
 
-      - run: |
+      - run: | # Just for deployment: create stdlib-fpm folder
           python config/fypp_deployment.py --deploy_stdlib_fpm
-          fpm test --profile release
+
+      - run: | # Use fpm gnu ci to check xdp and qp 
+          python config/fypp_deployment.py --with_xdp --with_qp
+          fpm test --profile release --flag '-DWITH_XDP -DWITH_QP'
 
     # Update and deploy the f90 files generated by github-ci to the `stdlib-fpm` branch.
       - name: Deploy 🚀
diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp
index 2307a6660..e91f14954 100644
--- a/src/stdlib_specialfunctions_gamma.fypp
+++ b/src/stdlib_specialfunctions_gamma.fypp
@@ -1,7 +1,7 @@
-#:set WITH_QP = False
-#:set WITH_XDP = False
 #:include "common.fypp"
-#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES
+#:set R_KINDS_TYPES = [KT for KT in REAL_KINDS_TYPES if KT[0] in ["sp","dp"]]
+#:set C_KINDS_TYPES = [KT for KT in CMPLX_KINDS_TYPES if KT[0] in ["sp","dp"]]
+#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES
 module stdlib_specialfunctions_gamma
     use iso_fortran_env, only : qp => real128
     use stdlib_kinds, only :  sp, dp, int8, int16, int32, int64
@@ -15,7 +15,7 @@ module stdlib_specialfunctions_gamma
     integer(int32), parameter :: max_fact_int32 = 13_int32
     integer(int64), parameter :: max_fact_int64 = 21_int64
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     ${t1}$, parameter :: tol_${k1}$ = epsilon(1.0_${k1}$)
     #:endfor
     real(qp), parameter :: tol_qp = epsilon(1.0_qp)
@@ -63,12 +63,12 @@ module stdlib_specialfunctions_gamma
     !! Lower incomplete gamma function
     !!
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure ingamma_low_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure ingamma_low_${t1[0]}$${k1}$
         #:endfor
     end interface lower_incomplete_gamma
@@ -79,12 +79,12 @@ module stdlib_specialfunctions_gamma
     !! Logarithm of lower incomplete gamma function
     !!
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure l_ingamma_low_${t1[0]}$${k1}$
         #:endfor
     end interface log_lower_incomplete_gamma
@@ -95,12 +95,12 @@ module stdlib_specialfunctions_gamma
     !! Upper incomplete gamma function
     !!
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure ingamma_up_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure ingamma_up_${t1[0]}$${k1}$
         #:endfor
     end interface upper_incomplete_gamma
@@ -111,12 +111,12 @@ module stdlib_specialfunctions_gamma
     !! Logarithm of upper incomplete gamma function
     !!
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure l_ingamma_up_${t1[0]}$${k1}$
         #:endfor
     end interface log_upper_incomplete_gamma
@@ -127,12 +127,12 @@ module stdlib_specialfunctions_gamma
     !! Regularized (normalized) lower incomplete gamma function, P
     !!
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure regamma_p_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure regamma_p_${t1[0]}$${k1}$
         #:endfor
     end interface regularized_gamma_p
@@ -143,12 +143,12 @@ module stdlib_specialfunctions_gamma
     !! Regularized (normalized) upper incomplete gamma function, Q
     !!
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure regamma_q_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure regamma_q_${t1[0]}$${k1}$
         #:endfor
     end interface regularized_gamma_q
@@ -159,12 +159,12 @@ module stdlib_specialfunctions_gamma
     ! Incomplete gamma G function.
     ! Internal use only
     !
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
         module procedure gpx_${t1[0]}$${k1}$         !for real p and x
         #:endfor
 
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
         module procedure gpx_${t1[0]}$${k1}$${k2}$   !for integer p and real x
           #:endfor
         #:endfor
@@ -177,7 +177,7 @@ module stdlib_specialfunctions_gamma
     ! Internal use only
     !
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
           module procedure l_gamma_${t1[0]}$${k1}$${k2}$
           #:endfor
         #:endfor
@@ -218,7 +218,7 @@ contains
 
 
 
-    #:for k1, t1 in CMPLX_KINDS_TYPES
+    #:for k1, t1 in C_KINDS_TYPES
       #:if k1 == "sp"
           #:set k2 = "dp"
       #:elif k1 == "dp"
@@ -373,7 +373,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
 
     impure elemental function l_gamma_${t1[0]}$${k1}$${k2}$(z, x) result(res)
     !
@@ -414,7 +414,7 @@ contains
 
 
 
-    #:for k1, t1 in CMPLX_KINDS_TYPES
+    #:for k1, t1 in C_KINDS_TYPES
       #:if k1 == "sp"
           #:set k2 = "dp"
       #:elif k1 == "dp"
@@ -556,7 +556,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
       #:if k1 == "sp"
           #:set k2 = "dp"
       #:elif k1 == "dp"
@@ -702,7 +702,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res)
     !
     ! Approximation of incomplete gamma G function with integer argument p.
@@ -841,7 +841,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
     !
     ! Approximation of lower incomplete gamma function with real p.
@@ -878,7 +878,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x)          &
         result(res)
     !
@@ -918,7 +918,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     impure elemental function l_ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
 
         ${t1}$, intent(in) :: p, x
@@ -955,7 +955,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x)        &
         result(res)
 
@@ -987,7 +987,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
     !
     ! Approximation of upper incomplete gamma function with real p.
@@ -1025,7 +1025,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x)           &
         result(res)
     !
@@ -1067,7 +1067,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     impure elemental function l_ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
 
         ${t1}$, intent(in) :: p, x
@@ -1105,7 +1105,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x)         &
         result(res)
 
@@ -1146,7 +1146,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res)
     !
     ! Approximation of regularized incomplete gamma function P(p,x) for real p
@@ -1181,7 +1181,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(p, x) result(res)
     !
     ! Approximation of regularized incomplete gamma function P(p,x) for integer p
@@ -1217,7 +1217,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res)
     !
     ! Approximation of regularized incomplete gamma function Q(p,x) for real p
@@ -1252,7 +1252,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
     impure elemental function regamma_q_${t1[0]}$${k1}$${k2}$(p, x) result(res)
     !
     ! Approximation of regularized incomplet gamma function Q(p,x) for integer p
diff --git a/test/quadrature/test_simps.fypp b/test/quadrature/test_simps.fypp
index 26814c1c1..cb6272612 100644
--- a/test/quadrature/test_simps.fypp
+++ b/test/quadrature/test_simps.fypp
@@ -7,15 +7,9 @@ module test_simps
 
     implicit none
 
-    real(sp), parameter :: tol_sp = 1000 * epsilon(1.0_sp)
-    real(dp), parameter :: tol_dp = 1000 * epsilon(1.0_dp)
-#:if WITH_XDP
-    real(xdp), parameter :: tol_xdp = 1000 * epsilon(1.0_xdp)
-#:endif
-#:if WITH_QP
-    real(qp), parameter :: tol_qp = 1000 * epsilon(1.0_qp)
-#:endif
-
+    #:for k1, t1 in REAL_KINDS_TYPES
+    ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$)
+    #:endfor
 
 contains
 
@@ -25,135 +19,72 @@ contains
         type(unittest_type), allocatable, intent(out) :: testsuite(:)
 
         testsuite = [ &
-            new_unittest("simps_sp", test_simps_sp), &
-            new_unittest("simps_dp", test_simps_dp), &
-            new_unittest("simps_qp", test_simps_qp), &
-            new_unittest("simps_weights_sp", test_simps_weights_sp), &
-            new_unittest("simps_weights_dp", test_simps_weights_dp), &
-            new_unittest("simps_weights_qp", test_simps_weights_qp), &
-            new_unittest("simps_zero_sp", test_simps_zero_sp), &
-            new_unittest("simps_zero_dp", test_simps_zero_dp), &
-            new_unittest("simps_zero_qp", test_simps_zero_qp), &
-            new_unittest("simps_even_sp", test_simps_even_sp), &
-            new_unittest("simps_even_dp", test_simps_even_dp), &
-            new_unittest("simps_even_qp", test_simps_even_qp), &
-            new_unittest("simps_weights_even_sp", test_simps_weights_even_sp), &
-            new_unittest("simps_weights_even_dp", test_simps_weights_even_dp), &
-            new_unittest("simps_weights_even_qp", test_simps_weights_even_qp), &
-            new_unittest("simps_six_sp", test_simps_six_sp), &
-            new_unittest("simps_six_dp", test_simps_six_dp), &
-            new_unittest("simps_six_qp", test_simps_six_qp) &
+            #:for k1, t1 in REAL_KINDS_TYPES[0:1] # set the first test independently to initialize the table
+            new_unittest("simps_${k1}$", test_simps_sp) &
+            #:endfor
+            #:for k1, t1 in REAL_KINDS_TYPES[1:]
+            , new_unittest("simps_${k1}$", test_simps_${k1}$) &
+            #:endfor
+            #:for k1, t1 in REAL_KINDS_TYPES
+            , new_unittest("simps_weights_${k1}$", test_simps_weights_${k1}$) &
+            #:endfor
+            #:for k1, t1 in REAL_KINDS_TYPES
+            , new_unittest("simps_zero_${k1}$", test_simps_zero_${k1}$) &
+            #:endfor
+            #:for k1, t1 in REAL_KINDS_TYPES
+            , new_unittest("simps_even_${k1}$", test_simps_even_${k1}$) &
+            #:endfor
+            #:for k1, t1 in REAL_KINDS_TYPES
+            , new_unittest("simps_weights_even_${k1}$", test_simps_weights_even_${k1}$) &
+            #:endfor
+            #:for k1, t1 in REAL_KINDS_TYPES
+            , new_unittest("simps_six_${k1}$", test_simps_six_${k1}$) &
+            #:endfor
             ]
     end subroutine collect_simps
 
-    subroutine test_simps_sp(error)
+    #:for k1, t1 in REAL_KINDS_TYPES
+    subroutine test_simps_${k1}$(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
 
         integer, parameter :: n = 13
-        real(sp), dimension(n) :: y
-        real(sp), dimension(n) :: x
-        real(sp) :: val
-        real(sp) :: ans
+        ${t1}$ :: y(n)
+        ${t1}$ :: x(n)
+        ${t1}$ :: val
+        ${t1}$ :: ans
         integer :: i
 
 
         y = [(real(i-1, sp)**2, i = 1, n)]
 
-        val = simps(y, 1.0_sp)
-        ans = 576.0_sp
-        call check(error, val, ans, thr=tol_sp)
-        if (allocated(error)) return
-
-        val = simps(y, 0.5_sp)
-        ans = 288.0_sp
-        call check(error, val, ans, thr=tol_sp)
-        if (allocated(error)) return
-
-        x = [(0.25_sp*(i-1), i = 1, n)]
-        val = simps(y, x)
-        ans = 144.0_sp
-        call check(error, val, ans, thr=tol_sp)
-    end subroutine test_simps_sp
-
-
-    subroutine test_simps_dp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 13
-        real(dp), dimension(n) :: y
-        real(dp), dimension(n) :: x
-        real(dp) :: val
-        real(dp) :: ans
-        integer :: i
-
-
-        y = [(real(i-1, dp)**2, i = 1, n)]
-
-        val = simps(y, 1.0_dp)
-        ans = 576.0_dp
-        call check(error, val, ans, thr=tol_dp)
-        if (allocated(error)) return
-
-        val = simps(y, 0.5_dp)
-        ans = 288.0_dp
-        call check(error, val, ans, thr=tol_dp)
-        if (allocated(error)) return
-
-        x = [(0.25_dp*(i-1), i = 1, n)]
-        val = simps(y, x)
-        ans = 144.0_dp
-        call check(error, val, ans, thr=tol_dp)
-    end subroutine test_simps_dp
-
-
-    subroutine test_simps_qp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-#:if WITH_QP
-        integer, parameter :: n = 13
-        real(qp), dimension(n) :: y
-        real(qp), dimension(n) :: x
-        real(qp) :: val
-        real(qp) :: ans
-        integer :: i
-
-
-        y = [(real(i-1, qp)**2, i = 1, n)]
-
-        val = simps(y, 1.0_qp)
-        ans = 576.0_qp
-        call check(error, val, ans, thr=tol_qp)
+        val = simps(y, 1.0_${k1}$)
+        ans = 576.0_${k1}$
+        call check(error, val, ans, thr=tol_${k1}$)
         if (allocated(error)) return
 
-        val = simps(y, 0.5_qp)
-        ans = 288.0_qp
-        call check(error, val, ans, thr=tol_qp)
+        val = simps(y, 0.5_${k1}$)
+        ans = 288.0_${k1}$
+        call check(error, val, ans, thr=tol_${k1}$)
         if (allocated(error)) return
 
-        x = [(0.25_qp*(i-1), i = 1, n)]
+        x = [(0.25_${k1}$*(i-1), i = 1, n)]
         val = simps(y, x)
-        ans = 144.0_qp
-        call check(error, val, ans, thr=tol_qp)
-#:else
-        call skip_test(error, "Quadruple precision is not enabled")
-#:endif
-    end subroutine test_simps_qp
-
+        ans = 144.0_${k1}$
+        call check(error, val, ans, thr=tol_${k1}$)
+    end subroutine test_simps_${k1}$
 
-    subroutine test_simps_weights_sp(error)
+    subroutine test_simps_weights_${k1}$(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
 
         integer, parameter :: n = 17
-        real(sp), dimension(n) :: y
-        real(sp), dimension(n) :: x
-        real(sp), dimension(n) :: w
+        ${t1}$ :: y(n)
+        ${t1}$ :: x(n)
+        ${t1}$ :: w(n)
         integer :: i
-        real(sp) :: val
-        real(sp) :: ans
+        ${t1}$ :: val
+        ${t1}$ :: ans
 
 
         y = [(real(i-1, sp), i = 1, n)]
@@ -162,233 +93,71 @@ contains
         w = simps_weights(x)
         val = sum(w*y)
         ans = simps(y, x)
-        call check(error, val, ans, thr=tol_sp)
-    end subroutine test_simps_weights_sp
-
-
-    subroutine test_simps_weights_dp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 17
-        real(dp), dimension(n) :: y
-        real(dp), dimension(n) :: x
-        real(dp), dimension(n) :: w
-        integer :: i
-        real(dp) :: val
-        real(dp) :: ans
-
-
-        y = [(real(i-1, dp), i = 1, n)]
-
-        x = y
-        w = simps_weights(x)
-        val = sum(w*y)
-        ans = simps(y, x)
-        call check(error, val, ans, thr=tol_dp)
-    end subroutine test_simps_weights_dp
-
-
-    subroutine test_simps_weights_qp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-#:if WITH_QP
-        integer, parameter :: n = 17
-        real(qp), dimension(n) :: y
-        real(qp), dimension(n) :: x
-        real(qp), dimension(n) :: w
-        integer :: i
-        real(qp) :: val
-        real(qp) :: ans
-
-
-        y = [(real(i-1, qp), i = 1, n)]
-
-        x = y
-        w = simps_weights(x)
-        val = sum(w*y)
-        ans = simps(y, x)
-        call check(error, val, ans, thr=tol_qp)
-#:else
-        call skip_test(error, "Quadruple precision is not enabled")
-#:endif
-    end subroutine test_simps_weights_qp
-
-
-    subroutine test_simps_zero_sp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        real(sp), dimension(0) :: a
-
-
-        call check(error, abs(simps(a, 1.0_sp)) < epsilon(0.0_sp))
-        if (allocated(error)) return
-        call check(error, abs(simps([1.0_sp], 1.0_sp)) < epsilon(0.0_sp))
-        if (allocated(error)) return
-        call check(error, abs(simps(a, a)) < epsilon(0.0_sp))
-        if (allocated(error)) return
-        call check(error, abs(simps([1.0_sp], [1.0_sp])) < epsilon(0.0_sp))
-    end subroutine test_simps_zero_sp
-
-
-    subroutine test_simps_zero_dp(error)
+        call check(error, val, ans, thr=tol_${k1}$)
+    end subroutine test_simps_weights_${k1}$
+    
+    subroutine test_simps_zero_${k1}$(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
 
-        real(dp), dimension(0) :: a
+        ${t1}$, dimension(0) :: a
 
 
-        call check(error, abs(simps(a, 1.0_dp)) < epsilon(0.0_dp))
+        call check(error, abs(simps(a, 1.0_${k1}$)) < epsilon(0.0_${k1}$))
         if (allocated(error)) return
-        call check(error, abs(simps([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
-        if (allocated(error)) return
-        call check(error, abs(simps(a, a)) < epsilon(0.0_dp))
-        if (allocated(error)) return
-        call check(error, abs(simps([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
-    end subroutine test_simps_zero_dp
-
-
-    subroutine test_simps_zero_qp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-#:if WITH_QP
-        real(qp), dimension(0) :: a
-
-
-        call check(error, abs(simps(a, 1.0_qp)) < epsilon(0.0_qp))
-        if (allocated(error)) return
-        call check(error, abs(simps([1.0_qp], 1.0_qp)) < epsilon(0.0_qp))
-        if (allocated(error)) return
-        call check(error, abs(simps(a, a)) < epsilon(0.0_qp))
-        if (allocated(error)) return
-        call check(error, abs(simps([1.0_qp], [1.0_qp])) < epsilon(0.0_qp))
-#:else
-        call skip_test(error, "Quadruple precision is not enabled")
-#:endif
-    end subroutine test_simps_zero_qp
-
-
-    subroutine test_simps_even_sp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 11
-        real(sp), dimension(n) :: y
-        real(sp), dimension(n) :: x
-        real(sp) :: val
-        real(sp) :: ans
-        integer :: i
-        integer :: even
-
-
-        y = [(3.0_sp*real(i-1, sp)**2, i = 1, n)]
-
-        do even = -1, 1
-
-            val = simps(y, 1.0_sp)
-            ans = 1000.0_sp
-            call check(error, val, ans, thr=tol_sp)
-            if (allocated(error)) return
-
-            val = simps(y, 0.5_sp)
-            ans = 500.0_sp
-            call check(error, val, ans, thr=tol_sp)
-            if (allocated(error)) return
-
-            x = [(0.25_sp*(i-1), i = 1, n)]
-            val = simps(y, x)
-            ans = 250.0_sp
-            call check(error, val, ans, thr=tol_sp)
-            if (allocated(error)) return
-        end do
-    end subroutine test_simps_even_sp
-
-
-    subroutine test_simps_even_dp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 11
-        real(dp), dimension(n) :: y
-        real(dp), dimension(n) :: x
-        real(dp) :: val
-        real(dp) :: ans
-        integer :: i
-
-
-        y = [(3.0_dp*real(i-1, dp)**2, i = 1, n)]
-
-        val = simps(y, 1.0_dp)
-        ans = 1000.0_dp
-        call check(error, val, ans, thr=tol_dp)
+        call check(error, abs(simps([1.0_${k1}$], 1.0_${k1}$)) < epsilon(0.0_${k1}$))
         if (allocated(error)) return
-
-        val = simps(y, 0.5_dp)
-        ans = 500.0_dp
-        call check(error, val, ans, thr=tol_dp)
+        call check(error, abs(simps(a, a)) < epsilon(0.0_${k1}$))
         if (allocated(error)) return
+        call check(error, abs(simps([1.0_${k1}$], [1.0_${k1}$])) < epsilon(0.0_${k1}$))
+    end subroutine test_simps_zero_${k1}$
 
-        x = [(0.25_dp*(i-1), i = 1, n)]
-        val = simps(y, x)
-        ans = 250.0_dp
-        call check(error, val, ans, thr=tol_dp)
-    end subroutine test_simps_even_dp
-
-
-    subroutine test_simps_even_qp(error)
+    subroutine test_simps_even_${k1}$(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
 
-#:if WITH_QP
         integer, parameter :: n = 11
-        real(qp), dimension(n) :: y
-        real(qp), dimension(n) :: x
-        real(qp) :: val
-        real(qp) :: ans
+        ${t1}$ :: y(n)
+        ${t1}$ :: x(n)
+        ${t1}$ :: val
+        ${t1}$ :: ans
         integer :: i
         integer :: even
 
 
-        y = [(3.0_qp*real(i-1, qp)**2, i = 1, n)]
+        y = [(3.0_${k1}$*real(i-1, sp)**2, i = 1, n)]
 
         do even = -1, 1
 
-            val = simps(y, 1.0_qp)
-            ans = 1000.0_qp
-            call check(error, val, ans, thr=tol_qp)
+            val = simps(y, 1.0_${k1}$)
+            ans = 1000.0_${k1}$
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
 
-            val = simps(y, 0.5_qp)
-            ans = 500.0_qp
-            call check(error, val, ans, thr=tol_qp)
+            val = simps(y, 0.5_${k1}$)
+            ans = 500.0_${k1}$
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
 
-            x = [(0.25_qp*(i-1), i = 1, n)]
+            x = [(0.25_${k1}$*(i-1), i = 1, n)]
             val = simps(y, x)
-            ans = 250.0_qp
-            call check(error, val, ans, thr=tol_qp)
+            ans = 250.0_${k1}$
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
         end do
-#:else
-        call skip_test(error, "Quadruple precision is not enabled")
-#:endif
-    end subroutine test_simps_even_qp
-
+    end subroutine test_simps_even_${k1}$
 
-    subroutine test_simps_weights_even_sp(error)
+    subroutine test_simps_weights_even_${k1}$(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
 
         integer, parameter :: n = 16
-        real(sp), dimension(n) :: y
-        real(sp), dimension(n) :: x
-        real(sp), dimension(n) :: w
+        ${t1}$ :: y(n)
+        ${t1}$ :: x(n)
+        ${t1}$ :: w(n)
         integer :: i
-        real(sp) :: val
-        real(sp) :: ans
+        ${t1}$ :: val
+        ${t1}$ :: ans
         integer :: even
 
 
@@ -399,176 +168,47 @@ contains
             w = simps_weights(x)
             val = sum(w*y)
             ans = simps(y, x)
-            call check(error, val, ans, thr=tol_sp)
-            if (allocated(error)) return
-        end do
-    end subroutine test_simps_weights_even_sp
-
-
-    subroutine test_simps_weights_even_dp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 16
-        real(dp), dimension(n) :: y
-        real(dp), dimension(n) :: x
-        real(dp), dimension(n) :: w
-        integer :: i
-        real(dp) :: val
-        real(dp) :: ans
-        integer :: even
-
-
-        y = [(real(i-1, dp), i = 1, n)]
-        x = y
-
-        do even = -1, 1
-            w = simps_weights(x)
-            val = sum(w*y)
-            ans = simps(y, x)
-            call check(error, val, ans, thr=tol_dp)
-            if (allocated(error)) return
-        end do
-    end subroutine test_simps_weights_even_dp
-
-
-    subroutine test_simps_weights_even_qp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-#:if WITH_QP
-        integer, parameter :: n = 16
-        real(qp), dimension(n) :: y
-        real(qp), dimension(n) :: x
-        real(qp), dimension(n) :: w
-        integer :: i
-        real(qp) :: val
-        real(qp) :: ans
-        integer :: even
-
-
-        y = [(real(i-1, qp), i = 1, n)]
-
-        x = y
-
-        do even = -1, 1
-            w = simps_weights(x)
-            val = sum(w*y)
-            ans = simps(y, x)
-            call check(error, val, ans, thr=tol_qp)
-            if (allocated(error)) return
-        end do
-#:else
-        call skip_test(error, "Quadruple precision is not enabled")
-#:endif
-    end subroutine test_simps_weights_even_qp
-
-
-    subroutine test_simps_six_sp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 6
-        real(sp), dimension(n) :: y
-        real(sp), dimension(n) :: x
-        real(sp) :: val
-        real(sp) :: ans
-        integer :: i
-        integer :: even
-
-
-        y = [(3.0_sp*real(i-1, sp)**2, i = 1, n)]
-
-        do even = -1, 1
-
-            val = simps(y, 1.0_sp)
-            ans = 125.0_sp
-            call check(error, val, ans, thr=tol_sp)
-            if (allocated(error)) return
-
-            val = simps(y, 0.5_sp)
-            ans = 62.5_sp
-            call check(error, val, ans, thr=tol_sp)
-            if (allocated(error)) return
-
-            x = [(0.25_sp*(i-1), i = 1, n)]
-            val = simps(y, x)
-            ans = 31.25_sp
-            call check(error, val, ans, thr=tol_sp)
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
         end do
-    end subroutine test_simps_six_sp
-
-
-    subroutine test_simps_six_dp(error)
-        !> Error handling
-        type(error_type), allocatable, intent(out) :: error
-
-        integer, parameter :: n = 6
-        real(dp), dimension(n) :: y
-        real(dp), dimension(n) :: x
-        real(dp) :: val
-        real(dp) :: ans
-        integer :: i
-
-
-        y = [(3.0_dp*real(i-1, dp)**2, i = 1, n)]
-
-        val = simps(y, 1.0_dp)
-        ans = 125.0_dp
-        call check(error, val, ans, thr=tol_dp)
-        if (allocated(error)) return
-
-        val = simps(y, 0.5_dp)
-        ans = 62.5_dp
-        call check(error, val, ans, thr=tol_dp)
-        if (allocated(error)) return
-
-        x = [(0.25_dp*(i-1), i = 1, n)]
-        val = simps(y, x)
-        ans = 31.25_dp
-        call check(error, val, ans, thr=tol_dp)
-    end subroutine test_simps_six_dp
-
+    end subroutine test_simps_weights_even_${k1}$
 
-    subroutine test_simps_six_qp(error)
+    subroutine test_simps_six_${k1}$(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
 
-#:if WITH_QP
         integer, parameter :: n = 6
-        real(qp), dimension(n) :: y
-        real(qp), dimension(n) :: x
-        real(qp) :: val
-        real(qp) :: ans
+        ${t1}$ :: y(n)
+        ${t1}$ :: x(n)
+        ${t1}$ :: val
+        ${t1}$ :: ans
         integer :: i
         integer :: even
 
 
-        y = [(3.0_qp*real(i-1, qp)**2, i = 1, n)]
+        y = [(3.0_${k1}$*real(i-1, sp)**2, i = 1, n)]
 
         do even = -1, 1
 
-            val = simps(y, 1.0_qp)
-            ans = 125.0_qp
-            call check(error, val, ans, thr=tol_qp)
+            val = simps(y, 1.0_${k1}$)
+            ans = 125.0_${k1}$
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
 
-            val = simps(y, 0.5_qp)
-            ans = 62.5_qp
-            call check(error, val, ans, thr=tol_qp)
+            val = simps(y, 0.5_${k1}$)
+            ans = 62.5_${k1}$
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
 
-            x = [(0.25_qp*(i-1), i = 1, n)]
+            x = [(0.25_${k1}$*(i-1), i = 1, n)]
             val = simps(y, x)
-            ans = 31.25_qp
-            call check(error, val, ans, thr=tol_qp)
+            ans = 31.25_${k1}$
+            call check(error, val, ans, thr=tol_${k1}$)
             if (allocated(error)) return
         end do
-#:else
-        call skip_test(error, "Quadruple precision is not enabled")
-#:endif
-    end subroutine test_simps_six_qp
+    end subroutine test_simps_six_${k1}$
+    
+    #:endfor
 
 end module
 
diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp
index a5853df1f..62ee4c1f9 100644
--- a/test/specialfunctions/test_specialfunctions_gamma.fypp
+++ b/test/specialfunctions/test_specialfunctions_gamma.fypp
@@ -1,7 +1,7 @@
-#:set WITH_QP = False
-#:set WITH_XDP = False
 #:include "common.fypp"
-#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES
+#:set R_KINDS_TYPES = [KT for KT in REAL_KINDS_TYPES if KT[0] in ["sp","dp"]]
+#:set C_KINDS_TYPES = [KT for KT in CMPLX_KINDS_TYPES if KT[0] in ["sp","dp"]]
+#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES
 module test_specialfunctions_gamma
     use testdrive, only : new_unittest, unittest_type, error_type, check
     use stdlib_kinds, only: sp, dp, int8, int16, int32, int64
@@ -18,13 +18,10 @@ module test_specialfunctions_gamma
 
     public :: collect_specialfunctions_gamma
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
     ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$)
     #:endfor
 
-
-
-
 contains
 
     subroutine collect_specialfunctions_gamma(testsuite)
@@ -46,7 +43,7 @@ contains
         #:endfor
 
         #:for k1, t1 in INT_KINDS_TYPES
-          #:for k2, t2 in REAL_KINDS_TYPES
+          #:for k2, t2 in R_KINDS_TYPES
             , new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$",     &
                            test_lincgamma_${t1[0]}$${k1}$${k2}$)               &
             , new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", &
@@ -62,7 +59,7 @@ contains
           #:endfor
         #:endfor
 
-        #:for k1, t1 in REAL_KINDS_TYPES
+        #:for k1, t1 in R_KINDS_TYPES
             , new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$",           &
                            test_lincgamma_${t1[0]}$${k1}$)                     &
             , new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$",       &
@@ -268,7 +265,7 @@ contains
 
 
     #:for k1, t1 in INT_KINDS_TYPES
-      #:for k2, t2 in REAL_KINDS_TYPES
+      #:for k2, t2 in R_KINDS_TYPES
 
     subroutine test_lincgamma_${t1[0]}$${k1}$${k2}$(error)
         type(error_type), allocatable, intent(out) :: error
@@ -417,7 +414,7 @@ contains
 
 
 
-    #:for k1, t1 in REAL_KINDS_TYPES
+    #:for k1, t1 in R_KINDS_TYPES
 
     subroutine test_lincgamma_${t1[0]}$${k1}$(error)
         type(error_type), allocatable, intent(out) :: error