From b2ff4bec6f00bf430134eeeb7af9a8b3d062c7c3 Mon Sep 17 00:00:00 2001 From: gururaj1512 Date: Wed, 12 Nov 2025 23:43:13 +0530 Subject: [PATCH 1/5] feat: add method rvs-normal-array-default --- src/stdlib_stats_distribution_normal.fypp | 80 +++++++++++++++++++++++ test/stats/test_distribution_normal.fypp | 40 ++++++++++++ 2 files changed, 120 insertions(+) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 3ee50ea6f..0bdff92cf 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/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 diff --git a/test/stats/test_distribution_normal.fypp b/test/stats/test_distribution_normal.fypp index 82e6faca7..89c0d711f 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(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 From c11973550304eb6abcc5f39996f2bd5b58a3f897 Mon Sep 17 00:00:00 2001 From: gururaj1512 Date: Sat, 15 Nov 2025 15:45:52 +0530 Subject: [PATCH 2/5] refactor as per suggested changes --- src/stdlib_stats_distribution_normal.fypp | 82 ++++++----------------- test/stats/test_distribution_normal.fypp | 2 +- 2 files changed, 20 insertions(+), 64 deletions(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 0bdff92cf..69201b648 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -33,9 +33,6 @@ 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 @@ -243,96 +240,55 @@ contains #:endfor #:for k1, t1 in REAL_KINDS_TYPES - impure function rvs_norm_array_default_${t1[0]}$${k1}$ (mold, array_size) result(res) + 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 ! - ${t1}$, intent(in) :: mold integer, intent(in) :: array_size + ${t1}$, intent(in) :: mold ${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 + res = rvs_norm_array_${t1[0]}$${k1}$ (0.0_${k1}$, 1.0_${k1}$, array_size) - 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 + 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 do i = 1, array_size - tr = rvs_norm_0_r${k1}$ () - ti = rvs_norm_0_r${k1}$ () + 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_default_${t1[0]}$${k1}$ + end function rvs_norm_array_${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 + 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 - integer :: i + ${t1}$, intent(in) :: mold ${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 + ! 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_${t1[0]}$${k1}$ + end function rvs_norm_array_default_${t1[0]}$${k1}$ #:endfor diff --git a/test/stats/test_distribution_normal.fypp b/test/stats/test_distribution_normal.fypp index 89c0d711f..8e959fab1 100644 --- a/test/stats/test_distribution_normal.fypp +++ b/test/stats/test_distribution_normal.fypp @@ -170,7 +170,7 @@ contains mold = (0.0_${k1}$, 0.0_${k1}$) #:endif - a2 = nor_rvs(mold, 10) + 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}$ From 6e2fd8ea01a76b7b0a4432ebe8db3724ffd32b22 Mon Sep 17 00:00:00 2001 From: gururaj1512 Date: Sun, 16 Nov 2025 12:23:26 +0530 Subject: [PATCH 3/5] chore: update specs --- doc/specs/stdlib_stats_distribution_normal.md | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) 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 From bb56acb7d6e3f90e166e0a23e2823b134931cee1 Mon Sep 17 00:00:00 2001 From: gururaj1512 Date: Fri, 21 Nov 2025 23:24:46 +0530 Subject: [PATCH 4/5] refactor: as per suggested changes --- doc/specs/stdlib_stats_distribution_normal.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 6cd76cf24..e14fa989d 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -52,7 +52,7 @@ Elemental function (passing both `loc` and `scale`). ### 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` (or same type and kind as `mold` when using the `array_size, mold` form). 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 From 04fe63aade930fc966102792df9ef718c6bb23af Mon Sep 17 00:00:00 2001 From: gururaj1512 Date: Thu, 4 Dec 2025 19:18:32 +0530 Subject: [PATCH 5/5] refactor: update implementation as suggested --- doc/specs/stdlib_stats_distribution_normal.md | 10 +++++----- src/stdlib_stats_distribution_normal.fypp | 4 ++-- test/stats/test_distribution_normal.fypp | 15 ++++++++++----- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index e14fa989d..b1d364891 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -23,7 +23,7 @@ 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. +With one or two arguments where the first is `array_size`, the function returns a rank-1 array of standard normal distributed random variates \(N(0,1)\). The `mold` argument determines the output type and kind; it is optional only for `real(dp)` (and defaults to `real(dp)` when omitted), but required for all other types. @note The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1] @@ -32,7 +32,7 @@ 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)` +`result = ` [[stdlib_stats_distribution_normal(module):rvs_normal(interface)]] `(array_size [, mold])` ### Class @@ -44,15 +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`. When used with `loc` and `scale`, specifies the size of the output array. When used with `mold`, must be provided as the first argument. +`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 alone or 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)\). +`mold`: optional argument (only for `real(dp)`; required for other types) 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 omitted (only allowed for `real(dp)`), defaults to `real(dp)`. When provided, generates standard normal variates \(N(0,1)\) of the specified type and kind. `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` (or same type and kind as `mold` when using the `(array_size, mold)` form). 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; defaults to `real(dp)` when `mold` is omitted). 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 69201b648..9059b32de 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -33,7 +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) + module procedure rvs_norm_array_default_${t1[0]}$${k1}$ !array_size, mold (mold optional for real(dp) only) #:endfor end interface rvs_normal @@ -246,7 +246,7 @@ contains ! The mold argument is used only to determine the type and is not referenced ! integer, intent(in) :: array_size - ${t1}$, intent(in) :: mold + ${t1}$, intent(in) #{if t1 == 'real(dp)'}#, optional #{endif}#:: mold ${t1}$ :: res(array_size) res = rvs_norm_array_${t1[0]}$${k1}$ (0.0_${k1}$, 1.0_${k1}$, array_size) diff --git a/test/stats/test_distribution_normal.fypp b/test/stats/test_distribution_normal.fypp index 8e959fab1..a20a569cc 100644 --- a/test/stats/test_distribution_normal.fypp +++ b/test/stats/test_distribution_normal.fypp @@ -164,14 +164,19 @@ contains call random_seed(seed, get) ! default mold form: mold used only to disambiguate kind - #:if t1[0] == "r" - mold = 0.0_${k1}$ + ! For real(dp), mold is optional; for other types (including complex), it's required + #:if t1[0] == "r" and k1 == "dp" + a2 = nor_rvs(10) ! mold optional for rdp only, defaults to real(dp) #:else - mold = (0.0_${k1}$, 0.0_${k1}$) + #! mold required for all other types including complex and non-dp kinds + #:if t1[0] == "r" + mold = 0.0_${k1}$ + #:else + mold = (0.0_${k1}$, 0.0_${k1}$) + #:endif + a2 = nor_rvs(10, mold) #: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}$