# Expected number of rounds for this labeling scheme

Here is an interesting math problem:

Let's say that there are $$\n\$$ indistinguishable unlabeled objects in a bin. For every "round", pull $$\k\$$ objects randomly out of the bin with equal probability, and apply a label on each object (Nothing happens to objects already having a label). Then put these $$\k\$$ objects, now all labelled, back into the bin and repeat the process until all objects have a label applied to them. What is the expected number of rounds needed to label all the objects, given $$\n\$$ and $$\k\$$?

Your task is to write a program that solves this exact math problem. Turns out, this is quite hard to calculate by hand and is better left for computers to bash out. I will list three methods of solving this problem below. You do not have to follow any of these three methods.

# Methods

### Method 1

I have asked this question on Math SE recently as I was also stumped on how to do this problem. One answer suggests the following formula: $$\sum_{j=1}^n (-1)^{j+1} \binom{n}{j} \frac{1}{1-\binom{n-j}{k}/\binom{n}{k}}$$ I have implemented this formula in Python: Try it online!

### Method 2

Another answer from the same question provides the following formula: $${\bf e_1}^T \left(I_n - \left({\mathcal{1}}_{j\geq i} {\frac{ {{n-i \choose j-i}} {i \choose k-j+i} }{n\choose k}} \right)_{i=0,j=0}^{n-1} \right)^{-1} {\bf 1}$$ where "$$\I_n\$$ is the $$\n \times n\$$ identity matrix, $$\\bf e_1\$$ is the first standard basis vector of $$\\mathbb{R}^n\$$ and $$\{\bf1}\$$ is the $$\n\$$-vector of all ones. And in the definition of the other matrix $$\\mathcal{1}\$$ is the indicator function, i.e. the $$\i,j\$$-element is $$\0\$$ if $$\j and otherwise that quotient of binomials...," quoting from the linked answer.

I have implemented this formula in Mathematica: Try it online!

### Method 3

(This is very likely not the golfiest way for any language, but this is what I used to generate the larger test cases.)

Another way is to utilize the recursive nature of this problem, as the results of each round depends on the results of all the previous rounds. The recursive formula is as follows ($$\E\$$ is the function that actually calculates the expected value):

\begin{align}E(n,k)&=f(n,k,n)\\f(n,k,u)&=\begin{cases}0&\text{if }u=0\\\frac{\displaystyle\binom nk+\sum_{i=1}^{\min(u,k)}{\binom ui\binom{n-u}{k-i}}f(n,k,u-i)}{\displaystyle\binom nk-\binom{n-u}k}&\text{otherwise}\end{cases}\end{align}

I have implemented this formula in Python: Try it online!

# I/O format and other rules

Your program or function should take in two positive integer values $$\n,k\$$ with $$\n\ge k\$$ and output the expected value corresponding to those values. The output can be a floating point value with an absolute/relative error under $$\10^{-6}\$$ from the true value, though the algorithm should theoretically work for arbitrarily large $$\n\$$ and $$\k\$$. In practice, it is acceptable if the program is limited by time, memory, or data-type size. This means that although your program may output inaccurate floating point values for certain inputs (as long as the outputs are within the allowed error as mentioned earlier), it should still theoretically be calculating the exact value of the true output.

Because the expected value is always going to be a rational number, you can also output as an exact fraction (does not need to be simplified and can be improper), if your language supports such a data type, or if you want to output a pair of numbers, the numerator and the denominator. If in doubt, always refer back to the default I/O methods.

Your program should also be deterministic, meaning it should always produce the same value when ran. Namely, no "running a gazillion simulations to get the experimental expected value."

# Test Cases

n,k -> output
1,1 -> 1
4,2 -> 3.8
6,2 -> 6.928571428571428571428571429
7,4 -> 3.779411764705882352941176471
7,7 -> 1
8,1 -> 21.74285714285714285714285714
8,5 -> 3.458181818181818181818181818
10,3 -> 9.046212999306305594338048699
12,7 -> 4.212148576244657154897357173
20,2 -> 35.30473485789483826080766262
100,7 -> 72.28020367501059456930677280
1000,69 -> 105.2447804224750397369209265
1256,5 -> 1934.893442395300917211652090

This is , so the shortest program in bytes wins!

• Sandbox Oct 24 at 3:36
• It seems that the output for 4,2 is 3.8 using any of the three methods. Oct 24 at 3:51
• Your program should also be deterministic, meaning it should always produce the same value when ran. Namely, no "running a gazillion simulations to get the experimental expected value." so is a deterministic simulation allowed?
– l4m2
Oct 24 at 3:57
• @alephalpha Ah yep, I accidentally put the output for 4,3. Will fix Oct 24 at 4:26
• @l4m2 The program should be theoretically correct for any pairs of $n,k$, so you would need to prove that your simulation would give the correct output for every single case. Oct 24 at 4:29

# 05AB1E, 13 10 bytes

-3 thanks to @alphalpha

Ýsc¤/<z¹F¥

Attempt This Online!

Uses a variant of the first method. Ports this Jelly answer for the binomial transform.

Ý      calculate [0,1,...,n]
s      swap so k is on top
c      calculate [0Ck, 1Ck, ..., nCk]
¤      push the last item without poping
/      divide by it
<      subtract one
z      and inverse. Luckily for us, 05AB1E have 1/0 = 0
¹F¥    calculate the binomial transform, by calculating the differences n times.
• Fü-}Ä -> Oct 24 at 6:52

# Charcoal, 48 42 bytes

ＮθＮηＩΣＥθ∕Π⁺±⊕…ιθ⁻θ…⁰η×Π…·¹⁻θι⁻Π⁻ι…⁰ηΠ⁻θ…⁰η

Attempt This Online! Link is to verbose version of code. Explanation: Implements a modification of a formula used by several other answers:

$$\sum_{i=1}^n\frac{-i\choose 1-i+n}{\left(\frac{i-1\choose k}{n\choose k}-1\right)}$$

$$= \sum_{i=1}^n\frac{{-i\choose n-i+1}{n\choose k}}{{i-1\choose k}-{n\choose k}}$$

$$= \sum_{i=0}^{n-1}\frac{{-1-i\choose n-i}{n\choose k}}{{i\choose k}-{n\choose k}}$$

$$= \sum_{i=0}^{n-1}\frac{\prod_{j=i+1}^n(-j)\prod_{j=0}^{k-1}(n-j)}{\prod_{j=1}^{n-i}j\left(\prod_{j=0}^{k-1}(i-j)-\prod_{j=0}^{k-1}(n-j)\right)}$$

ＮθＮηＩΣＥθ∕×Ｘ±¹ιΠ⁻θ⁺…⁰η…⁰⊕ι×Π…·¹⊕ι⁻Π⁻θ…⁰ηΠ⁻⁻θ⊕ι…⁰η

Attempt This Online! Link is to verbose version of code. Explanation: Implements a modification of Method 1:

$$\sum_{j=1}^n (-1)^{j-1} \binom n j \frac 1 {1 - \binom {n-j} k / \binom n k}$$

$$= \sum_{j=1}^n \frac {(-1)^{j-1} \binom n j \binom n k} {\binom n k - \binom {n-j} k}$$

$$= \sum_{j=1}^n \frac {(-1)^{j-1} \frac {n!} {j! (n-j)!} \frac {n!} {(n-k)!}} {\frac {n!} {(n-k)!} - \frac {(n-j)!} {(n-j-k)!}}$$

$$= \sum_{j=1}^n \frac {(-1)^{j-1} \prod_{i=n-j+1}^n i \prod_{i=n-k+1}^n i} {\prod_{i=1}^j i (\prod_{i=n-k+1}^n i - \prod_{i=n-j-k+1}^{n-j} i)}$$

$$= \sum_{j=1}^n \frac {(-1)^{j-1} \prod_{i=0}^{j-1} (n-i) \prod_{i=0}^{k-1} (n-i)} {\prod_{i=1}^j i (\prod_{i=0}^{k-1} (n-i) - \prod_{i=0}^{k-1} (n-j-i))}$$

$$= \sum_{j=0}^{n-1} \frac {(-1)^j \prod_{i=0}^j (n-i) \prod_{i=0}^{k-1} (n-i)} {\prod_{i=1}^{j+1} i (\prod_{i=0}^{k-1} (n-i) - \prod_{i=0}^{k-1} (n-j-i-1))}$$

It's then possible to factor the $$\ (-1)^j \$$ into the following product and then switch the operands of the subtraction in the denominator to correct the sign:

$$= \sum_{j=0}^{n-1} \frac {\prod_{i=0}^j (i-n) \prod_{i=0}^{k-1} (n-i)} {\prod_{i=1}^{j+1} i (\prod_{i=0}^{k-1} (n-j-i-1) - \prod_{i=0}^{k-1} (n-i))}$$

However this only reduces the code to 45 bytes:

ＮθＮηＩΣＥθ∕Π⁺⁻θ…⁰η⁻…⁰⊕ιθ×Π…·¹⊕ι⁻Π⁻⁻θ⊕ι…⁰ηΠ⁻θ…⁰η

Attempt This Online! Link is to verbose version of code.

Substituting $$\ i = n - j - 1 \$$ then results in the final formula used by the 42-byte solution.

Tried all three methods, with some minor changes to each for golfiness. In order:

# Wolfram Language (Mathematica), 52 bytes

b[i=-Range@#,++i+#].(b[-i,#2]/b@##-1)^-1&
b=Binomial

Try it online!

Shortest by a long shot. $$E(n,k)=\sum_{i=1}^n{-i\choose 1-i+n}\left(\frac{i-1\choose k}{n\choose k}-1\right)^{-1}.$$

### Wolfram Language (Mathematica), 82 bytes

Tr@Last@Inverse@Table[Boole[i==j]-b[#-i,#2-i+j]i~b~j/b@##,{i,#},{j,#}]&
b=Binomial

Try it online!

Reversed both axes of the matrix, summing the last row instead.

### Wolfram Language (Mathematica), 8583 82 bytes

{n,k}Check[1/Tr[m=b[#,i=Range@#]b[n-#,k-i]],0](n~b~k+m.#0/@(#-i))&@n
b=Binomial

Try it online!

Reversed the domain.

• Interesting, I really thought the 2nd method would benefit a language like this, given all its builtins and such for matrices, but I guess the first method is just that much simpler than the other two methods. Oct 24 at 23:13
• @AidenChow The built-ins are even nicer for sums of 2-vectorizeable-term products :P
– att
Oct 25 at 5:50

# Nekomata, 13 bytes

→r$ÇƆ/←ŗ0ɔ$ᵑ∆

Attempt This Online!

A port of @Command Master's 05AB1E answer. Unfortunately, Nekomata doesn't have 1/0=0.

→r$ÇƆ/←ŗ0ɔ$ᵑ∆
→r              [0,1,...,n]
$Ç [0Ck, 1Ck, ..., nCk] Ɔ/ [0Ck/nCk, 1Ck/nCk, ..., (n-1)Ck/nCk] ←ŗ decrement and reciprocal 0ɔ Append 0$ᵑ∆   Take the delta n times

# JavaScript (Node.js), 81 bytes

n=>g=(k,j=n,c=(a,b=k)=>!b||c(a,--b)*(a-b)/~b)=>j&&c(n,j)/(c(n-j)/c(n)-1)+g(k,j-1)

Try it online!

Method 1

-5 bytes from tsh

-1 from Arnauld

• n=>g=(k,j=0,c=(a,b=k)=>a--?c(a,b)+c(a,b-1):!b)=>j++<n&&c(n,j)/(1-c(n-j)/c(n))-g(k,j)
– tsh
Oct 24 at 5:57
• 83 bytes if you don't mind introducing a bitwise op. Oct 24 at 8:29

# R, 53 bytes

\(n,k,i=1:n-1)sum((^=choose)(-i-1,n-i)/(i^k/n^k-1))

Attempt This Online!

Last few test cases produce incorrect results due to numeric overflow.

# Pari/GP, 74 64 bytes

Saved 10 bytes thanks to @alephalpha

Use the formula

$$E(n,k)=\sum_{i=1}^n{-i\choose 1-i+n}\left(\frac{i-1\choose k}{n\choose k}-1\right)^{-1}.$$

Try it online!

b=binomial
E(n,k)=sum(i=1,n,b(-i,1-i+n)*(1/(b(i-1,k)/b(n,k)-1)))
• b(n,k)=binomial(n,k) -> b=binomial, *(1/(...)) -> /(...). Oct 24 at 11:46

# APL (Dyalog APL), 61 bytes

f←{+/(¯1*j+1)×(j!⍺)×÷1-(⍵!⍺-j←⍳⍺)÷⍵!⍺}

Attempt This Online!

not as efficient in bytes as possible, just a translation of formula 1.

# Maxima, 71 70 bytes

Saved 1 byte thanks to @att

In Maxima, binomial(-i,1-i+n) can be written in binomial(-i,-1-n) to save byte(s).

Use the formula

$$E(n,k)=\sum_{i=1}^n{-i\choose n+1-i}\left(\frac{i-1\choose k}{n\choose k}-1\right)^{-1}.$$

Try it online!

E(n,k):=sum(binomial(-i,-1-n)/(binomial(i-1,k)/binomial(n,k)-1),i,1,n)
• 1-i+n -> -1-n
– att
Oct 25 at 8:18

# Jelly, 12 bytes

Ḷc⁹÷c’İ;0I⁸¡

Try it online!

A simple port of @CommandMaster’s 05AB1E answer.