diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 7217171e3..6cd76cf24 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -23,6 +23,8 @@ With two arguments, the function returns a normal distributed random variate \(N With three arguments, the function returns a rank-1 array of normal distributed random variates. +With two arguments `array_size` and `mold`, the function returns a rank-1 array of standard normal distributed random variates \(N(0,1)\), where the `mold` argument is used only to determine the output type and kind. + @note The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1] @@ -30,6 +32,8 @@ The algorithm used for generating exponential random variates is fundamentally l `result = ` [[stdlib_stats_distribution_normal(module):rvs_normal(interface)]] `([loc, scale] [[, array_size]])` +`result = ` [[stdlib_stats_distribution_normal(module):rvs_normal(interface)]] `(array_size, mold)` + ### Class Elemental function (passing both `loc` and `scale`). @@ -40,13 +44,15 @@ Elemental function (passing both `loc` and `scale`). `scale`: optional argument has `intent(in)` and is a positive scalar of type `real` or `complex`. -`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. When used with `loc` and `scale`, specifies the size of the output array. When used with `mold`, must be provided as the first argument. + +`mold`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. Used only to determine the type and kind of the output; its value is not referenced. When provided, generates standard normal variates \(N(0,1)\). `loc` and `scale` arguments must be of the same type. ### Return value -The result is a scalar or rank-1 array, with a size of `array_size`, and the same type as `scale` and `loc`. If `scale` is non-positive, the result is `NaN`. +The result is a scalar or rank-1 array, with a size of `array_size`, and the same type as `scale` and `loc` (or same type and kind as `mold` when using the `array_size, mold` form). If `scale` is non-positive, the result is `NaN`. ### Example diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 3ee50ea6f..69201b648 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -33,6 +33,7 @@ module stdlib_stats_distribution_normal #:for k1, t1 in RC_KINDS_TYPES module procedure rvs_norm_array_${t1[0]}$${k1}$ !3 dummy variables + module procedure rvs_norm_array_default_${t1[0]}$${k1}$ !2 dummy variables (mold, array_size) #:endfor end interface rvs_normal @@ -238,6 +239,22 @@ contains #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + impure function rvs_norm_array_default_${t1[0]}$${k1}$ (array_size, mold) 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 + ! + integer, intent(in) :: array_size + ${t1}$, intent(in) :: mold + ${t1}$ :: res(array_size) + + res = rvs_norm_array_${t1[0]}$${k1}$ (0.0_${k1}$, 1.0_${k1}$, array_size) + + 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 @@ -256,6 +273,25 @@ contains #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + impure function rvs_norm_array_default_${t1[0]}$${k1}$ (array_size, mold) 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 + ! + integer, intent(in) :: array_size + ${t1}$, intent(in) :: mold + ${t1}$ :: res(array_size) + + ! Call the full procedure with default loc=(0,0), scale=(1,1) + res = rvs_norm_array_${t1[0]}$${k1}$ (cmplx(0.0_${k1}$, 0.0_${k1}$, kind=${k1}$), & + cmplx(1.0_${k1}$, 1.0_${k1}$, kind=${k1}$), & + array_size) + + end function rvs_norm_array_default_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES elemental function pdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res) ! diff --git a/test/stats/test_distribution_normal.fypp b/test/stats/test_distribution_normal.fypp index 82e6faca7..8e959fab1 100644 --- a/test/stats/test_distribution_normal.fypp +++ b/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(10, mold) + + 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