diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp
index de4f4b275..35788ef1e 100644
--- a/src/stdlib_stats_distribution_normal.fypp
+++ b/src/stdlib_stats_distribution_normal.fypp
@@ -1,10 +1,11 @@
 #:include "common.fypp"
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 module stdlib_stats_distribution_normal
-    use stdlib_kinds, only : sp, dp, xdp, qp, int32
-    use stdlib_error, only : error_stop
-    use stdlib_random, only : dist_rand
-    use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform
+    use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+    use stdlib_kinds, only: sp, dp, xdp, qp, int32
+    use stdlib_error, only: error_stop
+    use stdlib_random, only: dist_rand
+    use stdlib_stats_distribution_uniform, only: uni => rvs_uniform
 
     implicit none
     private
@@ -18,8 +19,6 @@ module stdlib_stats_distribution_normal
     public :: pdf_normal
     public :: cdf_normal
 
-
-
     interface rvs_normal
     !! version: experimental
     !!
@@ -30,16 +29,14 @@ module stdlib_stats_distribution_normal
         module procedure rvs_norm_0_rsp                  !0 dummy variable
 
         #:for k1, t1 in RC_KINDS_TYPES
-        module procedure rvs_norm_${t1[0]}$${k1}$        !2 dummy variables
+            module procedure rvs_norm_${t1[0]}$${k1}$        !2 dummy variables
         #:endfor
 
         #:for k1, t1 in RC_KINDS_TYPES
-        module procedure rvs_norm_array_${t1[0]}$${k1}$  !3 dummy variables
+            module procedure rvs_norm_array_${t1[0]}$${k1}$  !3 dummy variables
         #:endfor
     end interface rvs_normal
 
-
-
     interface pdf_normal
     !! version: experimental
     !!
@@ -48,12 +45,10 @@ module stdlib_stats_distribution_normal
     !! pdf_normal-normal-distribution-probability-density-function))
     !!
         #:for k1, t1 in RC_KINDS_TYPES
-        module procedure pdf_norm_${t1[0]}$${k1}$
+            module procedure pdf_norm_${t1[0]}$${k1}$
         #:endfor
     end interface pdf_normal
 
-
-
     interface cdf_normal
     !! version: experimental
     !!
@@ -62,272 +57,264 @@ module stdlib_stats_distribution_normal
     !! cdf_normal-normal-distribution-cumulative-distribution-function))
     !!
         #:for k1, t1 in RC_KINDS_TYPES
-        module procedure cdf_norm_${t1[0]}$${k1}$
+            module procedure cdf_norm_${t1[0]}$${k1}$
         #:endfor
     end interface cdf_normal
 
-
-
-
-
 contains
 
-    subroutine zigset
-    ! Marsaglia & Tsang generator for random normals & random exponentials.
-    ! Translated from C by Alan Miller (amiller@bigpond.net.au), released as public
-    ! domain (https://jblevins.org/mirror/amiller/)
-    !
-    ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
-    ! random variables', J. Statist. Software, v5(8).
-    !
-    ! This is an electronic journal which can be downloaded from:
-    ! http://www.jstatsoft.org/v05/i08
-    !
-    ! Latest version - 1 January 2001
-    !
+    impure subroutine zigset
+        ! Marsaglia & Tsang generator for random normals & random exponentials.
+        ! Translated from C by Alan Miller (amiller@bigpond.net.au), released as public
+        ! domain (https://jblevins.org/mirror/amiller/)
+        !
+        ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
+        ! random variables', J. Statist. Software, v5(8).
+        !
+        ! This is an electronic journal which can be downloaded from:
+        ! http://www.jstatsoft.org/v05/i08
+        !
+        ! Latest version - 1 January 2001
+        !
         real(dp), parameter :: M1 = 2147483648.0_dp, vn = 0.00991256303526217_dp
         real(dp)            :: dn, tn, q
         integer :: i
 
         dn = 3.442619855899_dp
         tn = dn
-    !tables for random normals
-        q = vn * exp(HALF * dn * dn)
-        kn(0) = int((dn / q) * M1, kind = int32)
+        !tables for random normals
+        q = vn*exp(HALF*dn*dn)
+        kn(0) = int((dn/q)*M1, kind=int32)
         kn(1) = 0
-        wn(0) = q / M1
-        wn(127) = dn / M1
+        wn(0) = q/M1
+        wn(127) = dn/M1
         fn(0) = ONE
-        fn(127) = exp( -HALF * dn * dn )
-        do  i = 126, 1, -1
-            dn = sqrt( -TWO * log( vn / dn + exp( -HALF * dn * dn ) ) )
-            kn(i+1) = int((dn / tn) * M1, kind = int32)
+        fn(127) = exp(-HALF*dn*dn)
+        do i = 126, 1, -1
+            dn = sqrt(-TWO*log(vn/dn + exp(-HALF*dn*dn)))
+            kn(i + 1) = int((dn/tn)*M1, kind=int32)
             tn = dn
-            fn(i) = exp(-HALF * dn * dn)
-            wn(i) = dn / M1
+            fn(i) = exp(-HALF*dn*dn)
+            wn(i) = dn/M1
         end do
         zig_norm_initialized = .true.
     end subroutine zigset
 
-
-
-    #:for k1, t1 in REAL_KINDS_TYPES
-    function rvs_norm_0_${t1[0]}$${k1}$( ) result(res)
-    !
-    ! Standard normal random vairate (0,1)
-    !
-        ${t1}$ :: res
-        ${t1}$, parameter  ::  r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r
-        ${t1}$ ::  x, y
-        integer :: hz, iz
-
-        if(.not. zig_norm_initialized) call zigset
-        iz = 0
-        hz = dist_rand(1_int32)          !32bit random integer
-        iz = iand( hz, 127 )             !random integer in [0, 127]
-        if( abs( hz ) < kn(iz) ) then
-            res = hz * wn(iz)
-        else
-            L1: do
-                L2: if( iz == 0 ) then
-                    do
-                        x = -log( uni(1.0_${k1}$) ) * rr
-                        y = -log( uni(1.0_${k1}$) )
-                        if( y + y >= x * x ) exit
-                    end do
-                    res = r + x
-                    if( hz <= 0 ) res = -res
-                    exit L1
-                end if L2
-                x = hz * wn(iz)
-                if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) <            &
-                    exp(-HALF * x * x) ) then
-                    res = x
-                    exit L1
-                end if
-                hz = dist_rand(1_int32)
-                iz = iand( hz, 127 )
-                if( abs( hz ) < kn(iz) ) then
-                    res = hz * wn(iz)
-                    exit L1
-                end if
-            end do L1
-        end if
-    end function rvs_norm_0_${t1[0]}$${k1}$
-
-    #:endfor
-
-
-
-    #:for k1, t1 in REAL_KINDS_TYPES
-    impure elemental &
-    function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
-    !
-    ! Normal random variate (loc, scale)
-    !
-        ${t1}$, intent(in) :: loc, scale
-        ${t1}$ :: res
-
-        if(scale == 0._${k1}$) call error_stop("Error(rvs_norm): Normal"       &
-            //" distribution scale parameter must be non-zero")
-        res = rvs_norm_0_${t1[0]}$${k1}$(  )
-        res = res * scale + loc
-    end function rvs_norm_${t1[0]}$${k1}$
-
-    #:endfor
-
-
-
-    #:for k1, t1 in CMPLX_KINDS_TYPES
-    impure elemental &
-    function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
-    !
-    ! Normally distributed complex. The real part and imaginary part are       &
-    ! independent of each other.
-    !
-        ${t1}$, intent(in) :: loc, scale
-        ${t1}$ :: res
-        real(${k1}$) :: tr, ti
-
-        tr = rvs_norm_r${k1}$(loc % re, scale % re)
-        ti = rvs_norm_r${k1}$(loc % im, scale % im)
-        res = cmplx(tr, ti, kind=${k1}$)
-    end function rvs_norm_${t1[0]}$${k1}$
-
-    #:endfor
-
-
-
     #:for k1, t1 in REAL_KINDS_TYPES
-    function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
-        ${t1}$, intent(in) :: loc, scale
-        integer, intent(in) :: array_size
-        ${t1}$ :: res(array_size)
-        ${t1}$, parameter  ::  r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r
-        ${t1}$ ::  x, y, re
-        integer :: hz, iz, i
-
-        if(scale == 0._${k1}$) call error_stop("Error(rvs_norm_array): Normal" &
-            //"distribution scale parameter must be non-zero")
-        if(.not. zig_norm_initialized) call zigset
-        do i = 1, array_size
+        impure function rvs_norm_0_${t1[0]}$${k1}$ () result(res)
+            !
+            ! Standard normal random variate (0,1)
+            !
+            ${t1}$ :: res
+            ${t1}$, parameter  ::  r = 3.442619855899_${k1}$, rr = 1.0_${k1}$/r
+            ${t1}$ ::  x, y
+            integer :: hz, iz
+
+            if (.not. zig_norm_initialized) call zigset
             iz = 0
-            hz = dist_rand(1_int32)
-            iz = iand( hz, 127 )
-            if( abs( hz ) < kn(iz) ) then
-                re = hz * wn(iz)
+            hz = dist_rand(1_int32)          !32bit random integer
+            iz = iand(hz, 127)             !random integer in [0, 127]
+            if (abs(hz) < kn(iz)) then
+                res = hz*wn(iz)
             else
                 L1: do
-                    L2: if( iz == 0 ) then
+                    L2: if (iz == 0) then
                         do
-                            x = -log( uni(1.0_${k1}$) ) * rr
-                            y = -log( uni(1.0_${k1}$) )
-                            if( y + y >= x * x ) exit
+                            x = -log(uni(1.0_${k1}$))*rr
+                            y = -log(uni(1.0_${k1}$))
+                            if (y + y >= x*x) exit
                         end do
-                        re = r + x
-                        if( hz <= 0 ) re = -re
+                        res = r + x
+                        if (hz <= 0) res = -res
                         exit L1
                     end if L2
-                    x = hz * wn(iz)
-                    if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) <        &
-                        exp(-HALF * x * x) ) then
-                        re = x
+                    x = hz*wn(iz)
+                    if (fn(iz) + uni(1.0_${k1}$)*(fn(iz - 1) - fn(iz)) < &
+                        exp(-HALF*x*x)) then
+                        res = x
                         exit L1
                     end if
-
                     hz = dist_rand(1_int32)
-                    iz = iand( hz, 127 )
-                    if( abs( hz ) < kn(iz) ) then
-                        re = hz * wn(iz)
+                    iz = iand(hz, 127)
+                    if (abs(hz) < kn(iz)) then
+                        res = hz*wn(iz)
                         exit L1
                     end if
                 end do L1
             end if
-            res(i) = re * scale + loc
-        end do
-    end function rvs_norm_array_${t1[0]}$${k1}$
+        end function rvs_norm_0_${t1[0]}$${k1}$
 
     #:endfor
 
+    #:for k1, t1 in REAL_KINDS_TYPES
+        impure elemental &
+            function rvs_norm_${t1[0]}$${k1}$ (loc, scale) result(res)
+            !
+            ! Normal random variate (loc, scale)
+            !
+            ${t1}$, intent(in) :: loc, scale
+            ${t1}$ :: res
+
+            if (scale <= 0._${k1}$) then
+                res = ieee_value(1._${k1}$, ieee_quiet_nan)
+            else
+                res = rvs_norm_0_${t1[0]}$${k1}$ ()
+                res = res*scale + loc
+            end if
 
-
-    #:for k1, t1 in CMPLX_KINDS_TYPES
-    function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
-        ${t1}$, intent(in) :: loc, scale
-        integer, intent(in) :: array_size
-        integer :: i
-        ${t1}$ :: res(array_size)
-        real(${k1}$) :: tr, ti
-
-        do i = 1, array_size
-            tr = rvs_norm_r${k1}$(loc % re, scale % re)
-            ti = rvs_norm_r${k1}$(loc % im, scale % im)
-            res(i) = cmplx(tr, ti, kind=${k1}$)
-        end do
-    end function rvs_norm_array_${t1[0]}$${k1}$
+        end function rvs_norm_${t1[0]}$${k1}$
 
     #:endfor
 
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+        impure elemental function rvs_norm_${t1[0]}$${k1}$ (loc, scale) result(res)
+            !
+            ! Normally distributed complex. The real part and imaginary part are       &
+            ! independent of each other.
+            !
+            ${t1}$, intent(in) :: loc, scale
+            ${t1}$ :: res
+            real(${k1}$) :: tr, ti
 
+            tr = rvs_norm_r${k1}$ (loc%re, scale%re)
+            ti = rvs_norm_r${k1}$ (loc%im, scale%im)
+            res = cmplx(tr, ti, kind=${k1}$)
 
-    #:for k1, t1 in REAL_KINDS_TYPES
-    impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
-    !
-    ! Normal distribution probability density function
-    !
-        ${t1}$, intent(in) :: x, loc, scale
-        ${t1}$ :: res
-        ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$))
-
-        if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal"       &
-            //"distribution scale parameter must be non-zero")
-        res = exp(- 0.5_${k1}$ * ((x - loc) / scale) * (x - loc) / scale) /    &
-              (sqrt_2_Pi * scale)
-    end function pdf_norm_${t1[0]}$${k1}$
+        end function rvs_norm_${t1[0]}$${k1}$
 
     #:endfor
 
+    #:for k1, t1 in REAL_KINDS_TYPES
+        impure function rvs_norm_array_${t1[0]}$${k1}$ (loc, scale, array_size) result(res)
+            ${t1}$, intent(in) :: loc, scale
+            integer, intent(in) :: array_size
+            ${t1}$ :: res(array_size)
+            ${t1}$, parameter  ::  r = 3.442619855899_${k1}$, rr = 1.0_${k1}$/r
+            ${t1}$ ::  x, y, re
+            integer :: hz, iz, i
+
+            if (.not. zig_norm_initialized) call zigset
+
+            if (scale <= 0._${k1}$) then
+                res = ieee_value(1._${k1}$, ieee_quiet_nan)
+                return
+            end if
+
+            do i = 1, array_size
+                iz = 0
+                hz = dist_rand(1_int32)
+                iz = iand(hz, 127)
+                if (abs(hz) < kn(iz)) then
+                    re = hz*wn(iz)
+                else
+                    L1: do
+                        L2: if (iz == 0) then
+                            do
+                                x = -log(uni(1.0_${k1}$))*rr
+                                y = -log(uni(1.0_${k1}$))
+                                if (y + y >= x*x) exit
+                            end do
+                            re = r + x
+                            if (hz <= 0) re = -re
+                            exit L1
+                        end if L2
+                        x = hz*wn(iz)
+                        if (fn(iz) + uni(1.0_${k1}$)*(fn(iz - 1) - fn(iz)) < &
+                            exp(-HALF*x*x)) then
+                            re = x
+                            exit L1
+                        end if
+
+                        hz = dist_rand(1_int32)
+                        iz = iand(hz, 127)
+                        if (abs(hz) < kn(iz)) then
+                            re = hz*wn(iz)
+                            exit L1
+                        end if
+                    end do L1
+                end if
+                res(i) = re*scale + loc
+            end do
+        end function rvs_norm_array_${t1[0]}$${k1}$
 
+    #:endfor
 
     #:for k1, t1 in CMPLX_KINDS_TYPES
-    impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
-        ${t1}$, intent(in) :: x, loc, scale
-        real(${k1}$) :: res
+        impure function rvs_norm_array_${t1[0]}$${k1}$ (loc, scale, array_size) result(res)
+            ${t1}$, intent(in) :: loc, scale
+            integer, intent(in) :: array_size
+            integer :: i
+            ${t1}$ :: res(array_size)
+            real(${k1}$) :: tr, ti
 
-        res = pdf_norm_r${k1}$(x % re, loc % re, scale % re)
-        res = res * pdf_norm_r${k1}$(x % im, loc % im, scale % im)
-    end function pdf_norm_${t1[0]}$${k1}$
+            do i = 1, array_size
+                tr = rvs_norm_r${k1}$ (loc%re, scale%re)
+                ti = rvs_norm_r${k1}$ (loc%im, scale%im)
+                res(i) = cmplx(tr, ti, kind=${k1}$)
+            end do
+
+        end function rvs_norm_array_${t1[0]}$${k1}$
 
     #:endfor
 
+    #:for k1, t1 in REAL_KINDS_TYPES
+        elemental function pdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
+            !
+            ! Normal distribution probability density function
+            !
+            ${t1}$, intent(in) :: x, loc, scale
+            ${t1}$ :: res
+            ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$*acos(-1.0_${k1}$))
+
+            if (scale <= 0._${k1}$) then
+                res = ieee_value(1._${k1}$, ieee_quiet_nan)
+            else
+                res = exp(-0.5_${k1}$*((x - loc)/scale)*(x - loc)/scale)/ &
+                      (sqrt_2_Pi*scale)
+            end if
 
+        end function pdf_norm_${t1[0]}$${k1}$
 
-    #:for k1, t1 in REAL_KINDS_TYPES
-    impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
-    !
-    ! Normal distribution cumulative distribution function
-    !
-        ${t1}$, intent(in) :: x, loc, scale
-        ${t1}$ :: res
-        ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$)
-
-        if(scale == 0._${k1}$) call error_stop("Error(cdf_norm): Normal"       &
-            //"distribution scale parameter must be non-zero")
-        res = erfc(- (x - loc) / (scale * sqrt_2)) / 2.0_${k1}$
-    end function cdf_norm_${t1[0]}$${k1}$
+    #:endfor
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+        elemental function pdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
+            ${t1}$, intent(in) :: x, loc, scale
+            real(${k1}$) :: res
+
+            res = pdf_norm_r${k1}$ (x%re, loc%re, scale%re)
+            res = res*pdf_norm_r${k1}$ (x%im, loc%im, scale%im)
+        end function pdf_norm_${t1[0]}$${k1}$
 
     #:endfor
 
+    #:for k1, t1 in REAL_KINDS_TYPES
+        elemental function cdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
+            !
+            ! Normal distribution cumulative distribution function
+            !
+            ${t1}$, intent(in) :: x, loc, scale
+            ${t1}$ :: res
+            ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$)
+
+            if (scale <= 0._${k1}$) then
+                res = ieee_value(1._${k1}$, ieee_quiet_nan)
+            else
+                res = erfc(-(x - loc)/(scale*sqrt_2))/2.0_${k1}$
+            end if
+
+        end function cdf_norm_${t1[0]}$${k1}$
 
+    #:endfor
 
     #:for k1, t1 in CMPLX_KINDS_TYPES
-    impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
-        ${t1}$, intent(in) :: x, loc, scale
-        real(${k1}$) :: res
+        elemental function cdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
+            ${t1}$, intent(in) :: x, loc, scale
+            real(${k1}$) :: res
 
-        res = cdf_norm_r${k1}$(x % re, loc % re, scale % re)
-        res = res * cdf_norm_r${k1}$(x % im, loc % im, scale % im)
-    end function cdf_norm_${t1[0]}$${k1}$
+            res = cdf_norm_r${k1}$ (x%re, loc%re, scale%re)
+            res = res*cdf_norm_r${k1}$ (x%im, loc%im, scale%im)
+        end function cdf_norm_${t1[0]}$${k1}$
 
     #:endfor