Introduction

As a rhythm game player, I often wonder what my max combo will be in my next play. This is a rather unpredictable outcome, and what I can do is to try to conclude a probability distribution of my max combo.

For those who are not familiar with rhythm games and also to make the question clearer, I state the problem in a more mathematical setting.

Consider a random bit string of length n∈Nn\in\mathbb N, where each bit is independent and has probability Y∈[0,1]Y\in[0,1] of being 11. Let Pn,k(Y)P_{n,k}(Y) be the probability that the length of the longest all-11 substring of the bit string is k∈Nk\in\mathbb N (where obviously Pn,k(Y)P_{n,k}(Y) is nonzero only when k≤nk\le n). What is the expression of Pn,k(Y)P_{n,k}(Y)?

A more interesting problem to consider is what the probability distribution tends to be when n→∞n\to\infty. Define the random variable κ≔k/n\kappa\coloneqq k/n where kk is the length of the longest all-11 substring. Define a parameter y≔Yny\coloneqq Y^n (this parameter is held constant while n→∞n\to\infty). Define the probability distribution function of κ\kappa as f(y,κ)≔lim⁡n→∞(n+1)Pn,κn ⁣(y1n).f(y,\kappa)\coloneqq\lim_{n\to\infty}\left(n+1\right)P_{n,\kappa n}\!\left(y^{\frac1n}\right). (1)(1) What is the expression of f(y,κ)f(y,\kappa)?

Notation

Notation for integer range: a…ba\ldots b denotes the integer range defined by the ends aa (inclusive) and bb (exclusive), or in other words {a,a+1,…,b−1}\left\{a,a+1,\ldots,b-1\right\}. It is defined to be empty if a≥ba\ge b. The operator …\ldots has a lower precedence than ++ and −- but a higher precedence than ∈\in.

The notation a ..ba\,..b denotes the inclusive integer range {a,a+1,…,b}\left\{a,a+1,\ldots,b\right\}. It is defined to be empty if a>ba>b.

The case for finite nn

A natural approach to find Pn,kP_{n,k} is to try to find a recurrence relation of Pn,kP_{n,k} for different nn and kk, and then use a dynamic programming (DP) algorithm to compute Pn,kP_{n,k} for any given nn and kk.

The first DP approach

For a rhythm game player, the most straightforward way of finding kk for a given bit string is to track the current combo, and update the max combo when the current combo is greater than the previous max combo.

To give the current combo a formal definition, denote each bit in the bit string as bib_i, where i∈0…ni\in0\ldots n. Define the current combo rir_i as the length of the longest all-11 substring of the bit string ending before (exclusive) ii (so ri=0r_i=0 if bi−1=0b_{i-1}=0, which is callled a combo break): ri≔max⁡{r∈0 ..i | ∀j∈i−r…i:bj=1},r_i\coloneqq\max\left\{r\in0\,..i\,\middle|\,\forall j\in i-r\ldots i:b_j=1\right\}, where i∈0 ..ni\in0\,..n.

Now, use three numbers (n,k,r)(n,k,r) to define a DP state. Denote Pn,k,rP_{n,k,r} to be the probability that the max combo is kk and the final combo (rnr_n) is rr. Then, consider a transition from state (n,k,r)(n,k,r) to state (n+1,k′,r′)(n+1,k',r') by adding a new bit bnb_n to the bit string. There are two cases:

  • If bn=0b_n=0 (has 1−Y1-Y probability), then this means a combo break, so we have r′=0r'=0 and k′=kk'=k.
  • If bn=1b_n=1 (has YY probability), then the combo continues, so we have r′=r+1r'=r+1. The max combo needs to be updated if needed, so we have k′=max⁡(k,r′)k'=\max(k,r').

However, in actual implementation of the DP algorithm, we need to reverse this transition by considering what state can lead to the current state (n,k,r)(n,k,r) (to use the bottom-up approach).

First, obviously in any possible case r∈0 ..kr\in0\,..k (currently we only consider the cases where n>k>0n>k>0). Divide all those cases into three groups:

  1. If r=0r=0, this is means a combo break, so the last bit is 00, and the previous state can have any possible final combo r′r'. Therefore, it can be transitioned from any (n−1,k,r′)(n-1,k,r') where r′∈0 ..kr'\in0\,..k. For each possible previous state, the probability of the transition to this new state is 1−Y1-Y.
  2. If r∈1 ..k−1r\in1\,..k-1, this means the last bit is 11, the previous final combo is r−1r-1, and the previous max combo is already kk. Therefore, the previous state is (n−1,k,r−1)(n-1,k,r-1), and the probability of the transition is YY.
  3. If r=kr=k, this means the max combo may (or may not) have been updated. In either case, the previous final combo is r−1=k−1r-1=k-1.
  • If the max combo is updated, the previous max combo must be k−1k-1 because it must not be less than the previous final combo k−1k-1 and must be less than the new max combo kk. Therefore, the previous state is (n−1,k−1,k−1)(n-1,k-1,k-1), and the probability of the transition is YY.
  • If the max combo is not updated, the previous max combo is the same as the new one, which is kk. Therefore, the previous state is (n−1,k,k−1)(n-1,k,k-1), and the probability of the transition is YY.

Therefore, we can write a recurrence relation that is valid when n>k>0n>k>0: Pn,k,r={(1−Y)∑r′=0kPn−1,k,r′,r=0YPn−1,k,r−1,r∈1 ..k−1Y(Pn−1,k−1,k−1+Pn−1,k,k−1),r=k.P_{n,k,r}=\begin{cases} \left(1-Y\right)\sum_{r'=0}^kP_{n-1,k,r'},&r=0\\ YP_{n-1,k,r-1},&r\in1\,..k-1\\ Y\left(P_{n-1,k-1,k-1}+P_{n-1,k,k-1}\right),&r=k. \end{cases}

However, there are also other cases (mostly edge cases) because we assumed n>k>0n>k>0. Actually, in the meaningfulness condition n≥k≥r≥0n\ge k\ge r\ge0 (necessary condition for Pn,k,rP_{n,k,r} to be nonzero), there are three inequality that can be altered between a less-than sign or an equal sign, so there are totally 23=82^3=8 cases. Considering all those cases (omitted in this article because of the triviality), we can write a recurrence relation that is valid for all n,k,rn,k,r, covering all the edge cases: Pn,k,r={1,n=k=r=0,YPn−1,n−1,n−1,n=k=r>0,0,n=k>r>0,0,n=k>r=0,Y(Pn−1,k−1,k−1+Pn−1,k,k−1),n>k=r>0,YPn−1,k,r−1,n>k>r>0,(1−Y)∑r′=0kPn−1,k,r′,n>k>r=0,(1−Y)Pn−1,0,0,n>k=r=0.P_{n,k,r}=\begin{cases} 1,& n=k=r=0,\\ YP_{n-1,n-1,n-1},& n=k=r>0,\\ 0,& n=k>r>0,\\ 0,& n=k>r=0,\\ Y\left(P_{n-1,k-1,k-1}+P_{n-1,k,k-1}\right),& n>k=r>0,\\ YP_{n-1,k,r-1},& n>k>r>0,\\ \left(1-Y\right)\sum_{r'=0}^kP_{n-1,k,r'},& n>k>r=0,\\ \left(1-Y\right)P_{n-1,0,0},& n>k=r=0. \end{cases} (2)(2)

Note that the probabilities related to note count nn only depend on those related to note count n−1n-1 and that the probabilities related to max combo kk and final combo rr only depend on those related to either less max combo than kk or less final combo than rr (except for the case n>k>r=0n>k>r=0, which can be specially treated before the current iteration of kk actually starts), so for the bottom-up DP we can reduce the spatial complexity from O ⁣(n3)O\!\left(n^3\right) to O ⁣(n2)O\!\left(n^2\right) by reducing the 3-dimensional DP to a 2-dimensional one. What needs to be taken care of is that the DP table needs to be updated from larger kk and rr to smaller kk and rr instead of the other way so that the numbers in the last iteration in nn are left untouched while we need to use them in the current iteration.

After the final iteration in nn finishes, we need to sum over the index rr to get the final answer: Pn,k=∑r=0kPn,k,r.P_{n,k}=\sum_{r=0}^kP_{n,k,r}.

Writing the code for the DP algorithm is then straightforward. Here is an implementation in Ruby. In the code, dp[k][r] means Pn,k,rP_{n,k,r} in the nnth iteration.

## Returns an array of size m+1,
## with the k-th element being the probability P_{m,k}.
def combo m
	(1..m).each_with_object [[1]] do |n, dp|
		dp[n] = [0]*n + [Y * dp[n-1][n-1]] # n = k > 0
		(n-1).downto 1 do |k| # n > k > 0
			dpk0 = (1-Y) * dp[k].sum
			dp[k][k] = Y * (dp[k-1][k-1] + dp[k][k-1])        # n > k = r > 0
			(k-1).downto(1) { |r| dp[k][r] = Y * dp[k][r-1] } # n > k > r > 0
			dp[k][0] = dpk0                                   # n > k > r = 0
		end
		dp[0][0] *= 1-Y # n > k = r = 0
	end.map &:sum
end

Because of the three nested loops, the time complexity of the DP algorithm is O ⁣(n3)O\!\left(n^3\right).

The second DP approach

Here is an alternative way to use DP to solve the problem. Instead of building a DP table with the k,rk,r indices, we can build a DP table with the n,kn,k indices.

First, we need to rewrite the recurrence relation of Pn,kP_{n,k} instead of that of Pn,k,rP_{n,k,r}. We then need to try to express Pn,k,rP_{n,k,r} in terms of Pn,kP_{n,k} terms. The easiest part is the case where n≥k=r=0n\ge k=r=0. By recursively applying Equation 2 to Pn,0,0P_{n,0,0}, we have Pn,0,0=(1−Y)Pn−1,0,0=(1−Y)2Pn−2,0,0=⋯=(1−Y)nP0,0,0.\begin{align*} P_{n,0,0}&=\left(1-Y\right)P_{n-1,0,0}\\ &=\left(1-Y\right)^2P_{n-2,0,0}\\ &=\cdots\\ &=\left(1-Y\right)^nP_{0,0,0}. \end{align*} Because P0,0,0=1P_{0,0,0}=1, we have Pn,0,0=(1−Y)n.P_{n,0,0}=\left(1-Y\right)^n. (3)(3)


For n>k>r>0n>k>r>0, we can recursively apply Equation 2 to get Pn,k,r=YPn−1,k,r−1=Y2Pn−2,k,r−2=⋯\begin{align*} P_{n,k,r}&=YP_{n-1,k,r-1}\\ &=Y^2P_{n-2,k,r-2}\\ &=\cdots \end{align*} This will finally either decend the note count to kk or decend the final combo to 00, determined by which comes first.

  • If n−r≤kn-r\le k, we will decend to the term Pk,k,r−(n−k)P_{k,k,r-(n-k)}, which must be zero according to the case n=k>r=0n=k>r=0 and the case n=k>r>0n=k>r>0 in Equation 2, so Pn,k,r=0P_{n,k,r}=0.
  • If n−r>kn-r>k, then we will decend to the term Pn−r,k,0P_{n-r,k,0}, which is equal to (1−Y)Pn−r−1,k\left(1-Y\right)P_{n-r-1,k} according to the case n>k>r=0n>k>r=0 in Equation 2.

Therefore, for n>k>r>0n>k>r>0, we have Pn,k,r={0,n−r≤k,Yr(1−Y)Pn−r−1,k,n−r>k.P_{n,k,r}=\begin{cases} 0,&n-r\le k,\\ Y^r\left(1-Y\right)P_{n-r-1,k},&n-r>k. \end{cases} (4)(4)


For the case n>k=r>0n>k=r>0, we can also recursively apply Equation 2 to get Pn,k,k=Y(Pn−1,k−1,k−1+Pn−1,k,k−1)=Y(Y(Pn−2,k−2,k−2+Pn−2,k−1,k−2)+Pn−1,k,k−1)=⋯=Y(Y(⋯Y(Pn−k,0,0+Pn−k,1,0)+⋯ )+Pn−1,k,k−1)=YkPn−k,0,0+∑j=1kYjPn−j,k−j+1,k−j.\begin{align*} P_{n,k,k}&=Y\left(P_{n-1,k-1,k-1}+P_{n-1,k,k-1}\right)\\ &=Y\left(Y\left(P_{n-2,k-2,k-2}+P_{n-2,k-1,k-2}\right)+P_{n-1,k,k-1}\right)\\ &=\cdots\\ &=Y\left(Y\left(\cdots Y\left(P_{n-k,0,0}+P_{n-k,1,0}\right)+\cdots\right)+P_{n-1,k,k-1}\right)\\ &=Y^kP_{n-k,0,0}+\sum_{j=1}^kY^jP_{n-j,k-j+1,k-j}. \end{align*} We can then substitute Equation 3 and 4 into the above equation. The substitution of Equation 3 can be done without a problem, but the substitution of Equation 4 requires some care because of the different cases.

  • If n−k>kn-k>k, then only the case n−r>kn-r>k in Equation 4 will be involved in the summation.
  • If n−k≤kn-k\le k, then both cases in Equation 4 will be involved in the summation. To be specific, for j∈1 .. 2k−n+1j\in1\,..\,2k-n+1, we need the case n−r≤kn-r\le k in Equation 4 (where the summed terms are just zero and can be omitted); for other terms in the summation, we need the other case.

Considering both cases, we may realize that we can just modify the range of the summation to j∈max⁡(1,2k−n+1) ..kj\in\max(1,2k-n+1)\,..k and adopt the case n−r>kn-r>k in Equation 4 for all terms in the summation. Therefore, we have Pn,k,k=Yk(1−Y)n−k+∑j=max⁡(1,2k−n+1)kYjYk−j(1−Y)Pn−j−(k−j)−1,k−j+1=Yk(1−Y)n−k+Yk(1−Y)∑k′=1min⁡(k,n−k−1)Pn−k−1,k′,\begin{align*} P_{n,k,k}&=Y^k\left(1-Y\right)^{n-k}+\sum_{j=\max(1,2k-n+1)}^{k}Y^jY^{k-j}\left(1-Y\right)P_{n-j-(k-j)-1,k-j+1}\\ &=Y^k\left(1-Y\right)^{n-k}+Y^k\left(1-Y\right)\sum_{k'=1}^{\min(k,n-k-1)}P_{n-k-1,k'}, \end{align*} where in the last line we changed the summation index to k′≔k−j+1k'\coloneqq k-j+1 to simplify it. Because Pn−k−1,0=Pn−k−1,0,0=(1−Y)n−k−1P_{n-k-1,0}=P_{n-k-1,0,0}=\left(1-Y\right)^{n-k-1} according to Equation 3, we can combine the two terms into one summation to get the final result for n>k=r>0n>k=r>0: Pn,k,k=Yk(1−Y)∑k′=0min⁡(k,n−k−1)Pn−k−1,k′.P_{n,k,k}=Y^k\left(1-Y\right)\sum_{k'=0}^{\min(k,n-k-1)}P_{n-k-1,k'}. (5)(5) Noticing the obvious fact that ∑k=0nPn,k=1\sum_{k=0}^nP_{n,k}=1, the above equation can be simplified, when k≥n−k−1k\ge n-k-1, to Pn,k,k=Yk(1−Y).P_{n,k,k}=Y^k\left(1-Y\right). (6)(6) This simplification is not specially useful, but it can be used to simplify the calculation in the program.


Then, for n>k>0n>k>0, express Pn,kP_{n,k} in terms of Pn,k,rP_{n,k,r} by summing over rr, and substitute previous results: Pn,k=∑r=0kPn,k,r=Pn,k,0+Pn,k,k+∑r=1k−1Pn,k,r=(1−Y)Pn−1,k+Yk(1−Y)∑k′=0min⁡(k,n−k−1)Pn−k−1,k′= +∑r=1min⁡(k−1,n−k−1)Yr(1−Y)Pn−r−1,k=(1−Y)(Yk∑k′=0min⁡(k,n−k−1)Pn−k−1,k′+∑r=0min⁡(k−1,n−k−1)YrPn−r−1,k)\begin{align*} P_{n,k}&=\sum_{r=0}^kP_{n,k,r}\\ &=P_{n,k,0}+P_{n,k,k}+\sum_{r=1}^{k-1}P_{n,k,r}\\ &=\left(1-Y\right)P_{n-1,k}+Y^k\left(1-Y\right)\sum_{k'=0}^{\min(k,n-k-1)}P_{n-k-1,k'}\\ &\phantom{=~}{}+\sum_{r=1}^{\min(k-1,n-k-1)}Y^r\left(1-Y\right)P_{n-r-1,k}\\ &=\left(1-Y\right)\left( Y^k\sum_{k'=0}^{\min(k,n-k-1)}P_{n-k-1,k'} +\sum_{r=0}^{\min(k-1,n-k-1)}Y^rP_{n-r-1,k} \right) \end{align*} where in the last term rr is summed to min⁡(k−1,n−k−1)\min(k-1,n-k-1) instead of k−1k-1 because of the different cases in Equation 4.

Finally, consider the edge cases where n=k≥0n=k\ge0 and n≥k=0n\ge k=0 (trivial), we have the complete resursive relation for Pn,kP_{n,k}: Pn,k={Yn,n=k≥0,(1−Y)n,n≥k=0,(1−Y) (Yk∑k′=0min⁡(k,n−k−1)Pn−k−1,k′+∑r=0min⁡(k−1,n−k−1)YrPn−r−1,k),n>k>0.P_{n,k}=\begin{cases} Y^n,&n=k\ge0,\\ \left(1-Y\right)^n,&n\ge k=0,\\ \displaystyle{\begin{split} \left(1-Y\right)&\,\left( Y^k\sum_{k'=0}^{\min(k,n-k-1)}P_{n-k-1,k'} \right.\\&\left. +\sum_{r=0}^{\min(k-1,n-k-1)}Y^rP_{n-r-1,k} \right), \end{split}}&n>k>0. \end{cases} (7)(7)

Then, we can write the program to calculate Pn,kP_{n,k}:

## Returns an array of size m+1,
## with the k-th element being the probability P_{m,k}.
def combo m
	(1..m).each_with_object [[1]] do |n, dp|
		dp[n] = (1..n-1).each_with_object [(1-Y)**n] do |k, dpn|
			dpn[k] = (1-Y) * (Y**k * (0..[k, n-k-1].min).sum { dp[n-k-1][_1] } + (0..[k-1, n-k-1].min).sum { Y**_1 * dp[n-_1-1][k] })
		end
		dp[n][n] = Y**n
	end.last
end

This algorithm has the same (asymptotic) space and time complexity as the previous one.

Polynomial coefficients

We have wrote programmes to calculate probabilities Pn,k(Y)P_{n,k}(Y) based on given YY, which we assumed to be a float number. However, float numbers have limited precision, and the calculation may be inaccurate. Actually, the calculation can be done symbolically.

The probability Pn,kP_{n,k} is a polynomial of degree (at most) nn in YY, and the coefficients of the polynomial are integers. This can be easily proven by using mathematical induction and utilizing Equation 7. Therefore, we can calculate the coefficients of the polynomial Pn,k(Y)P_{n,k}(Y) instead of calculate the value directly so that we get a symbolic but accurate result.

Both the two DP algorithms above can be modified to calculate the coefficients of the polynomial. Actually, we can define Y to be a polynomial object that can do arithmetic operations with other polynomials or numbers, and then the programmes can run without any modification. Here, I will modify the second DP algorithm to calculate the coefficients of the polynomial.

We can also utilize Equation 6 to simplify the calculation. Considering the edge cases involved in min⁡(k,n−k−1)\min(k,n-k-1) and min⁡(k−1,n−k−1)\min(k-1,n-k-1), there are three cases we need to consider:

  1. Case k>n−k−1k>n-k-1: Equation 6 can be applied, and rr is summed to n−k−1n-k-1.
  2. Case k=n−k−1k=n-k-1 (can only happen when nn is odd): Equation 6 can be applied, and rr is summed to k−1k-1.
  3. Case k<n−k−1k<n-k-1: Equation 6 cannot be applied, and rr is summed to k−1k-1.

Then, use arrays to store the coefficients of the polynomial Pn,k(Y)P_{n,k}(Y), and we can write the program to calculate the coefficients:

## Returns a nested array of size m+1 times m+1,
## with the j-th element of the k-th element being the coefficient of Y^j in P_{m,k}(Y).
def combo_pc m
	(1..m).each_with_object [[[1]]] do |n, dp|
		dp[n] = Array.new(n+1) { Array.new n+1, 0 }

		# dp[n][0] = (1-Y)**n
		0.upto(n-1) { dp[n][0][_1] = dp[n-1][0][_1] } # will be multiplied by 1-Y later

		1.upto n/2-1 do |k|
			# dp[n][k] = (1-Y) * (Y**k * (0..k).sum { |j| dp[n-k-1][j] } + (0..k-1).sum { |r| Y**r * dp[n-r-1][k] })
			0.upto(k) { |j| 0.upto(n-k-1) { dp[n][k][_1+k] += dp[n-k-1][j][_1] } }
			0.upto(k-1) { |r| 0.upto(n-r-1) { dp[n][k][_1+r] += dp[n-r-1][k][_1] } }
		end

		if n % 2 == 1
			k = n/2
			# dp[n][k] = (1-Y) * (Y**k + (0..k-1).sum { |r| Y**r * dp[n-r-1][k] })
			dp[n][k][k] = 1
			0.upto(k-1) { |r| 0.upto(n-r-1) { dp[n][k][_1+r] += dp[n-r-1][k][_1] } }
		end

		((n+1)/2).upto n-1 do |k|
			# dp[n][k] = (1-Y) * (Y**k + (0..n-k-1).sum { |r| Y**r * dp[n-r-1][k] })
			dp[n][k][k] = 1
			0.upto(n-k-1) { |r| 0.upto(n-r-1) { dp[n][k][_1+r] += dp[n-r-1][k][_1] } }
		end

		0.upto(n-1) { |k| n.downto(1) { dp[n][k][_1] -= dp[n][k][_1-1] } } # multiply by 1-Y

		# dp[n][n] = Y**n
		dp[n][n][n] = 1
	end.last
end

Here I list first few polynomials Pn,k(Y)P_{n,k}(Y) calculated by the above program: k=0123⋯n=0111−YY21−2Y+Y22Y−2Y2Y231−3Y+3Y2−Y33Y−5Y2+2Y32Y2−2Y3Y3⋮\begin{array}{r|llllc} & k=0 & 1 & 2 & 3 & \cdots\\ \hline n=0 & 1\\ 1 & 1-Y & Y\\ 2 & 1-2Y+Y^2 & 2Y-2Y^2 & Y^2\\ 3 & 1-3Y+3Y^2-Y^3 & 3Y-5Y^2+2Y^3 & 2Y^2-2Y^3 & Y^3\\ \vdots \end{array}

When evaluating the polynomials for large nn, the result is inaccurate for YY that is not close to 00 because of the limited precision of floating numbers. If YY is closer to 11, we can first find the coefficients of Pn,k(1−X)P_{n,k}(1-X) and then substitute X≔1−YX\coloneqq1-Y.

Plots of the probability distributions

Here are some plots of the probability distribution of max combo kk when n=50n=50:

Probability distribution of  when  for different

The plots are intuitive as they show that one has higher probability to get a higher max combo when they have a higher success rate.

There is a suspicious jump in Pn,k(Y)P_{n,k}(Y) near k=n/2k=n/2 when YY is close to 11. We can look at it closer:

Probability distribution of  when  for different

In the zoomed-in plot, we can also see a jump in first derivative (w.r.t. kk) of Pn,k(Y)P_{n,k}(Y) near k=n/3k=n/3. Actually, the jumps can be modeled in later sections when we talk about the case when n→∞n\to\infty.

The case when n→∞n\to\infty

A natural approach is to try substituting Equation 7 into Equation 1 to get a function w.r.t. the unknown function f(y,κ)f(y,\kappa). First, we can easily write the case when y=0y=0 because it means zero success rate, and the only possible max combo is zero: f(y=0,κ)=δ(κ).f(y=0,\kappa)=\delta(\kappa). (8)(8) Similarly, we can easily write the case when y=1y=1: f(y=1,κ)=δ(κ−1).f(y=1,\kappa)=\delta(\kappa-1). (9)(9)

From now on, we only consider the case when 0<y<10<y<1. First, for the case κ=0\kappa=0, according to Equation 7, f(y,κ=0)=lim⁡n→∞(n+1)(1−y1n)n={0,0<y≤1,∞,y=0.\begin{align*} f(y,\kappa=0)&=\lim_{n\to\infty}\left(n+1\right)\left(1-y^{\frac1n}\right)^n\\ &=\begin{cases}0,&0<y\le1,\\\infty,&y=0.\end{cases} \end{align*} The ∞\infty means that there is a Dirac δ\delta function (shown in Equation 8).

Then, for the case κ=1\kappa=1, according to Equation 7, f(y,κ=1)=lim⁡n→∞(n+1)y=∞.f(y,\kappa=1)=\lim_{n\to\infty}\left(n+1\right)y=\infty. The ∞\infty means that there is a Dirac δ\delta function. Actually, it is easy to see that there must be a yδ(κ−1)y\delta(\kappa-1) term in the expression of f(y,κ)f(y,\kappa) because the probability of getting a max combo (κ=1\kappa=1) is yy.

Define h(y,κ)≔f(y,κ)−yδ(κ−1),h(y,\kappa)\coloneqq f(y,\kappa)-y\delta(\kappa-1), (10)(10) and then we can get rid of the infinity here.

From now on, we only consider the case when 0<y<10<y<1 and 0<κ<10<\kappa<1. According to Equation 7, = f ⁣(y∈(0,1),κ∈(0,1))=lim⁡n→∞(n+1)(1−y1n)(yκ∑k′=0min⁡(κn,n−κn−1)Pn−κn−1,k′ ⁣(y1n)=lim⁡n→∞(n+1)(1−y1n)+∑r=0min⁡(κn−1,n−κn−1)yrnPn−r−1,κn ⁣(y1n))=lim⁡n→∞n(1−y1n)⋅lim⁡n→∞(yκ∑t=0,Δt=1nmin⁡(κ,1−κ)P(1−κ)n,tn ⁣(y1n)+∑t=0,Δt=1nmin⁡(κ,1−κ)ytP(1−t)n,κn ⁣(y1n))=−ln⁡ylim⁡n→∞(yκ∑t=0,Δt=1nmin⁡(κ,1−κ)1(1−κ)nf ⁣(y1−κ,t1−κ)=−ln⁡ylim⁡n→∞+∑t=0,Δt=1nmin⁡(κ,1−κ)yt1(1−t)nf ⁣(y1−t,κ1−t))=−ln⁡ylim⁡n→∞∑t=0,Δt=1nmin⁡(κ,1−κ)(yκ1−κf ⁣(y1−κ,t1−κ)+yt1−tf ⁣(y1−t,κ1−t))Δt=−ln⁡y∫t=0min⁡(κ,1−κ)(yκ1−κf ⁣(y1−κ,t1−κ)+yt1−tf ⁣(y1−t,κ1−t))dt.\begin{align*} &\phantom{=~}f\!\left(y\in\left(0,1\right),\kappa\in\left(0,1\right)\right)\\ &=\lim_{n\to\infty}\left(n+1\right)\left(1-y^{\frac1n}\right)\left( y^\kappa\sum_{k'=0}^{\min(\kappa n,n-\kappa n-1)}P_{n-\kappa n-1,k'}\!\left(y^{\frac1n}\right) \right.\\&\phantom{=\lim_{n\to\infty}\left(n+1\right)\left(1-y^{\frac1n}\right)}\left. +\sum_{r=0}^{\min(\kappa n-1,n-\kappa n-1)}y^{\frac rn}P_{n-r-1,\kappa n}\!\left(y^{\frac1n}\right) \right)\\ &=\lim_{n\to\infty}n\left(1-y^{\frac1n}\right)\cdot\lim_{n\to\infty}\left( y^\kappa\sum_{t=0,\Delta t=\frac1n}^{\min(\kappa,1-\kappa)}P_{(1-\kappa)n,tn}\!\left(y^{\frac1n}\right) +\sum_{t=0,\Delta t=\frac1n}^{\min(\kappa,1-\kappa)}y^tP_{(1-t)n,\kappa n}\!\left(y^{\frac1n}\right) \right)\\ &=-\ln y\lim_{n\to\infty}\left( y^\kappa\sum_{t=0,\Delta t=\frac1n}^{\min(\kappa,1-\kappa)} \frac1{\left(1-\kappa\right)n}f\!\left(y^{1-\kappa},\frac t{1-\kappa}\right) \right.\\&\phantom{=-\ln y\lim_{n\to\infty}}\left. +\sum_{t=0,\Delta t=\frac1n}^{\min(\kappa,1-\kappa)}y^t \frac1{\left(1-t\right)n}f\!\left(y^{1-t},\frac{\kappa}{1-t}\right) \right)\\ &=-\ln y\lim_{n\to\infty}\sum_{t=0,\Delta t=\frac 1n}^{\min(\kappa,1-\kappa)}\left( \frac{y^\kappa}{1-\kappa}f\!\left(y^{1-\kappa},\frac t{1-\kappa}\right) +\frac{y^t}{1-t}f\!\left(y^{1-t},\frac\kappa{1-t}\right) \right)\Delta t\\ &=-\ln y\int_{t=0}^{\min(\kappa,1-\kappa)}\left( \frac{y^\kappa}{1-\kappa}f\!\left(y^{1-\kappa},\frac t{1-\kappa}\right) +\frac{y^t}{1-t}f\!\left(y^{1-t},\frac\kappa{1-t}\right) \right)\mathrm dt. \end{align*} Add back the delta function at κ=1\kappa=1, and we have the integral equation f ⁣(y∈(0,1),κ)=−ln⁡y∫t=0min⁡(κ,1−κ) (yκ1−κf ⁣(y1−κ,t1−κ)+yt1−tf ⁣(y1−t,κ1−t))dt+yδ(κ−1).\begin{split} f\!\left(y\in\left(0,1\right),\kappa\right)=-\ln y\int_{t=0}^{\min(\kappa,1-\kappa)}&\,\left( \frac{y^\kappa}{1-\kappa}f\!\left(y^{1-\kappa},\frac t{1-\kappa}\right) \right.\\&\left. +\frac{y^t}{1-t}f\!\left(y^{1-t},\frac\kappa{1-t}\right) \right)\mathrm dt+y\delta(\kappa-1). \end{split}

There are two terms in the integral. Substitute u≔t1−κu\coloneqq\frac t{1-\kappa} in the first term, and we have ∫t=0min⁡(κ,1−κ)yκ1−κf ⁣(y1−κ,t1−κ)dt=∫0min⁡(κ1−κ,1)yκf ⁣(y1−κ,u)du=∫0min⁡(κ1−κ,1)(yκh ⁣(y1−κ,u)+yδ(u−1))du.\begin{align*} \int_{t=0}^{\min\left(\kappa,1-\kappa\right)}\frac{y^\kappa}{1-\kappa}f\!\left(y^{1-\kappa},\frac t{1-\kappa}\right)\mathrm dt &=\int_0^{\min(\frac\kappa{1-\kappa},1)}y^\kappa f\!\left(y^{1-\kappa},u\right)\mathrm du\\ &=\int_0^{\min(\frac\kappa{1-\kappa},1)}\left(y^\kappa h\!\left(y^{1-\kappa},u\right)+y\delta(u-1)\right)\mathrm du. \end{align*} Substitute v≔κ1−tv\coloneqq\frac\kappa{1-t} in the second term, and we have ∫t=0min⁡(κ,1−κ)yt1−tf ⁣(y1−t,κ1−t)dt=∫κmin⁡(κ1−κ,1)y1−κvf ⁣(yκv,v)dvv=∫κmin⁡(κ1−κ,1)(y1−κvh ⁣(yκv,v)+yδ(v−1))dvv.\begin{align*} \int_{t=0}^{\min(\kappa,1-\kappa)}\frac{y^t}{1-t}f\!\left(y^{1-t},\frac\kappa{1-t}\right)\mathrm dt &=\int_\kappa^{\min\left(\frac\kappa{1-\kappa},1\right)}y^{1-\frac\kappa v}f\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v\\ &=\int_\kappa^{\min\left(\frac\kappa{1-\kappa},1\right)}\left(y^{1-\frac\kappa v}h\!\left(y^{\frac\kappa v},v\right)+y\delta(v-1)\right)\frac{\mathrm dv}v. \end{align*}

Further, let (we only consider y∈(0,1)y\in\left(0,1\right) from now on) g(y,κ)≔h(y,κ)y,g(y,\kappa)\coloneqq\frac{h(y,\kappa)}{y}, (11)(11) then the integral equation becomes g(y,κ)=−ln⁡y (∫0min⁡(κ1−κ,1)(g ⁣(y1−κ,u)+δ(u−1))du+∫κmin⁡(κ1−κ,1)(g ⁣(yκv,v)+δ(v−1))dvv).\begin{split} g(y,\kappa)=-\ln y&\,\left( \int_0^{\min\left(\frac\kappa{1-\kappa},1\right)}\left(g\!\left(y^{1-\kappa},u\right)+\delta(u-1)\right)\mathrm du \right.\\&\left. +\int_\kappa^{\min\left(\frac\kappa{1-\kappa},1\right)}\left(g\!\left(y^{\frac\kappa v},v\right)+\delta(v-1)\right)\frac{\mathrm dv}v \right). \end{split} (12)(12)

There is another integral equation for gg. Because ∫01f ⁣(y,κ)dκ=1\int_0^1f\!\left(y,\kappa\right)\mathrm d\kappa=1, we have ∫01g ⁣(y,κ)dκ=1y−1.\int_0^1g\!\left(y,\kappa\right)\mathrm d\kappa=\frac1y-1. (13)(13)

Equation 12 and 13 are the equations that we are going to utilize to get the expression for g(y,κ)g(y,\kappa).

The case κ∈(12,1)\kappa\in\left(\frac12,1\right)

In this case, we have min⁡ ⁣(κ1−κ,1)=1,\min\!\left(\frac\kappa{1-\kappa},1\right)=1, so the Dirac delta functions in Equation 12 should be considered. In this case, it simplifies to g1(y,κ)≔g ⁣(y,κ∈(12,1))=−ln⁡y(yκ−1+∫κ1g ⁣(yκv,v)dvv+1),g_1(y,\kappa)\coloneqq g\!\left(y,\kappa\in\left(\frac12,1\right)\right)= -\ln y\left(y^{\kappa-1}+\int_\kappa^1g\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v+1\right), (14)(14) where Equation 13 is utilized when finding the first term.

We can try to solve Equation 14 by using Adomian decomposition method (ADM). Suppose g1g_1 can be written in a series g1=g1(0)+g1(1)+⋯ ,g_1=g_1^{(0)}+g_1^{(1)}+\cdots, and substitute it into Equation 14, and we have g1(0)(y,κ)+⋯=−ln⁡y(yκ−1+1+∫κ1(g(0) ⁣(yκv,v)+⋯ )dvv).g_1^{(0)}(y,\kappa)+\cdots =-\ln y\left(y^{\kappa-1}+1+\int_\kappa^1\left( g^{(0)}\!\left(y^{\frac\kappa v},v\right)+\cdots \right)\frac{\mathrm dv}v\right). Assume we may interchange integration and summation (which is OK here because we can verify the solution after we find it using ADM). Then, = g1(0)(y,κ)+g1(1)(y,κ)+⋯=−ln⁡y(yκ−1+1)−ln⁡y∫κ1g(0) ⁣(yκv,v)dvv−⋯ .\begin{align*} &\phantom{=~}g_1^{(0)}(y,\kappa)+g_1^{(1)}(y,\kappa)+\cdots\\ &=-\ln y\left(y^{\kappa-1}+1\right) -\ln y\int_\kappa^1g^{(0)}\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v-\cdots. \end{align*} If we let g1(0)(y,κ)≔−ln⁡y(yκ−1+1),g1(i+1)(y,κ)≔−ln⁡y∫κ1g(i) ⁣(yκv,v)dvv,i∈N,\begin{split} g_1^{(0)}(y,\kappa)&\coloneqq-\ln y\left(y^{\kappa-1}+1\right),\\ g_1^{(i+1)}(y,\kappa)&\coloneqq-\ln y\int_\kappa^1g^{(i)}\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v,\quad i\in\mathbb N, \end{split} (15)(15) then we can equate each term in the two series. If the sum g1=∑i=0∞g1(i)g_1=\sum_{i=0}^\infty g_1^{(i)} converges, then this is a guess of the solution to Equation 14, which we can verify whether it is correct or not.

Using Equation 15, we can find first few terms in the series by directly integrating. The first few terms are g1(0)(y,κ)=−ln⁡y(yκ−1+1),g1(1)(y,κ)=−ln⁡y(yκ−1−1+ln⁡yκ−1),g1(2)(y,κ)=−ln⁡y(yκ−1−1−ln⁡yκ−1+12(ln⁡yκ−1)2),g1(3)(y,κ)=−ln⁡y(yκ−1−1−ln⁡yκ−1−12(ln⁡yκ−1)2+16(ln⁡yκ−1)3),⋮\begin{split} g_1^{(0)}(y,\kappa)&=-\ln y\left(y^{\kappa-1}+1\right),\\ g_1^{(1)}(y,\kappa)&=-\ln y\left(y^{\kappa-1}-1+\ln y^{\kappa-1}\right),\\ g_1^{(2)}(y,\kappa)&=-\ln y\left(y^{\kappa-1}-1-\ln y^{\kappa-1}+\frac12\left(\ln y^{\kappa-1}\right)^2\right),\\ g_1^{(3)}(y,\kappa)&=-\ln y\left(y^{\kappa-1}-1-\ln y^{\kappa-1}-\frac12\left(\ln y^{\kappa-1}\right)^2+\frac16\left(\ln y^{\kappa-1}\right)^3\right),\\ \vdots& \end{split} We may then guess that the terms have general formula g1(i)(y,κ)=−ln⁡y(yκ−1+1i!(ln⁡yκ−1)i−∑j=0i−11j!(ln⁡yκ−1)j).g_1^{(i)}(y,\kappa)=-\ln y\left(y^{\kappa-1}+\frac1{i!}\left(\ln y^{\kappa-1}\right)^i -\sum_{j=0}^{i-1}\frac1{j!}\left(\ln y^{\kappa-1}\right)^j\right). Sum up the terms, and we have g1(y,κ)=∑i=0∞g1(i)(y,κ)=lim⁡q→∞∑i=0q−ln⁡y(yκ−1+1i!(ln⁡yκ−1)i−∑j=0i−11j!(ln⁡yκ−1)j)=−ln⁡y(exp⁡ln⁡yκ−1+lim⁡q→∞((q+1)yκ−1−∑j=0q∑i=j+1q1j!(ln⁡yκ−1)j))=−ln⁡y(yκ−1+lim⁡q→∞((q+1)yκ−1−∑j=0qq−jj!(ln⁡yκ−1)j))=−ln⁡y(yκ−1+lim⁡q→∞(qyκ−1−q∑j=0q1j!(ln⁡yκ−1)j)=−ln⁡y∑jq+yκ−1+ln⁡yκ−1exp⁡ln⁡yκ−1)=−ln⁡y(2+ln⁡yκ−1)yκ−1.\begin{align*} g_1(y,\kappa)&=\sum_{i=0}^\infty g_1^{(i)}(y,\kappa)\\ &=\lim_{q\to\infty}\sum_{i=0}^q-\ln y\left( y^{\kappa-1}+\frac1{i!}\left(\ln y^{\kappa-1}\right)^i -\sum_{j=0}^{i-1}\frac1{j!}\left(\ln y^{\kappa-1}\right)^j \right)\\ &=-\ln y\left(\exp\ln y^{\kappa-1}+\lim_{q\to\infty}\left( \left(q+1\right)y^{\kappa-1} -\sum_{j=0}^q\sum_{i=j+1}^q\frac1{j!}\left(\ln y^{\kappa-1}\right)^j \right)\right)\\ &=-\ln y\left(y^{\kappa-1}+\lim_{q\to\infty}\left( \left(q+1\right)y^{\kappa-1} -\sum_{j=0}^q\frac{q-j}{j!}\left(\ln y^{\kappa-1}\right)^j \right)\right)\\ &=-\ln y\left( y^{\kappa-1}+\lim_{q\to\infty}\left( qy^{\kappa-1}-q\sum_{j=0}^q\frac1{j!}\left(\ln y^{\kappa-1}\right)^j \right) \right.\\&\phantom{=-\ln y}\left.\vphantom{\sum_j^q} +y^{\kappa-1} +\ln y^{\kappa-1}\exp\ln y^{\kappa-1} \right)\\ &=-\ln y\left(2+\ln y^{\kappa-1}\right)y^{\kappa-1}. \end{align*}

Therefore, we have the final guess of solution g1(y,κ)=−ln⁡y(2+ln⁡yκ−1)yκ−1.g_1(y,\kappa)=-\ln y\left(2+\ln y^{\kappa-1}\right)y^{\kappa-1}. (16)(16) We can substitute it into Equation 14 to verify that it is indeed the solution.

The case κ∈(13,12)\kappa\in\left(\frac13,\frac12\right)

In this case, we have min⁡ ⁣(κ1−κ,1)=κ1−κ∈(12,1).\min\!\left(\frac\kappa{1-\kappa},1\right)=\frac\kappa{1-\kappa}\in\left(\frac12,1\right). We can then use the same method as in the previous case to find the solution.

First, by Equation 14, g2(y,κ)≔g ⁣(y,κ∈(13,12))=−ln⁡y(∫0κ1−κg ⁣(y1−κ,u)du+∫κκ1−κg ⁣(yκv,v)dvv)=−ln⁡y(∫01g ⁣(y1−κ,u)du−∫κ1−κ1g1 ⁣(y1−κ,u)du=−ln⁡y+∫κ12g2 ⁣(yκv,v)dvv+∫12κ1−κg1 ⁣(yκv,v)dvv).\begin{align*} g_2(y,\kappa)&\coloneqq g\!\left(y,\kappa\in\left(\frac13,\frac12\right)\right)\\ &=-\ln y\left( \int_0^{\frac\kappa{1-\kappa}}g\!\left(y^{1-\kappa},u\right)\mathrm du +\int_\kappa^{\frac\kappa{1-\kappa}}g\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v \right)\\ &=-\ln y\left( \int_0^1g\!\left(y^{1-\kappa},u\right)\mathrm du -\int_{\frac\kappa{1-\kappa}}^1g_1\!\left(y^{1-\kappa},u\right)\mathrm du \right.\\&\phantom{=-\ln y}\left. +\int_\kappa^{\frac12}g_2\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v +\int_{\frac12}^{\frac\kappa{1-\kappa}}g_1\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v \right). \end{align*} Substitute Equation 13 and 16 into the above equation, and we have g2(y,κ)=−ln⁡y(yκ−1+y−κ(1+ln⁡y−κ)−2y2κ−1(1+ln⁡y2κ−1)∫κ12=−ln⁡y+∫κ12g2 ⁣(yκv,v)dvv).\begin{split} g_2(y,\kappa)&=-\ln y\left( y^{\kappa-1} +y^{-\kappa}\left(1+\ln y^{-\kappa}\right) -2y^{2\kappa-1}\left(1+\ln y^{2\kappa-1}\right) \vphantom{\int_\kappa^{\frac12}}\right.\\&\phantom{=-\ln y}\left. +\int_\kappa^{\frac12}g_2\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v \right). \end{split} (17)(17)

Equation 17 can again be solved by ADM though the calculation is much more complicated than the previous case. We may guess g2=∑i=0∞g2(i)g_2=\sum_{i=0}^\infty g_2^{(i)} is the solution if the series converges, where g2(0)(y,κ)≔−ln⁡y(yκ−1+y−κ(1+ln⁡y−κ)−2y2κ−1(1+ln⁡y2κ−1)),g1(i+1)(y,κ)≔−ln⁡y∫κ12g(i) ⁣(yκv,v)dvv,i∈N.\begin{split} g_2^{(0)}(y,\kappa)&\coloneqq-\ln y\left( y^{\kappa-1} +y^{-\kappa}\left(1+\ln y^{-\kappa}\right) -2y^{2\kappa-1}\left(1+\ln y^{2\kappa-1}\right) \right),\\ g_1^{(i+1)}(y,\kappa)&\coloneqq-\ln y\int_\kappa^{\frac12}g^{(i)}\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v, \quad i\in\mathbb N. \end{split} The first few terms go too long to be written here before one may find the pattern, so they are omitted here. If you want to see them, use a mathematical software to help you, and you should be able to find the pattern after calculating first six (or so) terms. After looking at first few terms, the guessed general term is g2(i)=−ln⁡y (yκ−1+2y2κ−1(i−1−ln⁡y2κ−1)−2∑j=0i−1i−j−1j!(ln⁡y2κ−1)j−y−κ∑j=0i−11j!(ln⁡y2κ−1)j+1i!y−κ(1+ln⁡y−κ)(ln⁡y2κ−1)i).\begin{align*} g_2^{(i)}=-\ln y&\,\left( y^{\kappa-1}+2y^{2\kappa-1}\left(i-1-\ln y^{2\kappa-1}\right) -2\sum_{j=0}^{i-1}\frac{i-j-1}{j!}\left(\ln y^{2\kappa-1}\right)^j \right.\\&\left. {}-y^{-\kappa}\sum_{j=0}^{i-1}\frac1{j!}\left(\ln y^{2\kappa-1}\right)^j +\frac1{i!}y^{-\kappa}\left(1+\ln y^{-\kappa}\right)\left(\ln y^{2\kappa-1}\right)^i \right). \end{align*} Then we can sum it to get a guess of g2g_2.

After some tedious calculation, we have g2(y,κ)=−ln⁡y((2+ln⁡yκ−1)yκ−1−(2+4ln⁡y2κ−1+(ln⁡y2κ−1)2)y2κ−1).g_2(y,\kappa)=-\ln y\left( \left(2+\ln y^{\kappa-1}\right)y^{\kappa-1} -\left(2+4\ln y^{2\kappa-1}+\left(\ln y^{2\kappa-1}\right)^2\right)y^{2\kappa-1} \right). (18)(18) On may verify that this is indeed the solution by substituting it into Equation 17.

The case κ∈(14,13)\kappa\in\left(\frac14,\frac13\right)

By using very similar methods but after very tedious calculation, the solution is g3(y,κ)≔g ⁣(y,κ∈(14,13))=−ln⁡y((2+ln⁡yκ−1)yκ−1−(2+4ln⁡y2κ−1+(ln⁡y2κ−1)2)y2κ−112=−ln⁡y+(3ln⁡y3κ−1+3(ln⁡y3κ−1)2+12(ln⁡y3κ−1)3)y3κ−1).\begin{split} g_3(y,\kappa)&\coloneqq g\!\left(y,\kappa\in\left(\frac14,\frac13\right)\right)\\ &=-\ln y\left( \left(2+\ln y^{\kappa-1}\right)y^{\kappa-1} -\left(2+4\ln y^{2\kappa-1}+\left(\ln y^{2\kappa-1}\right)^2\right)y^{2\kappa-1}\vphantom{\frac12} \right.\\&\phantom{=-\ln y}\left. {}+\left(3\ln y^{3\kappa-1}+3\left(\ln y^{3\kappa-1}\right)^2+\frac12\left(\ln y^{3\kappa-1}\right)^3\right)y^{3\kappa-1} \right). \end{split} (19)(19)

Other cases

After seeing Equation 16, 18, and 19, one may guess the form of solution for other cases.

Guess the form of solution for κ∈(1q+1,1q)\kappa\in\left(\frac1{q+1},\frac1q\right), where q∈1…∞q\in1\ldots\infty, is gq(y,κ)≔g ⁣(y,κ∈(1q,1q+1))=∑s=1qΔgs(y,κ),g_q(y,\kappa)\coloneqq g\!\left(y,\kappa\in\left(\frac1q,\frac1{q+1}\right)\right)=\sum_{s=1}^q\Delta g_s(y,\kappa), where Δgs(y,κ)≔(−1)sysκ−1ln⁡y∑j=0sAs,jj!(ln⁡ysκ−1)j,\Delta g_s(y,\kappa)\coloneqq\left(-1\right)^sy^{s\kappa-1}\ln y\sum_{j=0}^s\frac{A_{s,j}}{j!}\left(\ln y^{s\kappa-1}\right)^j, where As,jA_{s,j} are coefficients to be determined.

Now, consider the cases q∈2…∞q\in2\ldots\infty. Because κ∈(1q+1,1q)\kappa\in\left(\frac1{q+1},\frac1q\right), min⁡ ⁣(κ1−κ,1)=κ1−κ∈(1q,1q−1).\min\!\left(\frac\kappa{1-\kappa},1\right)=\frac\kappa{1-\kappa}\in\left(\frac1q,\frac1{q-1}\right). Therefore, = ∫0min⁡(κ1−κ,1)(g ⁣(y1−κ,u)+δ(u−1))du=∫01g ⁣(y1−κ,u)du−∑p=1q−2∫1p+11pgp ⁣(y1−κ,u)du−∫κ1−κ1q−1gq−1 ⁣(y1−κ,u)du=yκ−1−1−∑p=1q−2∫1p+11p∑s=1pΔgs ⁣(y1−κ,u)du−∫κ1−κ1q−1∑s=1q−1Δgs ⁣(y1−κ,u)du=yκ−1−1−∑s=1q−1(∑p=sq−2∫1p+11p+∫κ1−κ1q−1)Δgs ⁣(y1−κ,u)du=yκ−1−1−∑s=1q−1∫κ1−κ1sΔgs ⁣(y1−κ,u)du,\begin{align*} &\phantom{=~}\int_0^{\min\left(\frac\kappa{1-\kappa},1\right)}\left(g\!\left(y^{1-\kappa},u\right)+\delta(u-1)\right)\mathrm du\\ &=\int_0^1g\!\left(y^{1-\kappa},u\right)\mathrm du -\sum_{p=1}^{q-2}\int_{\frac1{p+1}}^{\frac1p}g_p\!\left(y^{1-\kappa},u\right)\mathrm du -\int_{\frac\kappa{1-\kappa}}^\frac1{q-1}g_{q-1}\!\left(y^{1-\kappa},u\right)\mathrm du\\ &=y^{\kappa-1}-1 -\sum_{p=1}^{q-2}\int_{\frac1{p+1}}^{\frac1p}\sum_{s=1}^p\Delta g_s\!\left(y^{1-\kappa},u\right)\mathrm du -\int_{\frac\kappa{1-\kappa}}^\frac1{q-1}\sum_{s=1}^{q-1}\Delta g_s\!\left(y^{1-\kappa},u\right)\mathrm du\\ &=y^{\kappa-1}-1 -\sum_{s=1}^{q-1}\left( \sum_{p=s}^{q-2}\int_{\frac1{p+1}}^{\frac1p}+\int_{\frac\kappa{1-\kappa}}^\frac1{q-1} \right)\Delta g_s\!\left(y^{1-\kappa},u\right)\mathrm du\\ &=y^{\kappa-1}-1 -\sum_{s=1}^{q-1}\int_{\frac\kappa{1-\kappa}}^{\frac1s}\Delta g_s\!\left(y^{1-\kappa},u\right)\mathrm du, \end{align*} and = ∫κmin⁡(κ1−κ,1)(g ⁣(yκv,v)+δ(v−1))dvv=∫κ1qgq ⁣(yκv,v)dvv+∫1qκ1−κgq−1 ⁣(yκv,v)dvv=∫κ1q∑s=1qΔgs ⁣(yκv,v)dvv+∫1qκ1−κ∑s=1q−1Δgs ⁣(yκv,v)dvv=∑s=1q−1∫κκ1−κΔgs ⁣(yκv,v)dvv+∫κ1qΔgq ⁣(yκv,v)dvv.\begin{align*} &\phantom{=~}\int_\kappa^{\min\left(\frac\kappa{1-\kappa},1\right)}\left(g\!\left(y^{\frac\kappa v},v\right)+\delta(v-1)\right)\frac{\mathrm dv}v\\ &=\int_\kappa^{\frac1q}g_q\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v +\int_{\frac1q}^{\frac\kappa{1-\kappa}}g_{q-1}\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v\\ &=\int_\kappa^{\frac1q}\sum_{s=1}^q\Delta g_s\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v +\int_{\frac1q}^{\frac\kappa{1-\kappa}}\sum_{s=1}^{q-1}\Delta g_s\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v\\ &=\sum_{s=1}^{q-1}\int_\kappa^{\frac\kappa{1-\kappa}}\Delta g_s\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v +\int_\kappa^{\frac1q}\Delta g_q\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v. \end{align*} Substitute into Equation 12, and we have ∑s=1qΔgs(y,κ)=−ln⁡y (yκ−1−1−∑s=1q−1∫κ1−κ1sΔgs ⁣(y1−κ,u)du+∑s=1q−1∫κκ1−κΔgs ⁣(yκv,v)dvv+∫κ1qΔgq ⁣(yκv,v)dvv).\begin{split} \sum_{s=1}^q\Delta g_s(y,\kappa)=-\ln y&\,\left( y^{\kappa-1}-1 -\sum_{s=1}^{q-1}\int_{\frac\kappa{1-\kappa}}^{\frac1s}\Delta g_s\!\left(y^{1-\kappa},u\right)\mathrm du \right.\\&\left. {}+\sum_{s=1}^{q-1}\int_\kappa^{\frac\kappa{1-\kappa}}\Delta g_s\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v +\int_\kappa^{\frac1q}\Delta g_q\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v \right). \end{split} (20)(20) To simplify later expressions, define Bs,l≔(−1)l∑j=ls(−1)jAs,j.B_{s,l}\coloneqq\left(-1\right)^l\sum_{j=l}^s\left(-1\right)^jA_{s,j}. (21)(21)

Now, calculate the integrals in Equation 20. Before that, first we introduce a handy integral formula: ∫(ln⁡w)j dw=(−1)jj! w∑l=0j(−1)l(ln⁡w)ll!+C.\int\left(\ln w\right)^j\,\mathrm dw =\left(-1\right)^jj!\,w\sum_{l=0}^j\left(-1\right)^l\frac{\left(\ln w\right)^l}{l!}+C. This formula can be proved by mathematical induction and integration by parts.

Then, we have = ∫κ1−κ1sΔgs ⁣(y1−κ,u)du=∫κ1−κ1s(−1)sy(su−1)(1−κ)ln⁡y1−κ∑j=0sAs,jj!(ln⁡y(su−1)(1−κ))j du=(−1)ss∑j=0sAs,jj!∫κ1−κ1s(ln⁡y(su−1)(1−κ))j d(y(su−1)(1−κ))=(−1)ss∑j=0sAs,jj!∫y(s+1)κ−11(ln⁡w)j dw=(−1)ss∑j=0sAs,jj!(−1)jj!(1−y(s+1)κ−1∑l=0j(−1)l(ln⁡y(s+1)κ−1)ll!)=(−1)ss∑j=0s(−1)jAs,j−(−1)ssy(s+1)κ−1∑l=0s(−1)l(ln⁡y(s+1)κ−1)ll!∑j=ls(−1)jAs,j=(−1)s(Bs,0s−1sy(s+1)κ−1∑l=0sBs,ll!(ln⁡y(s+1)κ−1)l).\begin{align*} &\phantom{=~}\int_{\frac\kappa{1-\kappa}}^{\frac1s}\Delta g_s\!\left(y^{1-\kappa},u\right)\mathrm du\\ &=\int_{\frac\kappa{1-\kappa}}^{\frac1s}\left(-1\right)^s y^{\left(su-1\right)\left(1-\kappa\right)}\ln y^{1-\kappa} \sum_{j=0}^s\frac{A_{s,j}}{j!}\left(\ln y^{\left(su-1\right)\left(1-\kappa\right)}\right)^j\,\mathrm du\\ &=\frac{\left(-1\right)^s}s\sum_{j=0}^s\frac{A_{s,j}}{j!} \int_{\frac\kappa{1-\kappa}}^{\frac1s}\left(\ln y^{\left(su-1\right)\left(1-\kappa\right)}\right)^j\, \mathrm d\left(y^{\left(su-1\right)\left(1-\kappa\right)}\right)\\ &=\frac{\left(-1\right)^s}s\sum_{j=0}^s\frac{A_{s,j}}{j!} \int_{y^{\left(s+1\right)\kappa-1}}^1\left(\ln w\right)^j\,\mathrm dw\\ &=\frac{\left(-1\right)^s}s\sum_{j=0}^s\frac{A_{s,j}}{j!} \left(-1\right)^jj!\left( 1 -y^{\left(s+1\right)\kappa-1} \sum_{l=0}^j\left(-1\right)^l\frac{\left(\ln y^{\left(s+1\right)\kappa-1}\right)^l}{l!} \right)\\ &=\frac{\left(-1\right)^s}s\sum_{j=0}^s\left(-1\right)^jA_{s,j} -\frac{\left(-1\right)^s}sy^{\left(s+1\right)\kappa-1} \sum_{l=0}^s\left(-1\right)^l\frac{\left(\ln y^{\left(s+1\right)\kappa-1}\right)^l}{l!} \sum_{j=l}^s\left(-1\right)^jA_{s,j}\\ &=\left(-1\right)^s\left( \frac{B_{s,0}}s -\frac1sy^{\left(s+1\right)\kappa-1} \sum_{l=0}^s\frac{B_{s,l}}{l!} \left(\ln y^{\left(s+1\right)\kappa-1}\right)^l \right). \end{align*} = ∫κκ1−κΔgs ⁣(yκv,v)dvv=∫κκ1−κ(−1)syκv(sv−1)ln⁡yκv∑j=0sAs,jj!(ln⁡yκv(sv−1))jdvv=(−1)s∑j=0sAs,jj!∫κκ1−κ(ln⁡yκv(sv−1))j d(yκv(sv−1))=(−1)s∑j=0sAs,jj!∫ysκ−1y(s+1)κ−1(ln⁡w)j dw=(−1)s∑j=0sAs,jj!(−1)jj!(y(s+1)κ−1∑l=0j(−1)l(ln⁡y(s+1)κ−1)ll!=(−1)s∑j=0sAs,jj!(−1)jj!−ysκ−1∑l=0j(−1)l(ln⁡ysκ−1)ll!)=(−1)s(y(s+1)κ−1∑l=0sBs,ll!(ln⁡y(s+1)κ−1)l−ysκ−1∑l=0sBs,ll!(ln⁡ysκ−1)l).\begin{align*} &\phantom{=~}\int_\kappa^{\frac\kappa{1-\kappa}}\Delta g_s\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v\\ &=\int_\kappa^{\frac\kappa{1-\kappa}}\left(-1\right)^s y^{\frac\kappa v\left(sv-1\right)}\ln y^{\frac\kappa v} \sum_{j=0}^s\frac{A_{s,j}}{j!}\left(\ln y^{\frac\kappa v\left(sv-1\right)}\right)^j\frac{\mathrm dv}v\\ &=\left(-1\right)^s\sum_{j=0}^s\frac{A_{s,j}}{j!} \int_\kappa^{\frac\kappa{1-\kappa}}\left(\ln y^{\frac\kappa v\left(sv-1\right)}\right)^j\, \mathrm d\left(y^{\frac\kappa v\left(sv-1\right)}\right)\\ &=\left(-1\right)^s\sum_{j=0}^s\frac{A_{s,j}}{j!} \int_{y^{s\kappa-1}}^{y^{\left(s+1\right)\kappa-1}}\left(\ln w\right)^j\,\mathrm dw\\ &=\left(-1\right)^s\sum_{j=0}^s\frac{A_{s,j}}{j!}\left(-1\right)^jj!\left( y^{\left(s+1\right)\kappa-1}\sum_{l=0}^j\left(-1\right)^l\frac{\left(\ln y^{\left(s+1\right)\kappa-1}\right)^l}{l!} \right.\\&\phantom{=\left(-1\right)^s\sum_{j=0}^s\frac{A_{s,j}}{j!}\left(-1\right)^jj!}\left. {}-y^{s\kappa-1}\sum_{l=0}^j\left(-1\right)^l\frac{\left(\ln y^{s\kappa-1}\right)^l}{l!} \right)\\ &=\left(-1\right)^s\left( y^{\left(s+1\right)\kappa-1}\sum_{l=0}^s\frac{B_{s,l}}{l!} \left(\ln y^{\left(s+1\right)\kappa-1}\right)^l -y^{s\kappa-1}\sum_{l=0}^s\frac{B_{s,l}}{l!} \left(\ln y^{s\kappa-1}\right)^l \right). \end{align*} = ∫κ1qΔgq ⁣(yκv,v)dvv=(−1)q(Bq,0−yqκ−1∑l=0qBq,ll!(ln⁡yqκ−1)l).\begin{align*} &\phantom{=~}\int_\kappa^{\frac1q}\Delta g_q\!\left(y^{\frac\kappa v},v\right)\frac{\mathrm dv}v\\ &=\left(-1\right)^q\left( B_{q,0} -y^{q\kappa-1}\sum_{l=0}^q\frac{B_{q,l}}{l!}\left(\ln y^{q\kappa-1}\right)^l \right). \end{align*} Substitute these results into Equation 20, and we have = ∑s=1q(−1)sysκ−1ln⁡y∑j=0sAs,jj!(ln⁡ysκ−1)j=−ln⁡y(∑lsyκ−1−1=−ln⁡y−∑s=1q−1(−1)s(Bs,0s−1sy(s+1)κ−1∑l=0sBs,ll!(ln⁡y(s+1)κ−1)l)=−ln⁡y+∑s=1q−1(−1)s(y(s+1)κ−1∑l=0sBs,ll!(ln⁡y(s+1)κ−1)l−ysκ−1∑l=0sBs,ll!(ln⁡ysκ−1)l)=−ln⁡y+(−1)q(Bq,0−yqκ−1∑l=0qBq,ll!(ln⁡yqκ−1)l)).\begin{split} &\phantom{=~}\sum_{s=1}^q \left(-1\right)^sy^{s\kappa-1}\ln y\sum_{j=0}^s\frac{A_{s,j}}{j!}\left(\ln y^{s\kappa-1}\right)^j\\ &=-\ln y\left(\vphantom{\sum_l^s} y^{\kappa-1}-1 \right.\\&\phantom{=-\ln y}\left. {}-\sum_{s=1}^{q-1}\left(-1\right)^s\left( \frac{B_{s,0}}s -\frac1sy^{\left(s+1\right)\kappa-1} \sum_{l=0}^s\frac{B_{s,l}}{l!} \left(\ln y^{\left(s+1\right)\kappa-1}\right)^l \right) \right.\\&\phantom{=-\ln y}\left. {}+\sum_{s=1}^{q-1}\left(-1\right)^s\left( y^{\left(s+1\right)\kappa-1}\sum_{l=0}^s\frac{B_{s,l}}{l!} \left(\ln y^{\left(s+1\right)\kappa-1}\right)^l -y^{s\kappa-1}\sum_{l=0}^s\frac{B_{s,l}}{l!} \left(\ln y^{s\kappa-1}\right)^l \right) \right.\\&\phantom{=-\ln y}\left. {}+\left(-1\right)^q\left( B_{q,0} -y^{q\kappa-1}\sum_{l=0}^q\frac{B_{q,l}}{l!}\left(\ln y^{q\kappa-1}\right)^l \right) \right). \end{split} Cancel factor ln⁡y\ln y on both sides, and we have = ∑s=1q(−1)sysκ−1∑j=0sAs,jj!(ln⁡ysκ−1)j=−yκ−1+1=+∑s=1q−1(−1)sBs,0s−∑s=2q(−1)s−1s−1ysκ−1∑l=0s−1Bs−1,l(ln⁡ysκ−1)ll!=−∑s=2q(−1)s−1ysκ−1∑l=0s−1Bs−1,l(ln⁡ysκ−1)ll!+∑s=1q−1(−1)sysκ−1∑l=0sBs,l(ln⁡ysκ−1)ll!=−(−1)qBq,0+(−1)qyqκ−1∑l=0qBq,l(ln⁡yqκ−1)ll!=1+∑s=1q−1(−1)ssBs,0−(−1)qBq,0(∗)=−yκ−1(1+B1,0−B1,1ln⁡yκ−1)(∗∗)=+∑s=2q(−1)sysκ−1(∑l=0s−1(ss−1Bs−1,l+Bs,l)(ln⁡ysκ−1)ll!+Bs,s(ln⁡ysκ−1)ss!)(∗∗∗)\begin{align*} &\phantom{=~}\sum_{s=1}^q \left(-1\right)^sy^{s\kappa-1}\sum_{j=0}^s\frac{A_{s,j}}{j!}\left(\ln y^{s\kappa-1}\right)^j\\ &=-y^{\kappa-1}+1\\ &\phantom{=}{}+\sum_{s=1}^{q-1}\left(-1\right)^s\frac{B_{s,0}}s -\sum_{s=2}^q\frac{\left(-1\right)^{s-1}}{s-1}y^{s\kappa-1}\sum_{l=0}^{s-1}B_{s-1,l}\frac{\left(\ln y^{s\kappa-1}\right)^l}{l!}\\ &\phantom{=}{}-\sum_{s=2}^q\left(-1\right)^{s-1}y^{s\kappa-1}\sum_{l=0}^{s-1}B_{s-1,l}\frac{\left(\ln y^{s\kappa-1}\right)^l}{l!} +\sum_{s=1}^{q-1}\left(-1\right)^sy^{s\kappa-1}\sum_{l=0}^sB_{s,l}\frac{\left(\ln y^{s\kappa-1}\right)^l}{l!}\\ &\phantom{=}{}-\left(-1\right)^qB_{q,0}+\left(-1\right)^qy^{q\kappa-1}\sum_{l=0}^qB_{q,l}\frac{\left(\ln y^{q\kappa-1}\right)^l}{l!}\\ &=1+\sum_{s=1}^{q-1}\frac{\left(-1\right)^s}sB_{s,0}-\left(-1\right)^qB_{q,0}&(*)\\ &\phantom{=}-y^{\kappa-1}\left(1+B_{1,0}-B_{1,1}\ln y^{\kappa-1}\right)&(**)\\ &\phantom{=}{}+\sum_{s=2}^q\left(-1\right)^sy^{s\kappa-1}\left( \sum_{l=0}^{s-1}\left(\frac s{s-1}B_{s-1,l}+B_{s,l}\right)\frac{\left(\ln y^{s\kappa-1}\right)^l}{l!} +B_{s,s}\frac{\left(\ln y^{s\kappa-1}\right)^s}{s!} \right)&(***) \end{align*}

Equate the coefficients in Line (*) with the corresponding ones on the LHS, and we have 0=1+∑s=1q−1(−1)ssBs,0−(−1)qBq,0.0=1+\sum_{s=1}^{q-1}\frac{\left(-1\right)^s}sB_{s,0}-\left(-1\right)^qB_{q,0}. This equation holds for any q∈2…∞q\in 2\ldots\infty, so {(−1)qBq,0=1+∑s=1q−1(−1)ssBs,0,(−1)q+1Bq+1,0=1+∑s=1q(−1)ssBs,0,q∈2…∞.\begin{cases} \left(-1\right)^qB_{q,0}=1+\sum_{s=1}^{q-1}\frac{\left(-1\right)^s}sB_{s,0},\\ \left(-1\right)^{q+1}B_{q+1,0}=1+\sum_{s=1}^{q}\frac{\left(-1\right)^s}sB_{s,0}, \end{cases}\quad q\in 2\ldots\infty. Subtract the two equations, and we have Bq+1,0=−q+1qBq,0,q∈2…∞.B_{q+1,0}=-\frac{q+1}qB_{q,0},\quad q\in 2\ldots\infty. (22)(22) Equation 22 can determine Bq,0B_{q,0} for all q∈2…∞q\in 2\ldots\infty once B2,0B_{2,0} is determined. The relationship between B1,0B_{1,0} and B2,0B_{2,0} cannot be described by Equation 22, but is given by B2,0=1−B1,0.B_{2,0}=1-B_{1,0}. (23)(23)

Equate the coefficients in Line (**) with the corresponding ones on the LHS, and we have A1,0=1+B1,0,A1,1=B1,1.A_{1,0}=1+B_{1,0},\quad A_{1,1}=B_{1,1}. By Equation 21, this is equivalent to A1,0=1+A1,0−A1,1,A1,1=A1,1.A_{1,0}=1+A_{1,0}-A_{1,1},\quad A_{1,1}=A_{1,1}. Therefore, A1,1=1A_{1,1}=1, and thus B1,1=1.B_{1,1}=1. (24)(24)

Equate the coefficients in Line (***) with the corresponding ones on the LHS, and we have As,l=ss−1Bs−1,l+Bs,l,As,s=Bs,s.A_{s,l}=\frac s{s-1}B_{s-1,l}+B_{s,l},\quad A_{s,s}=B_{s,s}. By Equation 21, Bs,l=As,l−Bs,l+1B_{s,l}=A_{s,l}-B_{s,l+1} for l∈0…sl\in 0\ldots s, and As,s=Bs,sA_{s,s}=B_{s,s} is always true. Therefore, 0=ss−1Bs−1,l−Bs,l+1.0=\frac s{s-1}B_{s-1,l}-B_{s,l+1}. This equation is true for any s∈2 ..qs\in 2\,..q and l∈0…sl\in0\ldots s. Because qq is arbitrary, we can change the variable ss to qq and the equation tells us exactly the same information. Therefore, Bq,l=qq−1Bq−1,l−1,q∈2…∞,l∈1 ..q.B_{q,l}=\frac q{q-1}B_{q-1,l-1},\quad q\in 2\ldots\infty,\quad l\in 1\,..q. (25)(25)

Equation 22, 23, 24, and 25 are sufficient to determine Bq,lB_{q,l} for all q∈1…∞q\in 1\ldots\infty and l∈0 ..ql\in 0\,..q up to one arbitrary parameter. Define the arbitrary parameter b≔1−B1,0,b\coloneqq1-B_{1,0}, then the first few Bq,lB_{q,l} are l=01234q=11−b122b2−2b23−3b3b3−3b344b−4b4b4−4b4⋮⋱\begin{array}{r|ccccc} &l=0&1&2&3&4\\ \hline q=1&1-b&1\\ 2&2b&2-2b&2\\ 3&-3b&3b&3-3b&3\\ 4&4b&-4b&4b&4-4b&4\\ \vdots&\ddots \end{array} The general formula for Bq,lB_{q,l} is Bq,l={q,l=q,q(1−b),l=q−1,(−1)q+lqb,l∈0 ..q−2,B_{q,l}=\begin{cases} q,&l=q,\\ q\left(1-b\right),&l=q-1,\\ \left(-1\right)^{q+l}qb,&l\in 0\,..q-2, \end{cases} which may be proved by mathematical induction.

Actually, one may find b=0b=0 by simply comparing with the results in Equation 16, 18, or 19. Another way to find bb is comparing with Eqution 13. Here I wil show the latter approach. = ∫01g ⁣(y,κ)dκ=∑q=1∞∫1q+11qgq ⁣(y,κ)dκ=∑q=1∞∫1q+11q∑s=1qΔgs ⁣(y,κ)dκ=∑s=1∞∑q=s∞∫1q+11qΔgs ⁣(y,κ)dκ=∑s=1∞∫01sΔgs ⁣(y,κ)dκ=∑s=1∞∫01s(−1)sysκ−1ln⁡y∑j=0sAs,jj!(ln⁡ysκ−1)j dκ=∑s=1∞(−1)ss∑j=0sAs,jj!∫01s(ln⁡ysκ−1)j d(ysκ−1)=∑s=1∞(−1)ss(Bs,0−y−1∑l=0sBs,l(ln⁡y)ll!)=−(1−b)+y−1((1−b)−ln⁡y)=+∑s=2∞(−1)ss((−1)ssb−y−1(∑l=0s−2(−1)s+lsb(ln⁡y)ll!=+∑s=2∞(−1)ss (∑s∞(−1)ssb−y−1+s(1−b)(−ln⁡y)s−1(s−1)!+s(−ln⁡y)ss!))=b−1+y−1(1−b−ln⁡y)+y−1∑s=2∞((ln⁡y)s−1(s−1)!−(ln⁡y)ss!)=+b∑s=2∞(1−y−1∑l=0s−1(ln⁡y)ll!)=b−1+y−1(1−b−ln⁡y)+y−1((exp⁡ln⁡y−1)−(exp⁡ln⁡y−ln⁡y−1))=+blim⁡q→∞∑s=2q(1−y−1∑l=0s−1(ln⁡y)l)=b−1+y−1(1−b)+blim⁡q→∞(q−1−y−1(q−1+∑l=1q−1q−ll!(ln⁡y)l))=b−1+y−1(1−b)+b(−1+y−1+ln⁡y)=y−1−1+bln⁡y.\begin{align*} &\phantom{=~}\int_0^1g\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{q=1}^\infty\int_{\frac1{q+1}}^{\frac1q}g_q\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{q=1}^\infty\int_{\frac1{q+1}}^{\frac1q}\sum_{s=1}^q\Delta g_s\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{s=1}^\infty\sum_{q=s}^\infty\int_{\frac1{q+1}}^{\frac1q}\Delta g_s\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{s=1}^\infty\int_0^{\frac1s}\Delta g_s\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{s=1}^\infty\int_0^{\frac1s}\left(-1\right)^sy^{s\kappa-1}\ln y \sum_{j=0}^s\frac{A_{s,j}}{j!}\left(\ln y^{s\kappa-1}\right)^j\,\mathrm d\kappa\\ &=\sum_{s=1}^\infty\frac{\left(-1\right)^s}s\sum_{j=0}^s\frac{A_{s,j}}{j!} \int_0^{\frac1s}\left(\ln y^{s\kappa-1}\right)^j\,\mathrm d\left(y^{s\kappa-1}\right)\\ &=\sum_{s=1}^\infty\frac{\left(-1\right)^s}s\left( B_{s,0}-y^{-1}\sum_{l=0}^sB_{s,l}\frac{\left(\ln y\right)^l}{l!} \right)\\ &=-\left(1-b\right)+y^{-1}\left(\left(1-b\right)-\ln y\right)\\ &\phantom{=}{}+\sum_{s=2}^\infty\frac{\left(-1\right)^s}s\left( \left(-1\right)^ssb-y^{-1}\left( \sum_{l=0}^{s-2}\left(-1\right)^{s+l}sb\frac{\left(\ln y\right)^l}{l!} \right.\right.\\&\phantom{ ={}+\sum_{s=2}^\infty\frac{\left(-1\right)^s}s~\left(\vphantom{\sum_s^\infty}\right.\left(-1\right)^ssb-y^{-1} }\left.\left. {}+s\left(1-b\right)\frac{\left(-\ln y\right)^{s-1}}{\left(s-1\right)!} +s\frac{\left(-\ln y\right)^s}{s!} \right) \right)\\ &=b-1+y^{-1}\left(1-b-\ln y\right) +y^{-1}\sum_{s=2}^\infty\left(\frac{\left(\ln y\right)^{s-1}}{\left(s-1\right)!}-\frac{\left(\ln y\right)^s}{s!}\right)\\ &\phantom{=}{}+b\sum_{s=2}^\infty\left(1-y^{-1}\sum_{l=0}^{s-1}\frac{\left(\ln y\right)^l}{l!}\right)\\ &=b-1+y^{-1}\left(1-b-\ln y\right)+y^{-1}\left(\left(\exp\ln y-1\right)-\left(\exp\ln y-\ln y-1\right)\right)\\ &\phantom{=}{}+b\lim_{q\to\infty}\sum_{s=2}^q\left(1-y^{-1}\sum_{l=0}^{s-1}\left(\ln y\right)^l\right)\\ &=b-1+y^{-1}\left(1-b\right) +b\lim_{q\to\infty}\left(q-1-y^{-1}\left(q-1+\sum_{l=1}^{q-1}\frac{q-l}{l!}\left(\ln y\right)^l\right)\right)\\ &=b-1+y^{-1}\left(1-b\right)+b\left(-1+y^{-1}+\ln y\right)\\ &=y^{-1}-1+b\ln y. \end{align*} Compare the result with Equation 13, we have b=0.b=0.

The table of Bq,lB_{q,l} is now k=01234q=111202230033400044⋮⋱\begin{array}{r|ccccc} &k=0&1&2&3&4\\ \hline q=1&1&1\\ 2&0&2&2\\ 3&0&0&3&3\\ 4&0&0&0&4&4\\ \vdots&\ddots \end{array} The table of As,jA_{s,j} is then j=01234s=121224230363400484⋮⋱\begin{array}{r|ccccc} &j=0&1&2&3&4\\ \hline s=1&2&1\\ 2&2&4&2\\ 3&0&3&6&3\\ 4&0&0&4&8&4\\ \vdots&\ddots \end{array} The general formula for As,jA_{s,j} is As,j={s,j∈{s,s−2},2s,j=s−1,0,j∈0 ..s−3.A_{s,j}=\begin{cases} s,&j\in\left\{s,s-2\right\},\\ 2s,&j=s-1,\\ 0,&j\in0\,..s-3. \end{cases} (26)(26)

Therefore, the functions Δgs\Delta g_s are Δgs(y,κ)={−yκ−1ln⁡y(2+ln⁡yκ−1),s=1s(−1)s(s−2)!ysκ−1ln⁡y(ln⁡ysκ−1)s−2(1(ln⁡ysκ−1)2+2s−1ln⁡ysκ−1+1s(s−1)(ln⁡ysκ−1)2),s∈2…∞.\Delta g_s(y,\kappa)=\begin{cases} -y^{\kappa-1}\ln y\left(2+\ln y^{\kappa-1}\right),&s=1\\ \begin{split} &\textstyle{\frac{s\left(-1\right)^s}{\left(s-2\right)!}y^{s\kappa-1}\ln y\left(\ln y^{s\kappa-1}\right)^{s-2} \left(1\vphantom{\left(\ln y^{s\kappa-1}\right)^2}\right.}\\ &\textstyle{\left.\quad{}+\frac2{s-1}\ln y^{s\kappa-1}+\frac1{s\left(s-1\right)}\left(\ln y^{s\kappa-1}\right)^2\right),} \end{split}&s\in2\ldots\infty. \end{cases} Therefore, the functions gqg_q are gq(y,κ)=−yκ−1ln⁡y(2+ln⁡yκ−1)+ln⁡y∑s=2qs(−1)s(s−2)!ysκ−1(ln⁡ysκ−1)s−2(1(−1)ss!gq(y,κ)=+2s−1ln⁡ysκ−1+1s(s−1)(ln⁡ysκ−1)2)\begin{split} &g_q(y,\kappa)=-y^{\kappa-1}\ln y\left(2+\ln y^{\kappa-1}\right) +\ln y\sum_{s=2}^q\frac{s\left(-1\right)^s}{\left(s-2\right)!}y^{s\kappa-1}\left(\ln y^{s\kappa-1}\right)^{s-2} \left(1\vphantom{\frac{\left(-1\right)^s}{s!}}\right.\\ &\phantom{g_q(y,\kappa)=\quad}\left.{}+\frac2{s-1}\ln y^{s\kappa-1}+\frac1{s\left(s-1\right)}\left(\ln y^{s\kappa-1}\right)^2\right) \end{split} (27)(27) (the formula is also applicable to q=1q=1).

Edge cases

Now we have covered almost all cases. The only cases that we have not covered are the cases when κ=1q\kappa=\frac1q, where q∈2…∞q\in2\ldots\infty. The discontinuity in gg at κ=1q\kappa=\frac1q is = g ⁣(y,1q+)−g ⁣(y,1q−)=−Δgq ⁣(y,1q)={−2ln⁡y,q=2,0,q∈3…∞.\begin{split} &\phantom{=~}g\!\left(y,\frac1q^+\right)-g\!\left(y,\frac1q^-\right)\\ &=-\Delta g_q\!\left(y,\frac1q\right)\\ &=\begin{cases} -2\ln y,&q=2,\\ 0,&q\in3\ldots\infty. \end{cases} \end{split} (28)(28) Therefore, for q∈3…∞q\in3\ldots\infty, the function gg has defined limit at κ=1q\kappa=\frac1q, and the value of gg here should just be the limit value. Now, the only problem is at κ=12\kappa=\frac12. We should determine whether the value of gg at κ=12\kappa=\frac12 is its left limit or right limit.

Looking at Equation 12, one may see that the discontinuity at κ=12\kappa=\frac12 is due to the Dirac δ\delta function in the integrand. Therefore, whether gg at κ=12\kappa=\frac12 is g1g_1 or g2g_2 depends on whether the Dirac δ\delta function is within the integrated interval. If it is, then gg at κ=12\kappa=\frac12 is g1g_1; otherwise, it is g2g_2.

The inclusion of the Dirac δ\delta function in the integrated interval corresponds to the inclusion of the highest term in the summation in Equation 7. Because both min⁡(k,n−k−1)\min(k,n-k-1) and min⁡(k−1,n−k−1)\min(k-1,n-k-1) equal n−k−1n-k-1 when n=2kn=2k, the highest term in the summation can be reached, so the Dirac δ\delta function is within the integrated interval. Therefore, gg at κ=12\kappa=\frac12 is g1g_1.

Therefore, we may conclude that for any κ∈(0,1)\kappa\in\left(0,1\right), g(y,κ)=g⌈1κ⌉−1(y,κ).g(y,\kappa)=g_{\lceil\frac1\kappa\rceil-1}(y,\kappa). (29)(29)


Another edge case that is interesting to consider is when κ→0+\kappa\to0^+. However, because the domain of gg does not include κ=0\kappa=0 by definition, so we do not need to consider this case. By some mathematical analysis techniques, one may prove that the limit of gg as κ→0+\kappa\to0^+ is 00.

The solution

Substitute Equation 27 into Equation 29, and we have g(y,κ)=−yκ−1ln⁡y(2+ln⁡yκ−1)+ln⁡y∑s=2⌈1κ⌉−1s(−1)s(s−2)!ysκ−1(ln⁡ysκ−1)s−2(1(−1)ss!g(y,κ)=+2s−1ln⁡ysκ−1+1s(s−1)(ln⁡ysκ−1)2)\begin{split} &g(y,\kappa)=-y^{\kappa-1}\ln y\left(2+\ln y^{\kappa-1}\right) +\ln y\sum_{s=2}^{\lceil\frac1\kappa\rceil-1}\frac{s\left(-1\right)^s}{\left(s-2\right)!}y^{s\kappa-1}\left(\ln y^{s\kappa-1}\right)^{s-2} \left(1\vphantom{\frac{\left(-1\right)^s}{s!}}\right.\\ &\phantom{g(y,\kappa)=\quad}\left.{}+\frac2{s-1}\ln y^{s\kappa-1}+\frac1{s\left(s-1\right)}\left(\ln y^{s\kappa-1}\right)^2\right) \end{split} Substitute the result into Equation 11 and then Equation 10, and also consider Equation 8 and 9, and we have f(y,κ)={δ(κ),y=0,κ∈[0,1],0,y∈(0,1],κ=0,yδ(κ−1)−yκln⁡y(2+ln⁡yκ−1)+ln⁡y∑s=2⌈1κ⌉−1s(−1)s(s−2)!ysκ(ln⁡ysκ−1)s−2(1(ln⁡ysκ−1)2+2s−1ln⁡ysκ−1+1s(s−1)(ln⁡ysκ−1)2),y∈(0,1],κ∈(0,1].f(y,\kappa)=\begin{cases} \delta(\kappa),&y=0,\kappa\in[0,1],\\ 0,&y\in(0,1],\kappa=0,\\ \begin{split} &\textstyle{y\delta(\kappa-1)-y^\kappa\ln y\left(2+\ln y^{\kappa-1}\right)}\\ &\quad\textstyle{ {}+\ln y\sum_{s=2}^{\lceil\frac1\kappa\rceil-1}\frac{s\left(-1\right)^s}{\left(s-2\right)!} y^{s\kappa}\left(\ln y^{s\kappa-1}\right)^{s-2} \left(1\vphantom{\left(\ln y^{s\kappa-1}\right)^2}\right.}\\ &\qquad\textstyle{\left.{}+\frac2{s-1}\ln y^{s\kappa-1}+\frac1{s\left(s-1\right)}\left(\ln y^{s\kappa-1}\right)^2\right)}, \end{split}&y\in(0,1],\kappa\in(0,1]. \end{cases} (30)(30)

Plots of the probability density functions

Here are plots of the function f(y,κ)f(y,\kappa) whose expression is given by Equation 30:

Probability distribution of  when

We can compare it with a plot of the distributions when nn is finite (say, 100100), and we may see that they are very close:

Probability distribution of  when  and when  compared

We have not investigated the asymptotic behavior of the error if we approximate the distribution with finite nn by the distribution with infinite nn, but we may expect that the error is small enough for applicational uses when nn is a usual note count in a rhythm game chart (usually at least 500500).

Moments

It may be interesting to calculate the moments of the distribution.

We need to evaluate μν ⁣(y)≔∫01κνf ⁣(y,κ)dκ.\mu_\nu\!\left(y\right)\coloneqq\int_0^1\kappa^\nu f\!\left(y,\kappa\right)\mathrm d\kappa. First, calculate = ∫01sκν(ln⁡ysκ−1)jysκ−1ln⁡y dκ=∫y−11(log⁡yw+1s)ν(ln⁡w)j dw=1sν+1∑p=0ν(νp)1(ln⁡y)p∫y−11(ln⁡w)j+p dw=1sν+1∑p=0ν(νp)(−1)j+p(j+p)!(ln⁡y)p(1−y−1∑l=0j+p(ln⁡y)ll!).\begin{align*} &\phantom{=~}\int_0^{\frac1s}\kappa^\nu\left(\ln y^{s\kappa-1}\right)^jy^{s\kappa-1}\ln y\,\mathrm d\kappa\\ &=\int_{y^{-1}}^1\left(\frac{\log_yw+1}{s}\right)^\nu\left(\ln w\right)^j\,\mathrm dw\\ &=\frac1{s^{\nu+1}}\sum_{p=0}^\nu\binom\nu p\frac1{\left(\ln y\right)^p} \int_{y^{-1}}^1\left(\ln w\right)^{j+p}\,\mathrm dw\\ &=\frac{1}{s^{\nu+1}}\sum_{p=0}^\nu\binom\nu p\frac{\left(-1\right)^{j+p}\left(j+p\right)!}{\left(\ln y\right)^p} \left(1-y^{-1}\sum_{l=0}^{j+p}\frac{\left(\ln y\right)^l}{l!}\right). \end{align*} = ∫01sκνΔgs ⁣(y,κ)dκ=∑j=0s(−1)sAs,jj!1sν+1∑p=0ν(νp)(−1)j+p(j+p)!(ln⁡y)p(1−y−1∑l=0j+p(ln⁡y)ll!)=(−1)ssν+1∑p=0ν(νp)(−1)p(ln⁡y)p(∑j=0s(j+p)!(−1)jj!As,j−y−1∑j=0s(j+p)!(−1)jj!As,j∑l=0j+p(ln⁡y)ll!).\begin{align*} &\phantom{=~}\int_0^{\frac1s}\kappa^\nu\Delta g_s\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{j=0}^s\left(-1\right)^s\frac{A_{s,j}}{j!} \frac{1}{s^{\nu+1}}\sum_{p=0}^\nu\binom\nu p\frac{\left(-1\right)^{j+p}\left(j+p\right)!}{\left(\ln y\right)^p} \left(1-y^{-1}\sum_{l=0}^{j+p}\frac{\left(\ln y\right)^l}{l!}\right)\\ &=\frac{\left(-1\right)^s}{s^{\nu+1}}\sum_{p=0}^\nu\binom\nu p\frac{\left(-1\right)^p}{\left(\ln y\right)^p} \left( \sum_{j=0}^s\frac{\left(j+p\right)!\left(-1\right)^j}{j!}A_{s,j} \right.\\&\qquad\qquad\qquad\qquad\qquad\qquad\left.{} -y^{-1}\sum_{j=0}^s\frac{\left(j+p\right)!\left(-1\right)^j}{j!}A_{s,j}\sum_{l=0}^{j+p}\frac{\left(\ln y\right)^l}{l!} \right). \end{align*}

Define Bs,l,p≔∑j=max⁡(0,l−p)s(j+p)!(−1)jj!As,j,B_{s,l,p}\coloneqq\sum_{j=\max(0,l-p)}^s\frac{\left(j+p\right)!\left(-1\right)^j}{j!}A_{s,j}, (31)(31) Dν,p,l≔∑s=max⁡(1,l−p)∞(−1)ssν+1Bs,l,p.D_{\nu,p,l}\coloneqq\sum_{s=\max(1,l-p)}^\infty\frac{\left(-1\right)^s}{s^{\nu+1}}B_{s,l,p}. (32)(32) Then, = ∫01sκνΔgs ⁣(y,κ)dκ=(−1)ssν+1∑p=0ν(νp)(−1)p(ln⁡y)p(Bs,0,p−y−1∑l=0s+pBs,l,p(ln⁡y)ll!).\begin{align*} &\phantom{=~}\int_0^{\frac1s}\kappa^\nu\Delta g_s\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\frac{\left(-1\right)^s}{s^{\nu+1}}\sum_{p=0}^\nu\binom\nu p\frac{\left(-1\right)^p}{\left(\ln y\right)^p} \left(B_{s,0,p}-y^{-1}\sum_{l=0}^{s+p}B_{s,l,p}\frac{\left(\ln y\right)^l}{l!}\right). \end{align*} = ∫01κνg ⁣(y,κ)dκ=∑s=1∞(−1)ssν+1∑p=0ν(νp)(−1)p(ln⁡y)p(Bs,0,p−y−1∑l=0s+pBs,l,p(ln⁡y)ll!)=∑p=0ν(νp)(−1)p(ln⁡y)p(Dν,p,0−y−1∑l=0∞Dν,p,l(ln⁡y)ll!).\begin{align*} &\phantom{=~}\int_0^1\kappa^\nu g\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\sum_{s=1}^\infty \frac{\left(-1\right)^s}{s^{\nu+1}}\sum_{p=0}^\nu\binom\nu p\frac{\left(-1\right)^p}{\left(\ln y\right)^p} \left(B_{s,0,p}-y^{-1}\sum_{l=0}^{s+p}B_{s,l,p}\frac{\left(\ln y\right)^l}{l!}\right)\\ &=\sum_{p=0}^\nu\binom\nu p\frac{\left(-1\right)^p}{\left(\ln y\right)^p} \left(D_{\nu,p,0}-y^{-1}\sum_{l=0}^\infty D_{\nu,p,l}\frac{\left(\ln y\right)^l}{l!}\right). \end{align*}

Now, the only problem is how to get Dν,p,lD_{\nu,p,l}. Substitute Equation 26 into Equation 31, and after some calculations, we can get the general formula of Bs,l,pB_{s,l,p}: Bs,l,p=(−1)smax⁡(l,s+p−2)!(s−1)!⋅{p(p−1),l∈0 ..s+p−2,p−s,l=s+p−1,1,l=s+p.B_{s,l,p}=\frac{\left(-1\right)^s\max(l,s+p-2)!}{\left(s-1\right)!}\cdot\begin{cases} p\left(p-1\right),&l\in0\,..s+p-2,\\ p-s,&l=s+p-1,\\ 1,&l=s+p. \end{cases} Substitute it into Equation 32, and notice the edge cases, we can get Dν,p,l={p(p−1)∑s=1∞(s+p−2)!sν+1(s−1)!,l∈0 ..p−1,p!(p−1)+p(p−1)∑s=2∞(s+p−2)!sν+1(s−1)!,l=p,l!(l−p)ν+1(l−p−1)!−l!(2p−l−1)(l−p+1)ν+1(l−p)!+p(p−1)∑s=l−p+2∞(s+p−2)!sν+1(s−1)!,l∈p+1…∞={p(p−1)((p−1)!+Sν,p),l∈0 ..p,l!(l−p)ν+1(l−p−1)!−l!(2p−l−1)(l−p+1)ν+1(l−p)!+p(p−1)(Sν,p−∑s=2l−p+2(s+p−2)!sν+1(s−1)!),l∈p+1…∞,\begin{align*} D_{\nu,p,l}&=\begin{cases} p\left(p-1\right)\sum_{s=1}^\infty\frac{\left(s+p-2\right)!}{s^{\nu+1}\left(s-1\right)!}, &l\in0\,..p-1,\\ p!\left(p-1\right)+p\left(p-1\right) \sum_{s=2}^\infty\frac{\left(s+p-2\right)!}{s^{\nu+1}\left(s-1\right)!}, &l=p,\\ \begin{split} &\textstyle{\frac{l!}{\left(l-p\right)^{\nu+1}\left(l-p-1\right)!} -\frac{l!\left(2p-l-1\right)}{\left(l-p+1\right)^{\nu+1}\left(l-p\right)!}}\\ &\qquad\textstyle{ {}+p\left(p-1\right)\sum_{s=l-p+2}^\infty\frac{\left(s+p-2\right)!}{s^{\nu+1}\left(s-1\right)!},} \end{split} &l\in p+1\ldots\infty \end{cases}\\ &=\begin{cases} p\left(p-1\right)\left(\left(p-1\right)!+S_{\nu,p}\right),& l\in0\,..p,\\ \begin{split} &\textstyle{\frac{l!}{\left(l-p\right)^{\nu+1}\left(l-p-1\right)!} -\frac{l!\left(2p-l-1\right)}{\left(l-p+1\right)^{\nu+1}\left(l-p\right)!}}\\ &\qquad\textstyle{ {}+p\left(p-1\right)\left(S_{\nu,p}-\sum_{s=2}^{l-p+2}\frac{\left(s+p-2\right)!}{s^{\nu+1}\left(s-1\right)!}\right),} \end{split} &l\in p+1\ldots\infty, \end{cases} \end{align*} where the infinite sum Sν,p≔∑s=2∞(s+p−2)!sν+1(s−1)!.S_{\nu,p}\coloneqq\sum_{s=2}^\infty\frac{\left(s+p-2\right)!}{s^{\nu+1}\left(s-1\right)!}. There is no closed form for Sν,pS_{\nu,p}, but we may express it in terms of Stirling numbers of the first kind and the Riemann ζ\zeta function. For p∈1 ..νp\in1\,..\nu, we have Sν,p=−(p−1)!+∑s=1∞s(s+1)⋯(s+p−2)sν+1=−(p−1)!+∑s=1∞1sν+1∑λ=0p−1[p−1λ]sλ=−(p−1)!+∑λ=0p−1[p−1λ]ζ ⁣(ν−λ+1),\begin{align*} S_{\nu,p}&=-\left(p-1\right)!+\sum_{s=1}^\infty\frac{s\left(s+1\right)\cdots\left(s+p-2\right)}{s^{\nu+1}}\\ &=-\left(p-1\right)!+\sum_{s=1}^\infty\frac1{s^{\nu+1}}\sum_{\lambda=0}^{p-1}\begin{bmatrix}p-1\\\lambda\end{bmatrix}s^\lambda\\ &=-\left(p-1\right)!+\sum_{\lambda=0}^{p-1}\begin{bmatrix}p-1\\\lambda\end{bmatrix}\zeta\!\left(\nu-\lambda+1\right), \end{align*} where [⋅⋅]\begin{bmatrix}\cdot\\\cdot\end{bmatrix} denotes (unsigned) Stirling numbers of the first kind. For p=0p=0, we have Sν,0=∑s=2∞1sν+1(s−1)=∑s=2∞1sν+1(s−1)−∑s=2∞1s(s−1)+1=1−∑s=2∞1sν+1sν−1s−1=ν+1−∑s=1∞1sν+1∑λ=0ν−1sλ=ν+1−∑λ=0ν−1ζ ⁣(ν−λ+1).\begin{align*} S_{\nu,0} &=\sum_{s=2}^\infty\frac1{s^{\nu+1}\left(s-1\right)}\\ &=\sum_{s=2}^\infty\frac1{s^{\nu+1}\left(s-1\right)}-\sum_{s=2}^\infty\frac1{s\left(s-1\right)}+1\\ &=1-\sum_{s=2}^\infty\frac1{s^{\nu+1}}\frac{s^\nu-1}{s-1}\\ &=\nu+1-\sum_{s=1}^\infty\frac1{s^{\nu+1}}\sum_{\lambda=0}^{\nu-1}s^\lambda\\ &=\nu+1-\sum_{\lambda=0}^{\nu-1}\zeta\!\left(\nu-\lambda+1\right). \end{align*}

Then, the following steps will be extremely tedious, and I doubt there will be a closed form for our final result, so I will not continue to find the general formula for the moments.


However, we may obtain the first moment (mean) analytically. We have D1,0,l={−1,l=0,1l−1l+1,l∈1…∞,D1,1,l={0,l=0,1,1l+1l−1,l∈2…∞.D_{1,0,l}=\begin{cases}-1,&l=0,\\\frac1l-\frac1{l+1},&l\in1\ldots\infty,\end{cases} \quad D_{1,1,l}=\begin{cases}0,&l=0,1,\\\frac1l+\frac1{l-1},&l\in2\ldots\infty.\end{cases} μ1 ⁣(y)≔∫01κf ⁣(y,κ)dκ=y+y∫01κg ⁣(y,κ)dκ=li⁡y−ln⁡(−ln⁡y)−γln⁡y,\begin{align*} \mu_1\!\left(y\right) &\coloneqq\int_0^1\kappa f\!\left(y,\kappa\right)\mathrm d\kappa\\ &=y+y\int_0^1\kappa g\!\left(y,\kappa\right)\mathrm d\kappa\\ &=\frac{\operatorname{li}y-\ln(-\ln y)-\gamma}{\ln y}, \end{align*} where li⁡\operatorname{li} is the logarithmic integral function, and γ\gamma is the Euler–Mascheroni constant. The function seems undefined when y=0y=0 or y=1y=1, but it has limits at these points: μ1 ⁣(y→0+)=0,μ1 ⁣(y→1−)=1,\mu_1\!\left(y\to0^+\right)=0,\quad\mu_1\!\left(y\to1^-\right)=1, which is intuitive. (This function tends to 00 very slowly when y→0+y\to0^+, so slowly that I almost did not believe that when I did the numerical calculation first.)

The plot:

The mean value of  vs. 

We should also be able to find other statistical quantities like the median, the mode, the variance, etc., but they seem do not have closed forms.

Some interesting observations

The probability distribution of κ\kappa seems to tend to be a uniform distribution plus a Dirac δ\delta distribution when yy is very close to 11. This phenomenon is very visible if we look at the plot of f(y=0.9,κ)f(y=0.9,\kappa).

In other words, the distribution seems like f(y≈1,κ)≈(1−y)U ⁣(12,1)+yδ(κ−1),f(y\approx1,\kappa)\approx \left(1-y\right)U\!\left(\frac12,1\right)+y\delta(\kappa-1), where U(a,b)U(a,b) denotes the uniform distribution on the interval [a,b][a,b].

This can be justified by expanding f(y,κ)f(y,\kappa) in Taylor series of 1−y1-y and retaining the first-order terms only. Note that ya(ln⁡y)b=(y−1)b(1+(b2−a)(1−y)+⋯ ),y^a\left(\ln y\right)^b= \left(y-1\right)^b\left(1+\left(\frac b2-a\right)\left(1-y\right)+\cdots\right), so the only case where the Taylor series has a non-zero first-order term is when b=1b=1 or b=0b=0. In Equation 30, we can see that the power on ln⁡y\ln y is at least one for each term (because of the general ln⁡y\ln y factor in front), so only the terms with no ln⁡y\ln y factors but the general one will have a first-order term. In this case, the first order term is proportional to y−1y-1, and the proportional coefficient is just the coefficient in the front of the term in ff, which is independent of κ\kappa because κ\kappa only appears in the power index of yy.

Therefore, we may see that only q=1q=1 and q=2q=2 terms have a non-zero first-order term, and they are respectivey −2(y−1)-2\left(y-1\right) and 2(y−1)2\left(y-1\right). This means that when yy is very close to 11, f(y≈1,κ)≈{2(1−y),κ∈(12,1),0,κ∈(0,12).f(y\approx 1,\kappa)\approx\begin{cases} 2\left(1-y\right),&\kappa\in\left(\frac12,1\right),\\ 0,&\kappa\in\left(0,\frac12\right). \end{cases} This is exactly the uniform distribution.

There is an intuitive way to explain the appearance of the uniform distribution. When yy is very close to 11, the probability of getting one combo break (1−Y1-Y) is already very small, so it is very unlikely that there are two or more combo breaks. Assuming there is only one combo break and it may appear anywhere with equal probability. The combo break will cut the string of notes into two pieces, and the length of the larger piece is the max combo, which is uniformly distributed between half note count and full note count.


Every rhythm game player knows: never celebrate too early. You never know whether you will miss near the end. It is then interesting to know what is the probability of getting almost a full combo, i.e. what is the probability of getting κ\kappa very close to 11.

If we find the limit of f(y,κ)f(y,\kappa) as κ→1−\kappa\to1^-, it is f ⁣(y,κ→1−)=−2yln⁡y.f\!\left(y,\kappa\to1^-\right)=-2y\ln y. There is a peak of this probability density at y=e−1y=\mathrm e^{-1}. Therefore, when y=e−1y=\mathrm e^{-1}, the probability of getting κ\kappa very close to 11 is the largest.

When does y=e−1y=\mathrm e^{-1}, then? Because y=Yn=(1−nbn)n,y=Y^n=\left(1-\frac{n_\mathrm b}n\right)^n, where nbn_\mathrm b is the average number of combo breaks, then it tends to e−nb\mathrm e^{-n_\mathrm b} when n→∞n\to\infty. Therefore, the probability of getting almost a full combo is the highest when your average number of combo breaks is exactly one.


From the plot, it seems that the probability of getting κ\kappa a little bit higher than 12\frac12 is always higher than the probability of getting κ\kappa a little bit lower than 12\frac12. According to Equation 28, the jump in f(y,κ)f(y,\kappa) at κ=12\kappa=\frac12 is f ⁣(y,κ→12+)−f ⁣(y,κ→12−)=−2yln⁡y.f\!\left(y,\kappa\to\frac12^+\right)-f\!\left(y,\kappa\to\frac12^-\right)=-2y\ln y. Interestingly, this coincides with f ⁣(y,κ→1−)f\!\left(y,\kappa\to1^-\right).


Define y0(κ)≔arg maxy∈[0,1] f(y,κ),y_0(\kappa)\coloneqq\mathop{\mathrm{arg\,max}}_{y\in[0,1]}\,f(y,\kappa), and then it seems that y0:[0,1]→[0,1]y_0:[0,1]\to[0,1] is injective but not surjective. It is strictly increasing, and there is a jump at κ=12\kappa=\frac12 and at κ=1\kappa=1.

It has an elementary expression on [12,1)\left[\frac12,1\right): y0 ⁣(κ∈[12,1))=exp⁡−2κ+1+2κ2−2κ+1κ(κ−1).y_0\!\left(\kappa\in\left[\frac12,1\right)\right) =\exp\frac{-2\kappa+1+\sqrt{2\kappa^2-2\kappa+1}}{\kappa\left(\kappa-1\right)}.

Some applications

In Phigros, one should combo at least 60%60\% of the notes to get a white V (white V) rank. If on average I have one combo break in a chart, which has 13001300 notes, what is the probability of comboing at least 60%60\% of the notes in the chart?

Solution. The success rate is Y=1300−11300,y=Y1300≈e−1.Y=\frac{1300-1}{1300},\quad y=Y^{1300}\approx\mathrm e^{-1}. The probability of comboing more than 6060\\% of the notes is = ∫60%1f ⁣(y,κ)dκ=y+∫0.61−yκln⁡y(2+ln⁡yκ−1)dκ≈e−1+∫0.61e−κ(3−κ)dκ=75e−35≈0.768.\begin{align*} &\phantom{=~}\int_{60\%}^1f\!\left(y,\kappa\right)\mathrm d\kappa\\ &=y+\int_{0.6}^1-y^\kappa\ln y\left(2+\ln y^{\kappa-1}\right)\mathrm d\kappa\\ &\approx\mathrm e^{-1}+\int_{0.6}^1\mathrm e^{-\kappa}\left(3-\kappa\right)\mathrm d\kappa\\ &=\frac75\mathrm e^{-\frac35}\\ &\approx0.768. \end{align*}


Oh, my god! It is hard to come up with application problems. I hope readers find out the applications themselves.