Compute the multinomial coefficient

27

Time for another easy challenge in which all can participate!

The multinomial theorem states: Formula for computing the nth power of a multinomial

The expression in parentheses is the multinomial coefficient, defined as:

Multinomial coefficient

Allowing the terms ki to range over all integer partitions of n gives the n-th level of Pascal's m-simplex. Your task is to compute this coefficient.

Task

Write a program or function which takes m numbers, n, k1, k2, ... ,km-1, and outputs or returns the corresponding multinomial coefficient. Your program may optionally take m as an additional argument if need be. Note that km is not in the input.

  • These numbers may be input in any format one likes, for instance grouped into lists or encoded in unary, or anything else, as long as the actual computation of the multinomial coefficient is performed by your code, and not the encoding process.

  • Output format is similarly flexible.

  • All code should run in less than one minute for n and m up to 1000.

  • Do not worry about integer overflow.

  • Built-ins designed to compute the multinomial coefficient are not allowed.

  • Standard loopholes apply.

Scoring

This is code golf: Shortest solution in bytes wins.

Test cases

Input: 3, [2, 0]
Output: 3

Input: 3, [1, 1]
Output: 6

Input: 11, [1, 4, 4]
Output: 34650

Input: 4, [1,2]
Output: 12

Input: 15, [5,4,3,2]
Output: 37837800

Input: 95, [65,4,4]
Output: 1934550571913396675776550070308250

Input: 32, [2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]
Output: 4015057936610313875842560000000

Input: 15, [3,3,3,3]
Output: 168168000

Input: 1000, [10,10,10,10,10,10,10,10,10,10,100,100,100,100,100,100,100,100]
Output: 1892260836114766064839886173072628322819837473493540916521650371620708316292211493005889278395285403318471457333959691477413845818795311980925098433545057962732816261282589926581281484274178579110373517415585990780259179555579119249444675675971136703240347768185200859583936041679096016595989605569764359198616300820217344233610087468418992008471158382363562679752612394898708988062100932765563185864346460326847538659268068471585720069159997090290904151003744735224635733011050421493330583941651019570222984959183118891461330718594645532241449810403071583062752945668937388999711726969103987467123014208575736645381474142475995771446030088717454857668814925642941036383273459178373839445456712918381796599882439216894107889251444932486362309407245949950539480089149687317762667940531452670088934094510294534762190299611806466111882595667632800995865129329156425174586491525505695534290243513946995156554997365435062121633281021210807821617604582625046557789259061566742237246102255343862644466345335421894369143319723958653232683916869615649006682399919540931573841920000000000000

Input: 33, [17]
Output: 1166803110

Input: 55, [28]
Output: 3824345300380220

quintopia

Posted 2016-01-14T14:39:55.137

Reputation: 3 899

Can we have imprecision errors? I.e., instead of 1934550571913396675776550070308250, can we output 1.9345505719133966e+33? – Conor O'Brien – 2016-01-14T15:09:53.070

@CᴏɴᴏʀO'Bʀɪᴇɴ If you used 64-bit floats, you won't be able to represent input [1000 {999 ones}] at all, because the exponent is way beyond what 64-bit floats can represent. (128-bit floats will probably suffice, but I'm assuming you want to use JavaScript's native number type?) – Martin Ender – 2016-01-14T15:15:07.493

@MartinBüttner Yes, that is a correct assumption. – Conor O'Brien – 2016-01-14T15:28:49.277

2@quintopia "Time for another easy challenge in which all can participate!". Everyone except me! (Since I have no idea what Pascals simplex and multinomials are D:) LOL. – Ashwin Gupta – 2016-01-14T16:56:19.933

@AshwinGupta Don't worry about it. You just compute the expression in the second image and you're good to go! – quintopia – 2016-01-14T17:01:55.370

@quintopia hey thanks :D! I'll give it a shot. Sorry I didn't respond all day I was at school then had an after school thing. – Ashwin Gupta – 2016-01-15T03:15:40.073

Answers

21

Jelly, 7 6 bytes

;_/!:/

Look ma, no Unicode! This program takes a single list as input, with n at its first index.

Try it online! or verify all test cases at once.

How it works

;_/!:/ Input: A (list)

 _/    Reduce A by subtraction. This subtracts all other elements from the first.
;      Concatenate A with the result to the right.
   !   Apply factorial to all numbers in the resulting list.
    :/ Reduce the result by division. This divides the first element by the others.

Dennis

Posted 2016-01-14T14:39:55.137

Reputation: 196 637

This is pretty much the algorithm I had in mind as being the simplest. – quintopia – 2016-01-14T16:58:10.840

9

CJam, 11 bytes

l~_:-+:m!:/

Input as a single list with n first:

[95 65 4 4]

This handles inputs up to n and m 1000 pretty much instantly.

Test it here.

Explanation

l~  e# Read a line of input and evaluate it.
_   e# Duplicate.
:-  e# Fold subtraction over the list. A fold is essentially a foreach loop that starts
    e# from the second element. Hence, this subtracts all the k_i from n, giving k_m.
+   e# Append k_m to the list.
:m! e# Compute the factorial of each element in the list.
:/  e# Fold division over the list. Again, this divides n! by each of the k_i!.

Martin Ender

Posted 2016-01-14T14:39:55.137

Reputation: 184 808

It looks like you will actually lose the byte-count competition, but I have to say I am impressed with CJam's insane terseness. – phord – 2016-01-14T15:39:54.880

@phord Well CJam is no match for Jelly (or Pyth for that matter). But I was quite surprised myself how compact it ended up. My first solution had 21 bytes, and while it didn't seem optimal, I didn't think I could almost cut that in half. – Martin Ender – 2016-01-14T15:41:20.100

4

MATL, 21 15 bytes

Let's put the log-gamma function to good use. This avoids internal overflowing by working with logarithms of factorials, not with factorials themselves.

1+ZgiO$Gs-h1+Zgs-ZeYo

This works in current version (9.2.2) of the language/compiler, which is earlier than this challenge.

Inputs are: first a number, then a numeric vector. The result is produced as a double, which limits maximum output to somewhere around 2^52.

Example

>> matl 1+ZgiO$Gs-h1+Zgs-ZeYo
> 15
> [5 4 3 2]
37837800

Explanation

1+       % implicit input (number). Add 1
Zg       % log-gamma function
i        % input (numeric vector).
0$G      % push both inputs
s-       % sum the second input (vector) and subtract from first
h1+      % append to vector. Add 1
Zg       % log-gamma function, element-wise on extended vector
s        % sum of results
-        % subtract from previous result of log-gamma
Ze       % exponential
Yo       % round. Implicit display

Luis Mendo

Posted 2016-01-14T14:39:55.137

Reputation: 87 464

4Try it online! now has experimental MATL support: http://matl.tryitonline.net/#code=MStaZ2lPJEdzLWgxK1pncy1aZVlv&input=MTUKWzUgNCAzIDJdCg Suggestions are welcome. – Dennis – 2016-01-14T19:03:26.837

1@Dennis Hey! What a surprise!!! How can I thank you?? I do have a suggestion: if you ever come to Madrid I owe you a good dinner and some drinks – Luis Mendo – 2016-01-14T19:24:24.387

I'm really grateful. It's great to have it online. How will we handle revisions? I'm still constantly updating the language, you know... – Luis Mendo – 2016-01-14T19:25:06.893

For now, I'm manually updating the interpreters. If you make an update, just ping me in The Nineteenth Byte and I'll pull it asap. -- I'll have to go to Madrid in the near future, so I'll keep your offer in mind. ;) – Dennis – 2016-01-14T19:36:04.517

@Dennis Great! That way we can meet in person! – Luis Mendo – 2016-01-14T21:09:29.617

4

PowerShell, 91 74 bytes

Woo! My 100th answer on PPCG!

param($n,$k)(1..$n-join'*'|iex)/(($k|%{$n-=$_;1..$_})+(1..$n)-join'*'|iex)

Whew. Not going to win shortest-code, that's for sure. Uses a couple neat tricks with ranges, though. And this is probably complete gibberish to anyone not familiar with PowerShell.

Explanation

First we take input with param($n,$k) and expect $k to be an array, e.g. .\compute-the-multinomial-coefficient.ps1 11 @(1,4,4).

We'll start with the numerator (everything to the left of /). That's simply a range from 1..$n that has been -joined together with * and then evaluated with iex to compute the factorial (i.e., 1*2*3*...*$n).

Next, we loop over $k|%{...} and each iteration we subtract the current value $_ from $n (which we don't care about anymore) to formulate $k_m later. Additionally, we generate the range 1..$k_i each iteration, which gets left on the pipeline. Those pipeline objects get array-concatenated with the second expression, range 1..$n (which is $k_m at this point). All of that is finally -joined together with * and evaluated with iex, similar to the numerator (this works because x! * y! = 1*2*3*...*x * 1*2*3*...*y, so we don't care about individual ordering).

Finally, the / happens, the numerator is divided by the denominator, and output.

Handles output correctly for larger numbers, since we're not explicitly casting any variables as any particular datatypes, so PowerShell will silently re-cast as different datatypes on the fly as needed. For the larger numbers, outputs via scientific notation to best preserve significant figures as the datatypes get re-cast. For example, .\compute-the-multinomial-coefficient.ps1 55 @(28) will output 3.82434530038022E+15. I'm presuming this to be OK given "Output format is similarly flexible" is specified in the challenge and quintopia's comments "If the final result can fit within the natively supported integer types, then the result must be accurate. If it cannot, there is no restriction on what may be output."


Alternatively

Depending upon output formatting decisions, the following at 92 bytes

param($n,$k)((1..$n-join'*'|iex)/(($k|%{$n-=$_;1..$_})+(1..$n)-join'*'|iex)).ToString('G17')

Which is the same as the above, just uses explicit output formatting with .ToString('G17') to achieve the desired number of digits. For 55 @(28) this will output 3824345300380220.5


Edit1 -- Saved 17 bytes by getting rid of $d and just calculating it directly, and getting rid of the calculation for $k_m by stringing it along while we loop $k
Edit2 -- Added alternate version with explicit formatting

AdmBorkBork

Posted 2016-01-14T14:39:55.137

Reputation: 41 581

3

APL (Dyalog Extended), 9 bytes

×/2!/+\⍛,

Try it online!

Using the idea from my APL answer on another challenge that involves multinomials.

A tacit function whose left argument is the list of k's, and right argument is n. The test cases check if it agrees with Adam's solution with left and right arguments flipped.

How it works

×/2!/+\⍛,
     +\    ⍝ Cumulative sum of k's (up to m-1'th element)
       ⍛,  ⍝ Append n (sum of k_1 to k_m)
  2!/      ⍝ Binomial of consecutive pairs
×/         ⍝ Product

$$ \frac{(k_1 + k_2 + \cdots + k_m)!}{k_1! k_2! \cdots k_m!} = \frac{(k_1 + k_2)!}{k_1! k_2!} \times \frac{(k_1 + k_2 + \cdots + k_m)!}{(k_1 + k_2)! k_3! \cdots k_m!} $$

$$ = \frac{(k_1 + k_2)!}{k_1! k_2!} \times \frac{(k_1 + k_2 + k_3)!}{(k_1 + k_2)! k_3!} \times \frac{(k_1 + k_2 + \cdots + k_m)!}{(k_1 + k_2 + k_3)! \cdots k_m!} $$

$$ = \cdots = \binom{k_1 + k_2}{k_1} \binom{k_1 + k_2 + k_3}{k_1 + k_2} \cdots \binom{k_1 + \cdots + k_m}{k_1 + \cdots + k_{m-1}} $$

Bubbler

Posted 2016-01-14T14:39:55.137

Reputation: 16 616

2

APL (Dyalog Unicode), 16 bytesSBCS

Entirely based on the mathematical skill of my colleague Marshall.

Anonymous infix function. Takes k as right argument and n as left argument.

{×/⍵!⍺-+\¯1↓0,⍵}

Try it online!

{} anonymous lambda; is left argument (n) and is right argument (k)

0,⍵ prepend a zero to k

¯1↓ drop the last item from that

+\ cumulative sum of that

⍺- subtract that from n

⍵! (k) that

×/ product of that

Adám

Posted 2016-01-14T14:39:55.137

Reputation: 37 779

2

Python 3, 93 91

Thanks to Dennis and FryAmTheEggman.

f=lambda x:0**x or x*f(x-1)
def g(n,k):
    r=f(n)
    for i in k:r//=f(i)
    return r//f(n-sum(k))

n as integer, k as iterable.

Ungolfed:

import functools #cache

@functools.lru_cache(maxsize=None) #cache results to speed up calculations
def factorial(x):
    if x <= 1: return 1
    else: return x * factorial(x-1)

def multinomial(n, k):
    ret = factorial(n)
    for i in k: ret //= factorial(i)
    km = n - sum(k)
    return ret//factorial(km)

Trang Oul

Posted 2016-01-14T14:39:55.137

Reputation: 656

1You can use a single space instead of four for the dynamic whitespace bit – Conor O'Brien – 2016-01-14T15:31:54.570

I used tabs, they got replaced in this post. Byte count seems to be OK. I'm not sure about float result and possible overflow. – Trang Oul – 2016-01-14T15:47:32.827

@TrangOul Oh, I thought / would yield an integer if the result was representable as an integer. It doesn't though. In that case, you're right. – Martin Ender – 2016-01-14T15:54:23.177

2>

  • This produces an incorrect for 95, [65, 4, 4]. Note that the input does not contain k_m. 2. You don't seem to be using from functools import* at all.
  • < – Dennis – 2016-01-14T16:01:59.503

    @Dennis: fixed. I need functools for reduce.

    – Trang Oul – 2016-01-14T20:01:23.140

    2>

  • Your golfed code doesn't use reduce. 2. import math;f=math.factorial saves a byte. 3. Python 2 would allow you to get rid of the second / in //.
  • < – Dennis – 2016-01-14T20:07:01.237

    Yeah, I realized that after writing previous comment... Sorry about that. Note to self: don't play with code golf with wife around babbling :/ . – Trang Oul – 2016-01-14T20:15:37.340

    1

    Defining f on your own saves some bytes: f=lambda x:0**x or x*f(x-1).

    – FryAmTheEggman – 2016-01-14T20:21:41.777

    Oh, I didn't realize that in Python 0**0 == 1. Thanks. – Trang Oul – 2016-01-15T07:08:17.400

    2

    Mathematica, 26 bytes

    #!/Times@@({#-+##2,##2}!)&
    

    Example:

    In[1]:= #!/Times@@({#-+##2,##2}!)&[95,65,4,4]
    
    Out[1]= 1934550571913396675776550070308250
    

    alephalpha

    Posted 2016-01-14T14:39:55.137

    Reputation: 23 988

    1

    Haskell, 59 58 bytes

    f n=product[1..n]
    n#k=foldl((.f).div)(f n)k`div`f(n-sum k)
    

    Try it online!

    Thanks to BMO for saving 1 byte!

    Cristian Lupascu

    Posted 2016-01-14T14:39:55.137

    Reputation: 8 369

    1

    05AB1E, 8 bytes

    Ƹ«!R.«÷
    

    Try it online! Explanation:

    Æ           Subtract all the elements from the first
     ¸«         Append to the original list
       !        Take the factorial of all the elements
        R.«÷    Reduce by integer division
    

    I can't seem to find better ways of performing step 2 or step 4.

    Neil

    Posted 2016-01-14T14:39:55.137

    Reputation: 95 035

    1

    APL (Dyalog Unicode), 17 bytes

    ⊢∘!÷(×/∘!⊣,⊢-1⊥⊣)
    

    Try it online!

    Infix tacit function (Thanks to @Adám for the 2 bytes it saves.)


    APL (Dyalog Unicode), 19 bytes

    {(!⍺)÷×/!(⍺-+/⍵),⍵}
    

    Try it online!

    Infix Dfn.

    Both functions above compute the given formula.

    J. Sallé

    Posted 2016-01-14T14:39:55.137

    Reputation: 3 233

    1

    PARI/GP, 43 bytes

    Pretty straightforward; aside from formatting, the ungolfed version might well be identical.

    m(n,v)=n!/prod(i=1,#v,v[i]!)/(n-vecsum(v))!
    

    Charles

    Posted 2016-01-14T14:39:55.137

    Reputation: 2 435

    1

    Matlab 48 bytes

    You need to set format to long in advance to get the higher precision. Then, it's quite straightforward:

    @(n,k)factorial(n)/prod(factorial([k,n-sum(k)]))
    
    ans(95, [65,4,4])
    ans =
    
     1.934550571913395e+33
    

    brainkz

    Posted 2016-01-14T14:39:55.137

    Reputation: 349

    1

    Pyth, 10 bytes

    /F.!MaQ-FQ
    

    Try it online: Demonstration

    Explanation:

    /F.!MaQ-FQ   implicit: Q = input list
           -FQ   reduce Q by subtraction
         aQ      append the result to Q
      .!M        compute the factorial for each number
    /F           reduce by division
    

    Jakube

    Posted 2016-01-14T14:39:55.137

    Reputation: 21 462

    1

    J, 16 bytes

    [(%*/)&:!],(-+/)
    

    Usage

    For larger values, a suffix of x is used to denote extended precision integers.

       f =: [(%*/)&:!],(-+/)
       11 f 1 4 4
    34650
       15x f 5 4 3 2
    37837800
    

    Explanation

    [(%*/)&:!],(-+/)  Input: n on LHS, A on RHS
                 +/   Reduce A using addition
                -     Subtract that sum from n, this is the missing term
             ]        Get A
              ,       Append the missing term to A to make A'
    [                 Get n
          &:!         Take the factorial of n and each value in A'
       */             Reduce using multiplication the factorials of A'
      %               Divide n! by that product and return
    

    miles

    Posted 2016-01-14T14:39:55.137

    Reputation: 15 654

    0

    Perl 6,  52  50 bytes

    ->\n,\k{[*](1..n)div[*] ([*] 1..$_ for |k,[-] n,|k)}
    

    Test it

    ->\n,\k{[*](1..n)/[*] ([*] 1..$_ for |k,[-] n,|k)}
    

    Test it (result is a Rational with denominator of 1)

    Expanded:

    ->     # pointy block lambda
      \n,
      \k
    {
        [*]( 1 .. n )   # factorial of 「n」
    
      /                 # divide (produces Rational)
    
        [*]             # reduce the following using &infix:«*»
    
          (
              [*] 1..$_ # the factorial of
    
            for         # each of the following
    
              |k,       # the values of 「k」 (slipped into list)
              [-] n,|k  # 「n」 minus the values in 「k」
          )
    }
    

    Brad Gilbert b2gills

    Posted 2016-01-14T14:39:55.137

    Reputation: 12 713

    0

    Clojure, 70 bytes

    #(let[a apply](a /(map(fn[x](a *(map inc(range x))))(conj %(a - %)))))
    

    Creates an anonymous function taking all of the arguments as a single list, with n first.

    30 characters are "wasted" just defining the damn factorial function. Oh well.

    MattPutnam

    Posted 2016-01-14T14:39:55.137

    Reputation: 521