Suppose we select dice \(D=\{i_1,\dots,i_K\}\). The expected value \(\mathbb{E}_D\) for this selection is

\[ \mathbb{E}_D=\frac{1}{6^K}\sum_{j_1,\ldots,j_K}(A_{i_1,j_1}+\dots+A_{i_K,j_K})^2-(C_{i_1}+\dots+C_{i_k}) \]

where the sum is taken over all tuples \((j_1,\dots,j_k)\) such that \(1\leq j_i\leq 6\) for each \(i\). There are \(6^K\) such tuples and each is equiprobable to roll.

Let \(S_i=A_{i,1}+\dots+A_{i,6}\) and \(Q_i=A_{i,1}^2+\dots+A_{i,6}^2\).

Using the binomial theorem to expand the square and merging equal terms:

\[ \mathbb{E}_D=\frac{1}{6^K}\left(\sum_{i\in D}6^{K-1}A_{i,j}^2+\sum_{i,k\in D}2\cdot 6^{K-2}A_{i,j}A_{k,l}\right)-\sum_{i\in D}C_i \]

\[ \mathbb{E}_D=\frac{1}{6}\sum_{i\in D}A_{i,j}^2+\frac{1}{18}\sum_{i,k\in D}A_{i,j}A_{k,l}-\sum_{i \in D} C_i \]

\[ \mathbb{E}_D=\frac{1}{36}(\sum_{i\in D}S_i)^2+\frac{1}{6}\sum_{i\in D}Q_i-\frac{1}{36}\sum_{i\in D}S_i^2-\sum_{i\in D} C_i \]

Let \(B_i=6Q_i-S_i^2-36C_i\). This simplifies to

\[ \mathbb{E}_D=\frac{1}{36}(\sum_{i \in D}S_i)^2+\frac{1}{36}\sum_{i\in D}B_i \]

Introduce a variable \(X\) and consider the function

\[ \mathbb{E}_D(X)=\frac{1}{36}(\sum_{i \in D}S_i)X+\frac{1}{36}\sum_{i \in D}B_i=\frac{1}{36}\sum_i(S_iX+B_i) \]

If the value of \(X\) is fixed, then to maximize \(\mathbb{E}_D(X)\) we greedily take the \(K\) dice with largest \(S_iX+B_i\). Call this set \(D_X\). If the global optimum dice selection is \(G\) and \(S_G=\sum_{i\in G}S_i\), then \(G\) is also an optimum of \(\mathbb{E}_D(S_G)\), so \(G=D_{S_G}\).

We don’t know \(S_G\), so let’s just visit every possible value: let \(L=\min_{D}\sum_{i\in D}S_i\) and \(R=\max_{D}\sum_{i\in D}S_i\). Run \(X\) through \([L,R]\), and for each \(X\) find the set \(D_X\) of the best dice that can be used for this \(X\). We will visit \(G\) when \(X=S_G\). This takes \(O(N(R-L))\) time.

To speed this up we must skip consecutive values of \(X\) that yield the same dice selection. Consider an ordering of all \(N\) dice by the value of \(S_iX+B_i\). If it is non-increasing the optimal selection is the \(K\) first dice on it.

Suppose that all \(S_i\) are distinct. Then the orderings for \(X=-\infty\) and \(X=+\infty\) are reversed. When \(X\) sweeps continuously from \(-\infty\) to \(+\infty\), at certain moments dice in this ordering swap relative positions. Two dice \(i\) and \(j\) such that \(S_i\neq S_j\) swap positions when \(S_iX+B_i=S_jX+B_j\), or equivalently when \(X=\tfrac{B_i-B_j}{S_j-S_i}\).

If at this moment \(i\) and \(j\) are not adjacent in the ordering, then all the dice between them will have the same value of \(S_{*}X+B_{*}\) at the moment of the swap, and thus the whole range gets reversed. Hence to simulate this process it is sufficient to consider only swaps of dice in adjacent positions.

If \(S_i=S_j\) the relative ordering of \(i\) and \(j\) does not change throughout the sweep. Each dice pair \((i,j)\) therefore swap positions at most once, so there are at most \(\binom{N}{2}=O(N^2)\) swaps in this sweep.

The sweep can be simulated directly in \(O(N^3)\) by checking every adjacent pair of dice at each step, and performing the swap that occurs earliest in time. This can be further optimized to \(O(N^2\log N)\) online by tracking the swap events in a set or heap, or offline by computing and sorting all the swaps upfront. The expected values \(\mathbb{E}_D\) can be tracked during the process in constant time.