Add default (mold, array_size) overloads for rvs_normal

src/stdlib_stats_distribution_normal.fypp

@@ -34,6 +34,10 @@ module stdlib_stats_distribution_normal
         #:for k1, t1 in RC_KINDS_TYPES
             module procedure rvs_norm_array_${t1[0]}$${k1}$  !3 dummy variables
         #:endfor
+
+        #:for k1, t1 in RC_KINDS_TYPES
+            module procedure rvs_norm_array_default_${t1[0]}$${k1}$  !2 dummy variables (mold, array_size)
+        #:endfor
     end interface rvs_normal
 
     interface pdf_normal
@@ -238,6 +242,82 @@ contains
 
     #:endfor
 
+    #:for k1, t1 in REAL_KINDS_TYPES
+        impure function rvs_norm_array_default_${t1[0]}$${k1}$ (mold, array_size) result(res)
+            !
+            ! Standard normal array random variate with default loc=0, scale=1
+            ! The mold argument is used only to determine the type and is not referenced
+            !
+            ${t1}$, intent(in) :: mold
+            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
+
+            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  ! Default: loc=0, scale=1, so re*1 + 0 = re
+            end do
+        end function rvs_norm_array_default_${t1[0]}$${k1}$
+
+    #:endfor
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+        impure function rvs_norm_array_default_${t1[0]}$${k1}$ (mold, array_size) result(res)
+            !
+            ! Standard normal complex array random variate with default loc=0, scale=1
+            ! The mold argument is used only to determine the type and is not referenced
+            !
+            ${t1}$, intent(in) :: mold
+            integer, intent(in) :: array_size
+            integer :: i
+            ${t1}$ :: res(array_size)
+            real(${k1}$) :: tr, ti
+
+            do i = 1, array_size
+                tr = rvs_norm_0_r${k1}$ ()
+                ti = rvs_norm_0_r${k1}$ ()
+                res(i) = cmplx(tr, ti, kind=${k1}$)
+            end do
+
+        end function rvs_norm_array_default_${t1[0]}$${k1}$
+
+    #:endfor
+
     #:for k1, t1 in CMPLX_KINDS_TYPES
         impure function rvs_norm_array_${t1[0]}$${k1}$ (loc, scale, array_size) result(res)
             ${t1}$, intent(in) :: loc, scale

test/stats/test_distribution_normal.fypp

@@ -26,6 +26,10 @@ program test_distribution_normal
     call test_nor_rvs_${t1[0]}$${k1}$
     #:endfor
 
+    #:for k1, t1 in RC_KINDS_TYPES
+    call test_nor_rvs_default_${t1[0]}$${k1}$
+    #:endfor
+
 
 
     #:for k1, t1 in RC_KINDS_TYPES
@@ -138,6 +142,42 @@ contains
     #:endfor
 
 
+    #:for k1, t1 in RC_KINDS_TYPES
+    subroutine test_nor_rvs_default_${t1[0]}$${k1}$
+        ${t1}$ :: a1(10), a2(10), mold
+        integer :: i
+        integer :: seed, get
+
+        print *, "Test normal_distribution_rvs_default_${t1[0]}$${k1}$"
+        seed = 25836914
+        call random_seed(seed, get)
+
+        ! explicit form with loc=0, scale=1
+        #:if t1[0] == "r"
+            a1 = nor_rvs(0.0_${k1}$, 1.0_${k1}$, 10)
+        #:else
+            a1 = nor_rvs((0.0_${k1}$, 0.0_${k1}$), (1.0_${k1}$, 1.0_${k1}$), 10)
+        #:endif
+
+        ! reset seed to reproduce same random sequence
+        seed = 25836914
+        call random_seed(seed, get)
+
+        ! default mold form: mold used only to disambiguate kind
+        #:if t1[0] == "r"
+            mold = 0.0_${k1}$
+        #:else
+            mold = (0.0_${k1}$, 0.0_${k1}$)
+        #:endif
+
+        a2 = nor_rvs(mold, 10)
+
+        call check(all(a1 == a2), msg="normal_distribution_rvs_default_${t1[0]}$${k1}$ failed", warn=warn)
+    end subroutine test_nor_rvs_default_${t1[0]}$${k1}$
+
+    #:endfor
+
+
 
 
     #:for k1, t1 in RC_KINDS_TYPES