#! /usr/bin/env python
#
def commul ( n, nfactor, iarray ):

#*****************************************************************************80
#
## COMMUL computes a multinomial combinatorial coefficient.
#
#  Discussion:
#
#    The multinomial coefficient is a generalization of the binomial
#    coefficient.  It may be interpreted as the number of combinations of
#    N objects, where IARRAY(1) objects are indistinguishable of type 1,
#    ... and IARRAY(K) are indistinguishable of type NFACT.
#
#    The formula is
#
#      NCOMB = N! / ( IARRAY(1)! IARRAY(2)! ... IARRAY(NFACT)! )
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 January 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, determines the numerator.
#
#    Input, integer NFACTOR, the number of factors in the numerator.
#
#    Input, integer IARRAY(NFACTOR).
#    IARRAY contains the NFACT values used in the denominator.
#    Note that the sum of these entries should be N,
#    and that all entries should be nonnegative.
#
#    Output, integer NCOMB, the value of the multinomial coefficient.
#
  import numpy as np
  from r8_gamma_log import r8_gamma_log
  from sys import exit

  for i in range ( 0, nfactor ):

    if ( iarray[i] < 0 ):
      print ( '' )
      print ( 'COMMUL - Fatal error!' )
      print ( '  Entry %d of IARRAY = %d' % ( i, iarray[i] ) )
      print ( '  But this value must be nonnegative.' )
      exit ( 'COMMUL - Fatal error!' )

  isum = np.sum ( iarray )

  if ( isum != n ):
    print ( '' )
    print ( 'COMMUL - Fatal error!' )
    print ( '  The sum of the IARRAY entries is %d' % (isum ) )
    print ( '  But it must equal N = %d' % ( n ) )
    exit ( 'COMMUL - Fatal error!' )
 
  facn = r8_gamma_log ( float ( n + 1 ) )
 
  for i in range ( 0, nfactor ):
 
    fack = r8_gamma_log ( float ( iarray[i] + 1 ) )
    facn = facn - fack

  ncomb = int ( np.round ( np.exp ( facn ) ) )

  return ncomb

def commul_test ( ):

#*****************************************************************************80
#
## COMMUL_TEST tests COMMUL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    26 January 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'COMMUL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  COMMUL computes a multinomial coefficient.' )

  n = 8
  nfactor = 2
  factor = np.array ( [ 6, 2 ] )
  ncomb = commul ( n, nfactor, factor )
  print ( '' )
  print ( '  N = %d' % ( n ) )
  print ( '  Number of factors = %d' % ( nfactor ) )
  for i in range ( 0, nfactor ):
    print ( '  %2d  %2d' % ( i, factor[i] ) )
  print ( '  Value of coefficient = %d' % ( ncomb ) )

  n = 8
  nfactor = 3
  factor = np.array ( [ 2, 2, 4 ] )
  ncomb = commul ( n, nfactor, factor )
  print ( '' )
  print ( '  N = %d' % ( n ) )
  print ( '  Number of factors = %d' % ( nfactor ) )
  for i in range ( 0, nfactor ):
    print ( '  %2d  %2d' % ( i, factor[i] ) )
  print ( '  Value of coefficient = %d' % ( ncomb ) )

  n = 13
  nfactor = 4
  factor = np.array ( [ 5, 3, 3, 2 ] )
  ncomb = commul ( n, nfactor, factor )
  print ( '' )
  print ( '  N = %d' % ( n ) )
  print ( '  Number of factors = %d' % ( nfactor ) )
  for i in range ( 0, nfactor ):
    print ( '  %2d  %2d' % ( i, factor[i] ) )
  print ( '  Value of coefficient = %d' % ( ncomb ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'COMMUL_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  commul_test ( )
  timestamp ( )