#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
module stdlib_stats_distribution_exponential
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
use stdlib_kinds, only : sp, dp, xdp, qp, int32
use stdlib_random, only : dist_rand
use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform
implicit none
private
integer :: ke(0:255)
real(dp) :: we(0:255), fe(0:255)
logical :: zig_exp_initialized = .false.
public :: rvs_exp
public :: pdf_exp
public :: cdf_exp
interface rvs_exp
!! Version experimental
!!
!! Exponential Distribution Random Variates
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
!! rvs_exp-exponential-distribution-random-variates))
!!
module procedure rvs_exp_0_rsp !0 dummy variable
#:for k1, t1 in RC_KINDS_TYPES
module procedure rvs_exp_${t1[0]}$${k1}$ !1 dummy variable
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
module procedure rvs_exp_array_${t1[0]}$${k1}$ !2 dummy variables
#:endfor
end interface rvs_exp
interface pdf_exp
!! Version experimental
!!
!! Exponential Distribution Probability Density Function
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
!! pdf_exp-exponential-distribution-probability-density-function))
!!
#:for k1, t1 in RC_KINDS_TYPES
module procedure pdf_exp_${t1[0]}$${k1}$
#:endfor
end interface pdf_exp
interface cdf_exp
!! Version experimental
!!
!! Exponential Cumulative Distribution Function
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
!! cdf_exp-exponential-distribution-cumulative-distribution-function))
!!
#:for k1, t1 in RC_KINDS_TYPES
module procedure cdf_exp_${t1[0]}$${k1}$
#:endfor
end interface cdf_exp
contains
impure subroutine zigset
! Marsaglia & Tsang generator for random normals & random exponentials.
! Translated from C by Alan Miller (amiller@bigpond.net.au)
!
! 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 :: M2 = 2147483648.0_dp, ve = 0.003949659822581572_dp
real(dp), parameter :: ONE = 1.0_dp
real(dp) :: de, te, q
integer :: i
de = 7.697117470131487_dp
te = de
! tables for random exponentials
q = ve * exp(de)
ke(0) = int((de / q) * M2, kind = int32)
ke(1) = 0
we(0) = q / M2
we(255) = de / M2
fe(0) = ONE
fe(255) = exp(- de)
do i = 254, 1, -1
de = -log(ve / de + exp(- de))
ke(i+1) = int(M2 * (de / te), kind = int32)
te = de
fe(i) = exp(- de)
we(i) = de / M2
end do
zig_exp_initialized = .true.
end subroutine zigset
#:for k1, t1 in REAL_KINDS_TYPES
impure function rvs_exp_0_${t1[0]}$${k1}$( ) result(res)
!
! Standard exponential random variate (lambda=1)
!
${t1}$ :: res, x
${t1}$, parameter :: r = 7.69711747013104972_${k1}$
integer :: jz, iz
if(.not. zig_exp_initialized ) call zigset
iz = 0
jz = dist_rand(1_int32) ! 32bit random integer
iz = iand( jz, 255 ) ! random integer in [0, 255]
if( abs( jz ) < ke(iz) ) then
res = abs(jz) * we(iz)
else
L1: do
if( iz == 0 ) then
res = r - log( uni(1.0_${k1}$) )
exit L1
end if
x = abs( jz ) * we(iz)
if(fe(iz) + uni(1.0_${k1}$) * (fe(iz-1) - fe(iz)) < exp(-x)) then
res = x
exit L1
end if
jz = dist_rand(1_int32)
iz = iand( jz, 255 )
if( abs( jz ) < ke(iz) ) then
res = abs( jz ) * we(iz)
exit L1
end if
end do L1
endif
end function rvs_exp_0_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res)
!
! Exponential distributed random variate
!
${t1}$, intent(in) :: lambda
${t1}$ :: res
if (lambda <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = rvs_exp_0_${t1[0]}$${k1}$( )
res = res / lambda
end if
end function rvs_exp_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res)
${t1}$, intent(in) :: lambda
${t1}$ :: res
real(${k1}$) :: tr, ti
tr = rvs_exp_r${k1}$(lambda % re)
ti = rvs_exp_r${k1}$(lambda % im)
res = cmplx(tr, ti, kind=${k1}$)
end function rvs_exp_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res)
${t1}$, intent(in) :: lambda
integer, intent(in) :: array_size
${t1}$ :: res(array_size), x, re
${t1}$, parameter :: r = 7.69711747013104972_${k1}$
integer :: jz, iz, i
if (lambda <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if
if(.not. zig_exp_initialized) call zigset
do i = 1, array_size
iz = 0
jz = dist_rand(1_int32)
iz = iand( jz, 255 )
if( abs( jz ) < ke(iz) ) then
re = abs(jz) * we(iz)
else
L1: do
if( iz == 0 ) then
re = r - log( uni(1.0_${k1}$) )
exit L1
end if
x = abs( jz ) * we(iz)
if(fe(iz) + uni(1.0_${k1}$)*(fe(iz-1)-fe(iz)) < exp(-x)) then
re = x
exit L1
end if
jz = dist_rand(1_int32)
iz = iand( jz, 255 )
if( abs( jz ) < ke(iz) ) then
re = abs( jz ) * we(iz)
exit L1
end if
end do L1
endif
res(i) = re / lambda
end do
end function rvs_exp_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res)
${t1}$, intent(in) :: lambda
integer, intent(in) :: array_size
${t1}$ :: res(array_size)
integer :: i
real(${k1}$) :: tr, ti
do i = 1, array_size
tr = rvs_exp_r${k1}$(lambda % re)
ti = rvs_exp_r${k1}$(lambda % im)
res(i) = cmplx(tr, ti, kind=${k1}$)
end do
end function rvs_exp_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
!
! Exponential Distribution Probability Density Function
!
${t1}$, intent(in) :: x, lambda
real(${k1}$) :: res
if (lambda <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else if (x < 0._${k1}$) then
res = 0._${k1}$
else
res = exp(- x * lambda) * lambda
end if
end function pdf_exp_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
${t1}$, intent(in) :: x, lambda
real(${k1}$) :: res
res = pdf_exp_r${k1}$(x % re, lambda % re)
res = res * pdf_exp_r${k1}$(x % im, lambda % im)
end function pdf_exp_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
!
! Exponential Distribution Cumulative Distribution Function
!
${t1}$, intent(in) :: x, lambda
real(${k1}$) :: res
if (lambda <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else if (x < 0._${k1}$) then
res = 0._${k1}$
else
res = 1.0_${k1}$ - exp(- x * lambda)
end if
end function cdf_exp_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
${t1}$, intent(in) :: x, lambda
real(${k1}$) :: res
res = cdf_exp_r${k1}$(x % re, lambda % re)
res = res * cdf_exp_r${k1}$(x % im, lambda % im)
end function cdf_exp_${t1[0]}$${k1}$
#:endfor
end module stdlib_stats_distribution_exponential
