diff --git a/doc/specs/index.md b/doc/specs/index.md
index 1378fa8b8..e78f95b44 100644
--- a/doc/specs/index.md
+++ b/doc/specs/index.md
@@ -32,6 +32,7 @@ This is an index/directory of the specifications (specs) for each new module/fea
  - [stats_distributions_uniform](./stdlib_stats_distribution_uniform.html) - Uniform Probability Distribution
  - [stats_distributions_normal](./stdlib_stats_distribution_normal.html) - Normal Probability Distribution
  - [stats_distributions_exponential](./stdlib_stats_distribution_exponential.html) - Exponential Probability Distribution
+ - [stats_distributions_gamma](./stdlib_stats_distribution_gamma.html) - Gamma Probability Distribution
  - [string\_type](./stdlib_string_type.html) - Basic string support
  - [stringlist_type](./stdlib_stringlist_type.html) - 1-Dimensional list of strings
  - [strings](./stdlib_strings.html) - String handling and manipulation routines
diff --git a/doc/specs/stdlib_stats_distribution_gamma.md b/doc/specs/stdlib_stats_distribution_gamma.md
new file mode 100644
index 000000000..d3f5939f0
--- /dev/null
+++ b/doc/specs/stdlib_stats_distribution_gamma.md
@@ -0,0 +1,267 @@
+---
+
+title: stats_distribution_gamma
+---
+
+# Statistical Distributions -- Gamma Distribution Module
+
+[TOC]
+
+## `rvs_gamma` - gamma distribution random variates
+
+### Status
+
+Experimental
+
+### Description
+
+With one argument for shape parameter, the function returns a random sample from the standard gamma distribution `Gam(shape)` with `rate = 1.0`. 
+
+With two arguments, the function returns a random sample from gamma distribution `Gam(shape, rate)`.
+
+With three arguments, the function returns a rank one array of gamma distributed random variates.
+
+For complex shape and rate parameters, the real and imaginary parts are sampled independently of each other.
+
+
+### Syntax
+
+`result = [[stdlib_stats_distribution_gamma(module):rvs_gamma(interface)]](shape [, rate] [[, array_size]])`
+
+### Class
+
+Function
+
+### Arguments
+
+`shape` : has `intent(in)` and is a scalar of type `real` or `complex`.
+
+`rate`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
+
+`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind.
+
+### Return value
+
+The result is a scalar or rank one array, with a size of `array_size`, and has the same type of `shape`.
+
+### Example
+
+```fortran
+program demo_gamma_rvs
+    use stdlib_random, only : random_seed
+    use stdlib_stats_distribution_gamma, only: rgamma => rvs_gamma
+
+    implicit none
+    real ::  g(2,3,4)
+    complex :: shape, scale
+    integer :: put, get
+
+    put = 1234567
+    call random_seed(put, get)
+
+    print *, rgamma(2.0)
+    !single standard gamma random variate with shape of 2.0, rate=1.0
+
+! 2.50538206
+
+    print *, rgamma(3.0,2.0)  !gamma random variate with shape=3.0, rate=2.0
+
+! 1.30591583
+
+    g(:,:,:) = 0.5
+    print *, rgamma(g)
+    !a rank 3 array of 60 standard gamma random variates with shape=0.5
+
+!  1.03841162  1.33044529  0.912742674  0.131288037  0.638593793 
+!  1.03565669E-02  0.624804378  1.12179172  4.91380468E-02  6.69969944E-03 
+!  6.67014271E-02  0.132111162  0.101102419  0.648416579  1.14922595 
+!  2.29003578E-02  1.85964716E-04  1.21213868E-02  1.69112933 
+!  7.30440915E-02  0.395139128  0.182758048  0.427981257  0.985665262
+
+    print *, rgamma(0.5,1.0,10)
+    ! an array of 10 random variates with shape=0.5, rate=1.0
+
+!  1.39297554E-04  0.296419382  0.352113068  2.80515051  3.65264394E-04 
+!  0.197743446  5.54569438E-02  9.30598825E-02  1.02596343  1.85311246
+
+    shape = (3.0, 4.0)
+    scale = (2.0, 0.7)
+    print *, rgamma(shape,scale)
+    !single complex gamma random variate with real part of shape = 3.0,
+    !rate=2.0; imagainary part of shape=4.0, rate=0.7
+
+! (0.826188326,3.54749799)
+
+end program demo_gamma_rvs
+```
+
+## `pdf_gamma` - gamma distribution probability density function
+
+### Status
+
+Experimental
+
+### Description
+
+The probability density function (pdf) of the single real variable gamma distribution:
+
+$$ f(x)= \frac{scale^{shape}}{\Gamma (shape)}x^{shape-1}e^{-scale \times x} , for \;  x>0, shape, scale>0$$
+
+For a complex variable (x + y i) with independent real x and imaginary y parts, the joint probability density function is the product of the corresponding marginal pdf of real and imaginary pdf (for more details, see
+"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
+
+$$f(x+\mathit{i}y)=f(x)f(y)$$
+
+### Syntax
+
+`result = [[stdlib_stats_distribution_gamma(module):pdf_gamma(interface)]](x, shape, rate)`
+
+### Class
+
+Elemental function
+
+### Arguments
+
+`x`: has `intent(in)` and is a scalar of type `real` or `complex`.
+
+`shape` has `intent(in)` and is a scalar of type `real` or `complex`.
+
+`rate`: has `intent(in)` and is a scalar of type `real` or `complex`.
+
+All arguments must have the same type.
+
+### Return value
+
+The result is a scalar or an array, with a shape conformable to arguments, of type `real`.
+
+### Example
+
+```fortran
+program demo_gamma_pdf
+    use stdlib_random, only : random_seed
+    use stdlib_stats_distribution_gamma, only: rgamma => rvs_gamma,&
+                                             gamma_pdf => pdf_gamma
+
+    implicit none
+    real :: x(2,3,4),g(2,3,4),s(2,3,4)
+    complex :: shape, scale
+    integer :: put, get
+
+    put = 1234567
+    call random_seed(put, get)
+
+    print *, gamma_pdf(1.0, 1.0, 1.0)
+    !a probability density at 1.0 with shape=1.0, rate=1.0
+
+! 0.367879450
+
+    g(:,:,:) = 2.0
+    s(:,:,:) = 1.0
+    x = reshape(rgamma(2.0, 1.0, 24),[2,3,4]) ! gamma random variates array
+    print *, gamma_pdf(x,g,s)     ! a rank 3 gamma probability density array
+
+!  0.204550430  0.320178866  0.274986655  0.348611295  0.101865448 
+!  0.102199331  0.358981341  0.223676488  0.254329354  0.356714427 
+!  0.267390072  0.305148095  0.367848188  7.26194456E-02  1.49471285E-02 
+!  0.246272027  0.360770017  0.339665830  0.101558588  0.358678699 
+!  0.224196941  0.359253854  7.56355673E-02  0.251869917
+
+    shape = (1.0, 1.5)
+    scale = (1.0, 2.)
+    print *, gamma_pdf((1.5,1.0), shape, scale)
+    ! a complex expon probability density function at (1.5,1.0) with real part
+    !of shape=1.0, rate=1.0 and imaginary part of shape=1.5, rate=2.0
+
+! 9.63761061E-02
+
+end program demo_gamma_pdf
+```
+
+## `cdf_gamma` - gamma distribution cumulative distribution function
+
+### Status
+
+Experimental
+
+### Description
+
+Cumulative distribution function (cdf) of the single real variable gamma distribution:
+
+$$ F(x)= \frac{\gamma (shape, scale \times x)}{\Gamma (shape)}, for \;  x>0, shape, scale>0 $$
+
+For a complex variable (x + y i) with independent real x and imaginary y parts, the joint cumulative distribution function is the product of corresponding marginal cdf of real and imaginary cdf (for more details, see
+"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
+
+$$F(x+\mathit{i}y)=F(x)F(y)$$
+
+### Syntax
+
+`result = [[stdlib_stats_distribution_gamma(module):cdf_gamma(interface)]](x, shape, rate)`
+
+### Class
+
+Elemental function
+
+### Arguments
+
+`x`: has `intent(in)` and is a scalar of type `real` or `complex`.
+
+`shape`: has `intent(in)` and is a scalar of type `real` or `complex`.
+
+`rate`: has `intent(in)` and is a scalar of type `real` or `complex`.
+
+All arguments must have the same type.
+
+### Return value
+
+The result is a scalar of type `real` with the same kind of input arguments.
+
+### Example
+
+```fortran
+program demo_gamma_cdf
+    use stdlib_random, only : random_seed
+    use stdlib_stats_distribution_gamma, only: rgamma => rvs_gamma,&
+                                             gamma_cdf => cdf_gamma
+
+    implicit none
+    real :: x(2,3,4),g(2,3,4),s(2,3,4)
+    complex :: shape, scale
+    integer :: seed_put, seed_get
+
+    seed_put = 1234567
+    call random_seed(seed_put, seed_get)
+
+    print *, gamma_cdf(1.0, 0.5,1.0)
+    ! a standard gamma cumulative at 1.0 with a shape=0.5, rate=1.0
+
+! 0.842700839
+
+    print *, gamma_cdf(2.0, 1.5,2.0)
+    ! a cumulative at 2.0 with a shape=1.5, rate=2.0
+
+! 0.953988254
+
+    g(:,:,:) = 1.0
+    s(:,:,:) = 1.0
+    x = reshape(rgamma(1.0, 1.0, 24),[2,3,4])
+    !gamma random variates array with a shape=1.0, rate=1.0
+    print *, gamma_cdf(x,g,s)        ! a rank 3 standard gamma cumulative array
+
+!  0.710880339  0.472411335  0.578345954  0.383050948  0.870905757 
+!  0.870430350  0.170215249  0.677347481  0.620089889  0.161825046 
+!  4.17549349E-02  0.510665894  0.252201647  0.911497891  0.984424412 
+!  0.635621786  0.177783430  0.414842933  0.871342421  0.338317066 
+!  2.06879266E-02  0.335232288  0.907408893  0.624871135
+
+    shape = (.7, 2.1)
+    scale = (0.5,1.0)
+    print *, gamma_cdf((0.5,0.5),shape,scale)
+    !complex gamma cumulative distribution at (0.5,0.5) with real part of
+    !shape=0.7,rate=0.5 and imaginary part of shape=2.1,rate=1.0
+
+! 2.87349485E-02
+
+end program demo_gamma_cdf
+
+```
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8a6fe66cc..c991b997f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -43,6 +43,7 @@ set(fppFiles
     stdlib_stats_distribution_uniform.fypp
     stdlib_stats_distribution_normal.fypp
     stdlib_stats_distribution_exponential.fypp
+    stdlib_stats_distribution_gamma.fypp
     stdlib_stats_var.fypp
     stdlib_quadrature.fypp
     stdlib_quadrature_trapz.fypp
diff --git a/src/stdlib_stats_distribution_gamma.fypp b/src/stdlib_stats_distribution_gamma.fypp
new file mode 100644
index 000000000..664eddeb5
--- /dev/null
+++ b/src/stdlib_stats_distribution_gamma.fypp
@@ -0,0 +1,314 @@
+#:set WITH_QP = False
+#:set WITH_XDP = False
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+Module stdlib_stats_distribution_gamma
+    use stdlib_kinds, only : sp, dp
+    use stdlib_error, only : error_stop
+    use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform
+    use stdlib_stats_distribution_normal, only : rnor=>rvs_normal
+    use stdlib_specialfunctions_gamma, only : lincgam => lower_incomplete_gamma
+
+    implicit none
+    private
+
+    public :: rvs_gamma
+    public :: pdf_gamma
+    public :: cdf_gamma
+
+
+    interface rvs_gamma
+    !! Version experimental
+    !!
+    !! Gamma Distribution Random Variates
+    !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html#
+    !! rvs_gamma-gamma-distribution-random-variates))
+    !!
+        #:for k1, t1 in RC_KINDS_TYPES
+        module procedure gamma_dist_rvs_1_${t1[0]}$${k1}$     ! 1 argument
+        #:endfor
+
+        #:for k1, t1 in RC_KINDS_TYPES
+        module procedure gamma_dist_rvs_${t1[0]}$${k1}$       ! 2 arguments
+        #:endfor
+
+        #:for k1, t1 in RC_KINDS_TYPES
+        module procedure gamma_dist_rvs_array_${t1[0]}$${k1}$ ! 3 arguments
+        #:endfor
+    end interface rvs_gamma
+
+
+    interface pdf_gamma
+    !! Version experimental
+    !!
+    !! Gamma Distribution Probability Density Function
+    !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html#
+    !! pdf_gamma-gamma-distribution-probability-density-function))
+    !!
+        #:for k1, t1 in RC_KINDS_TYPES
+        module procedure gamma_dist_pdf_${t1[0]}$${k1}$
+        #:endfor
+    end interface pdf_gamma
+
+
+    interface cdf_gamma
+    !! Version experimental
+    !!
+    !! Gamma Distribution Cumulative Distribution Function
+    !! ([Specification](../page/specs/stdlib_stats_distribution_gamma.html#
+    !! cdf_gamma_gamma-distribution-cumulative-density-function))
+    !!
+        #:for k1, t1 in RC_KINDS_TYPES
+        module procedure gamma_dist_cdf_${t1[0]}$${k1}$
+        #:endfor
+    end interface cdf_gamma
+
+
+
+
+contains
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    impure elemental function gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) result(res)
+    ! Gamma distribution random variate. "A Simple Method for Generating Gamma
+    ! Variables", G. Marsaglia & W. W. Tsang, ACM Transactions on Mathematical
+    ! Software, 26(3), 2000, p. 363
+    !
+        ${t1}$, intent(in) :: shape
+        ${t1}$ :: res
+        ${t1}$ :: x, v, u, zz
+        ${t1}$, save :: alpha = 0._${k1}$, d, c
+        ${t1}$, parameter :: sq = 0.0331_${k1}$, tol = 1000 * epsilon(1.0_${k1}$)
+
+
+        if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs): Gamma"  &
+            //" distribution shape parameter must be greater than zero")
+
+        zz = shape
+
+        if(zz < 1._${k1}$) zz = 1._${k1}$ + zz
+        !shift shape parameter > 1
+        if(abs(zz - alpha) > tol) then
+        !initial run
+            alpha = zz
+            d = alpha - 1._${k1}$ / 3._${k1}$
+            c = 1._${k1}$ / (3._${k1}$ * sqrt(d))
+
+        endif
+
+        do
+            do
+                x = rnor(0.0_${k1}$, 1.0_${k1}$)
+                v = 1._${k1}$ + c * x
+                v = v * v * v
+
+                if(v > 0._${k1}$) exit
+
+            end do
+
+            x = x * x
+            u = uni(1.0_${k1}$)
+
+            if(u < (1._${k1}$ - sq * x * x)) exit
+
+            if(log(u) < 0.5_${k1}$ * x + d * (1._${k1}$ - v + log(v))) exit
+
+        end do
+
+        res = d * v
+
+        if(shape < 1._${k1}$) then
+        !restore shape parameter < 1
+            u = uni(1.0_${k1}$)
+            res = res * u ** (1._${k1}$ / shape)
+
+        endif
+    end function gamma_dist_rvs_1_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    impure elemental function gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) result(res)
+    ! Complex parameter gamma distributed. The real part and imaginary part are
+    ! independent of each other.
+    !
+        ${t1}$, intent(in) :: shape
+        ${t1}$ :: res
+
+        res = cmplx(gamma_dist_rvs_1_r${k1}$(shape%re),                        &
+                    gamma_dist_rvs_1_r${k1}$(shape%im), kind=${k1}$)
+    end function gamma_dist_rvs_1_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    impure elemental function gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate)      &
+        result(res)
+    !
+        ${t1}$, intent(in) :: shape, rate
+        ${t1}$ :: res
+
+        if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_rvs): Gamma"  &
+        //" distribution rate parameter must be greater than zero")
+
+        res = gamma_dist_rvs_1_${t1[0]}$${k1}$(shape) / rate
+    end function gamma_dist_rvs_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    impure elemental function gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate)      &
+        result(res)
+    ! Complex parameter gamma distributed. The real part and imaginary part are           &
+    ! independent of each other.
+    !
+        ${t1}$, intent(in) :: shape, rate
+        ${t1}$ :: res
+
+        res = cmplx(gamma_dist_rvs_r${k1}$(shape%re, rate%re),                 &
+                    gamma_dist_rvs_r${k1}$(shape%im, rate%im), kind=${k1}$)
+    end function gamma_dist_rvs_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    function gamma_dist_rvs_array_${t1[0]}$${k1}$(shape, rate, array_size)     &
+        result(res)
+    !
+        ${t1}$, intent(in) :: shape, rate
+        integer, intent(in) :: array_size
+        ${t1}$ :: res(array_size)
+        integer :: i
+
+        do i = 1, array_size
+
+            res(i) = gamma_dist_rvs_${t1[0]}$${k1}$(shape, rate)
+
+        end do
+    end function gamma_dist_rvs_array_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    function gamma_dist_rvs_array_${t1[0]}$${k1}$(shape, rate, array_size)     &
+        result(res)
+    ! Complex parameter gamma distributed. The real part and imaginary part are           &
+    ! independent of each other.
+    !
+        ${t1}$, intent(in) :: shape, rate
+        integer, intent(in) :: array_size
+        ${t1}$ :: res(array_size)
+        integer :: i
+
+        do i = 1, array_size
+
+            res(i) = cmplx(gamma_dist_rvs_r${k1}$(shape%re, rate%re),          &
+                           gamma_dist_rvs_r${k1}$(shape%im, rate%im),          &
+                           kind=${k1}$)
+
+        end do
+    end function gamma_dist_rvs_array_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    impure elemental function gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate)   &
+        result(res)
+    ! Gamma distribution probability density function
+    !
+        ${t1}$, intent(in) :: x, shape, rate
+        real(${k1}$) :: res
+
+        if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma"  &
+            //" distribution rate parameter must be greaeter than zero")
+
+        if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" &
+            //" distribution shape parameter must be greater than zero")
+
+        if(x <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma"     &
+            //" distribution variate x must be greater than zero")
+
+        if(x == 0.0_${k1}$) then
+
+            if(shape <= 1.0_${k1}$) then
+
+                res = huge(1.0) + 1.0
+
+            else
+
+                res = 0.0_${k1}$
+
+            endif
+
+        else
+
+            res = exp((shape - 1._${k1}$) * log(x) - x * rate + shape *        &
+                log(rate) - log_gamma(shape))
+
+        endif
+    end function gamma_dist_pdf_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    impure elemental function gamma_dist_pdf_${t1[0]}$${k1}$(x, shape, rate)    &
+        result(res)
+    ! Complex parameter gamma distributed. The real part and imaginary part are           &
+    ! independent of each other.
+    !
+        ${t1}$, intent(in) :: x, shape, rate
+        real(${k1}$) :: res
+
+        res = gamma_dist_pdf_r${k1}$(x%re, shape%re, rate%re)
+        res = res * gamma_dist_pdf_r${k1}$(x%im, shape%im, rate%im)
+    end function gamma_dist_pdf_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    impure elemental function gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate)   &
+        result(res)
+    ! Gamma distribution cumulative distribution function
+    !
+        ${t1}$, intent(in) :: x, shape, rate
+        real(${k1}$) :: res
+
+        if(rate <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma"  &
+            //" distribution rate parameter must be greaeter than zero")
+
+        if(shape <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma" &
+            //" distribution shape parameter must be greater than zero")
+
+        if(x <= 0.0_${k1}$) call error_stop("Error(gamma_dist_pdf): Gamma"     &
+            //" distribution variate x must be greater than zero")
+
+        res = lincgam(shape, rate * x) / gamma(shape)
+    end function gamma_dist_cdf_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    impure elemental function gamma_dist_cdf_${t1[0]}$${k1}$(x, shape, rate)    &
+        result(res)
+    ! Complex parameter gamma distributed. The real part and imaginary part are           &
+    ! independent of each other.
+    !
+        ${t1}$, intent(in) :: x, shape, rate
+        real(${k1}$) :: res
+
+        res = gamma_dist_cdf_r${k1}$(x%re, shape%re, rate%re)
+        res = res * gamma_dist_cdf_r${k1}$(x%im, shape%im, rate%im)
+    end function gamma_dist_cdf_${t1[0]}$${k1}$
+
+    #:endfor
+
+end module stdlib_stats_distribution_gamma
diff --git a/test/stats/CMakeLists.txt b/test/stats/CMakeLists.txt
index ff9d45063..94990ccf1 100644
--- a/test/stats/CMakeLists.txt
+++ b/test/stats/CMakeLists.txt
@@ -8,6 +8,7 @@ set(fppFiles
     test_distribution_uniform.fypp
     test_distribution_normal.fypp
     test_distribution_exponential.fypp
+    test_distribution_gamma.fypp
 )
 
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
@@ -24,6 +25,7 @@ ADDTEST(random)
 ADDTEST(distribution_uniform)
 ADDTEST(distribution_normal)
 ADDTEST(distribution_exponential)
+ADDTEST(distribution_gamma)
 
 if(DEFINED CMAKE_MAXIMUM_RANK)
     if(${CMAKE_MAXIMUM_RANK} GREATER 7)
diff --git a/test/stats/test_distribution_gamma.fypp b/test/stats/test_distribution_gamma.fypp
new file mode 100644
index 000000000..a49cfb7b1
--- /dev/null
+++ b/test/stats/test_distribution_gamma.fypp
@@ -0,0 +1,337 @@
+#:set WITH_QP = False
+#:set WITH_XDP = False
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+module test_stats_distribution_gamma
+    use testdrive, only : new_unittest, unittest_type, error_type, check
+    use stdlib_kinds, only : sp, dp
+    use stdlib_random, only : random_seed
+    use stdlib_stats_distribution_gamma, only : rgamma => rvs_gamma,           &
+                              gamma_pdf => pdf_gamma, gamma_cdf => cdf_gamma
+
+    implicit none
+    private
+
+    public :: collect_stats_distribution_gamma
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$)
+    #:endfor
+
+
+
+
+contains
+
+    subroutine collect_stats_distribution_gamma(testsuite)
+        type(unittest_type), allocatable, intent(out) :: testsuite(:)
+
+        testsuite = [                                                          &
+            new_unittest("gamma_random_generator", test_gamma_random_generator)&
+
+        #:for k1, t1 in RC_KINDS_TYPES
+            , new_unittest("gamma_rvs_${t1[0]}$${k1}$",                        &
+                           test_gamma_rvs_${t1[0]}$${k1}$)                     &
+        #:endfor
+
+        #:for k1, t1 in RC_KINDS_TYPES
+            , new_unittest("gamma_pdf_${t1[0]}$${k1}$",                        &
+                           test_gamma_pdf_${t1[0]}$${k1}$)                     &
+        #:endfor
+
+        #:for k1, t1 in RC_KINDS_TYPES
+            , new_unittest("gamma_cdf_${t1[0]}$${k1}$",                        &
+                           test_gamma_cdf_${t1[0]}$${k1}$)                     &
+        #:endfor
+            ]
+    end subroutine collect_stats_distribution_gamma
+
+
+
+    subroutine test_gamma_random_generator(error)
+        type(error_type), allocatable, intent(out) :: error
+        integer, parameter :: num = 10000000, array_size = 1000
+        integer :: i, j, freq(0 : array_size), put, get
+        real(dp) :: chisq, expct
+        character(80) :: ch1
+
+        print *, ""
+        print *, "Test gamma random generator with chi-squared"
+
+        put = 1234567
+        call random_seed(put, get)
+        freq = 0
+
+        do i = 1, num
+
+            j = 1000 * gamma_cdf(rgamma(2.0,1.5),2.0,1.5)
+            freq(j) = freq(j) + 1
+
+        end do
+
+        chisq = 0.0_dp
+        expct = num / array_size
+
+        do i = 0, array_size - 1
+
+            chisq = chisq + (freq(i) - expct) ** 2 / expct
+
+        end do
+
+        write(ch1, '(f10.1)') chisq
+
+        call check(error, chisq < 1143.9,                                       &
+                   "gamma randomness failed chi-squared test.",                 &
+                   "The critical values for chi-squared with 1000 dof is"//     &
+                   " 1143.92. Chi-squared for gamma random generator is :"//    &
+                   trim(ch1))
+
+    end subroutine test_gamma_random_generator
+
+
+
+    #:for k1, t1 in RC_KINDS_TYPES
+    subroutine test_gamma_rvs_${t1[0]}$${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        integer, parameter :: k = 5, n = 10
+        ${t1}$ :: res(n), gshape, scale
+        integer :: i
+        integer :: seed, get
+    #:if t1[0] == "r"
+      #! for real type
+        ${t1}$ :: ans(n) = [0.85758907497718884_${k1}$,                        &
+                            1.0206623865526090_${k1}$,                         &
+                            0.99753931024198650_${k1}$,                        &
+                            0.97653359790345839_${k1}$,                        &
+                            0.41853482638322043_${k1}$,                        &
+                            2.2012288073086310_${k1}$,                         &
+                            2.0639542613306592_${k1}$,                         &
+                            3.1794669730880192_${k1}$,                         &
+                            1.9329744662223280_${k1}$,                         &
+                            1.0257959670932111_${k1}$]
+    #:else
+      #! for complex type
+        ${t1}$ :: ans(n) =                                                     &
+                     [(1.0719863437214860_${k1}$, 0.46775532101393819_${k1}$), &
+                     (0.42382516926807201_${k1}$, 0.96340496644915230_${k1}$), &
+                      (2.7515360091357888_${k1}$, 0.14837198853150388_${k1}$), &
+                      (1.4536367104245524_${k1}$, 0.56852736336951559_${k1}$), &
+                (0.34559143458416125_${k1}$, 4.96217685362488267E-002_${k1}$), &
+                      (1.9657884897696516_${k1}$,  3.1124314799641013_${k1}$), &
+                 (3.4155160623540453_${k1}$, 5.04948933894018709E-002_${k1}$), &
+                     (0.94594398345216302_${k1}$, 0.45691588305890624_${k1}$), &
+                      (1.1493158751025965_${k1}$, 0.12944763723941669_${k1}$), &
+                      (2.9691469633592282_${k1}$, 1.1617408197125874_${k1}$)]
+    #:endif
+
+        print *, "Test gamma_distribution_rvs_${t1[0]}$${k1}$"
+        seed = 639741825
+        call random_seed(seed, get)
+
+    #:if t1[0] == "r"
+      #! for real type
+        gshape = 2.0_${k1}$; scale = 1.0_${k1}$
+    #:else
+      #! for complex type
+        gshape = (2.0_${k1}$, 0.7_${k1}$); scale = (0.8_${k1}$, 1.2_${k1}$)
+    #:endif
+
+        do i = 1, k
+
+            res(i) = rgamma(gshape, scale)
+
+        end do
+
+        res(k + 1 : n) = rgamma(gshape, scale, k)
+
+        do i = 1, n
+
+            call check(error, res(i), ans(i), "gamma_distribution_rvs_"//      &
+                      "${t1[0]}$${k1}$ failed", thr = tol_${k1}$)
+
+        end do
+    end subroutine test_gamma_rvs_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+
+
+    #:for k1, t1 in RC_KINDS_TYPES
+    subroutine test_gamma_pdf_${t1[0]}$${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ${t1}$ :: x1, x2(3,4), gshape, scale
+        integer :: i
+        integer :: seed, get
+        real(${k1}$) :: res(15)
+    #:if t1[0] == "r"
+      #! for real type
+        real(${k1}$), parameter :: ans(15) =                                   &
+                                   [3.4495412572168718E-002_${k1}$,            &
+                                    3.4495412572168718E-002_${k1}$,            &
+                                    3.4495412572168718E-002_${k1}$,            &
+                                    0.29116634347089576_${k1}$,                &
+                                    0.28338290850731412_${k1}$,                &
+                                    0.27922270935613586_${k1}$,                &
+                                    0.36440665523348270_${k1}$,                &
+                                    0.24379209619143699_${k1}$,                &
+                                    6.3815638087140858E-002_${k1}$,            &
+                                    0.25844600948718588_${k1}$,                &
+                                    0.17268118913523497_${k1}$,                &
+                                    0.31181223194308200_${k1}$,                &
+                                    0.24027095040543087_${k1}$,                &
+                                    0.36765502365831570_${k1}$,                &
+                                    9.9011714088769673E-002_${k1}$]
+    #:else
+      #! for complex type
+        real(${k1}$), parameter :: ans(15) =                                   &
+                                   [0.11554282574059289_${k1}$,                &
+                                    0.11554282574059289_${k1}$,                &
+                                    0.11554282574059289_${k1}$,                &
+                                    9.2682318951901529E-002_${k1}$,            &
+                                    0.40166849087286088_${k1}$,                &
+                                    0.37468980496232701_${k1}$,                &
+                                    0.14712363446345342_${k1}$,                &
+                                    0.22561628567985184_${k1}$,                &
+                                    0.12765403024301181_${k1}$,                &
+                                    3.9182498867847360E-002_${k1}$,            &
+                                    2.5873533461032859E-003_${k1}$,            &
+                                    0.10105832622792968_${k1}$,                &
+                                    0.24044091896609490_${k1}$,                &
+                                    4.9885356046115948E-003_${k1}$,            &
+                                    0.11085827028639164_${k1}$]
+    #:endif
+
+        print *, "Test gamma_distribution_pdf_${t1[0]}$${k1}$"
+        seed = 345987126
+        call random_seed(seed, get)
+    #:if t1[0] == "r"
+      #! for real type
+        gshape = 2.0_${k1}$; scale = 1.0_${k1}$
+    #:else
+      #! for complex type
+        gshape = (2.0_${k1}$, 0.7_${k1}$); scale = (0.8_${k1}$, 1.2_${k1}$)
+    #:endif
+
+        x1 = rgamma(gshape, scale)
+        x2 = reshape(rgamma(gshape, scale, 12), [3,4])
+        res(1:3) = gamma_pdf(x1, gshape, scale)
+        res(4:15) = reshape(gamma_pdf(x2, gshape, scale), [12])
+
+        do i = 1, 15
+
+                call check(error, res(i), ans(i), "gamma_distribution"// &
+                           "_pdf_${t1[0]}$${k1}$ failed", thr = tol_${k1}$)
+
+        end do
+    end subroutine test_gamma_pdf_${t1[0]}$${k1}$
+
+    #:endfor
+
+
+
+
+    #:for k1, t1 in RC_KINDS_TYPES
+    subroutine test_gamma_cdf_${t1[0]}$${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ${t1}$ :: x1, x2(3,4), gshape, scale
+        integer :: i
+        integer :: seed, get
+        real(${k1}$) :: res(15)
+    #:if t1[0] == "r"
+      #! for real type
+        real(${k1}$), parameter :: ans(15) =                                   &
+                                   [5.4876256610822634E-002_${k1}$,            &
+                                    5.4876256610822634E-002_${k1}$,            &
+                                    5.4876256610822634E-002_${k1}$,            &
+                                    0.31541195839514946_${k1}$,                &
+                                    0.38568161497244058_${k1}$,                &
+                                    0.23220859761573376_${k1}$,                &
+                                    0.39336687687155714_${k1}$,                &
+                                    0.80559422971604655_${k1}$,                &
+                                    0.88631934249921673_${k1}$,                &
+                                    0.37667963185005432_${k1}$,                &
+                                    0.14176369124149241_${k1}$,                &
+                                    0.45590880930769767_${k1}$,                &
+                                    0.27856937500418372_${k1}$,                &
+                                    0.18103305981618728_${k1}$,                &
+                                    0.72986385036366463_${k1}$]
+    #:else
+      #! for complex type
+        real(${k1}$), parameter :: ans(15) =                                   &
+                                   [3.2122120543510921E-002_${k1}$,            &
+                                    3.2122120543510921E-002_${k1}$,            &
+                                    3.2122120543510921E-002_${k1}$,            &
+                                    0.20931149671160035_${k1}$,                &
+                                    0.77957028981310350_${k1}$,                &
+                                    0.17082639598330887_${k1}$,                &
+                                    2.7594977080807291E-002_${k1}$,            &
+                                    2.3794072479821078E-002_${k1}$,            &
+                                    5.2298181677386930E-002_${k1}$,            &
+                                    0.22327051157236336_${k1}$,                &
+                                    0.27365315981967359_${k1}$,                &
+                                    3.4968870668437825E-002_${k1}$,            &
+                                    0.58026010546190465_${k1}$,                &
+                                    0.23090402450867176_${k1}$,                &
+                                    0.25072609802292339_${k1}$]
+    #:endif
+
+        print *, "Test gamma_distribution_cdf_${t1[0]}$${k1}$"
+        seed = 567985123
+        call random_seed(seed, get)
+
+    #:if t1[0] == "r"
+      #! for real type
+        gshape = 2.0_${k1}$; scale = 2.0_${k1}$
+    #:else
+      #! for complex type
+        gshape = (2.0_${k1}$, 0.7_${k1}$); scale = (0.8_${k1}$, 1.2_${k1}$)
+    #:endif
+
+        x1 = rgamma(gshape, scale)
+        x2 = reshape(rgamma(gshape, scale, 12), [3,4])
+        res(1:3) = gamma_cdf(x1, gshape, scale)
+        res(4:15) = reshape(gamma_cdf(x2, gshape, scale), [12])
+
+        do i = 1, 15
+
+                call check(error, res(i), ans(i), "gamma_distribution"// &
+                           "_cdf_${t1[0]}$${k1}$ failed", thr = tol_${k1}$)
+
+        end do
+    end subroutine test_gamma_cdf_${t1[0]}$${k1}$
+
+    #:endfor
+
+end module test_stats_distribution_gamma
+
+
+
+program tester
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_stats_distribution_gamma, only : collect_stats_distribution_gamma
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [new_testsuite("Stats_distribution_gamma",                    &
+                  collect_stats_distribution_gamma)]
+
+    do is = 1, size(testsuites)
+
+        write(error_unit, fmt) "Testing:", testsuites(is) % name
+        call run_testsuite(testsuites(is) % collect, error_unit, stat)
+
+    end do
+
+    if(stat > 0) then
+
+        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        error stop
+
+    end if
+end program tester