223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

multinomial tested+simplified

This commit is contained in:
Martin Diehl 2021-07-22 15:31:30 +02:00
parent d9aa638ad7
commit 4fab078d18
1 changed files with 10 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/math.f90
View File

@ -1068,6 +1068,7 @@ integer pure function math_factorial(n)
integer, intent(in) :: n integer, intent(in) :: n
math_factorial = product(math_range(n)) math_factorial = product(math_range(n))
end function math_factorial end function math_factorial
@ -1081,6 +1082,7 @@ integer pure function math_binomial(n,k)
integer, intent(in) :: n, k integer, intent(in) :: n, k
integer :: i, k_, n_ integer :: i, k_, n_
k_ = min(k,n-k) k_ = min(k,n-k)
n_ = n n_ = n
math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n
@ -1095,15 +1097,13 @@ end function math_binomial
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief multinomial coefficient !> @brief multinomial coefficient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer pure function math_multinomial(alpha) integer pure function math_multinomial(k)
integer, intent(in), dimension(:) :: alpha integer, intent(in), dimension(:) :: k
integer :: i integer :: i
math_multinomial = 1
do i = 1, size(alpha) math_multinomial = product([(math_binomial(sum(k(1:i)),k(i)),i=1,size(k))])
math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
enddo
end function math_multinomial end function math_multinomial
@ -1116,6 +1116,7 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4)
real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4 real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4
real(pReal), dimension (3,3) :: m real(pReal), dimension (3,3) :: m
m(1:3,1) = v1-v2 m(1:3,1) = v1-v2
m(1:3,2) = v1-v3 m(1:3,2) = v1-v3
m(1:3,3) = v1-v4 m(1:3,3) = v1-v4
@ -1289,6 +1290,9 @@ subroutine selfTest
if(math_binomial(49,6) /= 13983816) & if(math_binomial(49,6) /= 13983816) &
error stop 'math_binomial' error stop 'math_binomial'
if(math_multinomial([1,2,3,4]) /= 12600) &
error stop 'math_multinomial'
ijk = cshift([1,2,3],int(r*1.0e2_pReal)) ijk = cshift([1,2,3],int(r*1.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
error stop 'math_LeviCivita(even)' error stop 'math_LeviCivita(even)'