diff --git a/src/math.f90 b/src/math.f90
index 191791ee2..36fe91bac 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -1068,6 +1068,7 @@ integer pure function math_factorial(n)
 
   integer, intent(in) :: n
 
+
   math_factorial = product(math_range(n))
 
 end function math_factorial
@@ -1081,6 +1082,7 @@ integer pure function math_binomial(n,k)
   integer, intent(in) :: n, k
   integer :: i, k_, n_
 
+
   k_ = min(k,n-k)
   n_ = 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
 !--------------------------------------------------------------------------------------------------
-integer pure function math_multinomial(alpha)
+integer pure function math_multinomial(k)
 
-  integer, intent(in), dimension(:) :: alpha
+  integer, intent(in), dimension(:) :: k
   integer :: i
 
-  math_multinomial = 1
-  do i = 1, size(alpha)
-    math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
-  enddo
+
+  math_multinomial = product([(math_binomial(sum(k(1:i)),k(i)),i=1,size(k))])
 
 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,3) :: m
 
+
   m(1:3,1) = v1-v2
   m(1:3,2) = v1-v3
   m(1:3,3) = v1-v4
@@ -1289,6 +1290,9 @@ subroutine selfTest
   if(math_binomial(49,6) /= 13983816) &
     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))
   if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
     error stop 'math_LeviCivita(even)'