The original question is here (The Board Football Problem (Probability)) and part II is here(The Board Football Problem (Part II)). I was told to segment the question in order to increase the chances of it being answered.
As kids we used to play the game "Board Football". I'm not particularly sure of the popularity of the game but the rules were pretty simple to follow. The basic objective of the game was to score as many goals as possible and you played it with nothing more than dice and a notebook.
It's a two player game, which goes on as long as the players don't bored. ie there isn't any way to win/finish the game.Both you (say A) and the opponent (say B) have a dice. Before starting the game you decide the minimum number of "passes" that need to be completed before a goal is scored(say p).
Let's assume that A starts first.A rolls the dice and gets a certain number, say 5. This means that A completes 5 passes in that round and has to complete p-5 more passes in order to progress to the "Final Stage". However, before the round is completed the opponent(B) has the opportunity to "intercept" you and obtain possession of the ball.
The way to do this is pretty simple. B rolls his dice and if he gets the same number as A, then B obtains possession of the ball and the round is completed. In the next round, B will start off from 0 passes.
We weren't particularly finicky back then so even if we rolled a number to exceed the number of passes made, the player in possession was still allowed to progress to the "Final Stage". It didn't matter in what order the passes in each round were made as long as the required passes was obtained.
Once the required passes were made and the player reached the "Final Stage", they had to flip a coin. If it turns up as heads,the player in possession scores and if it turns up as tails, he doesn't score. In both the cases, the possession is reverted back to B for the next round. We considered the coin-flipping part to take place in the same round as the round in which the player in possession completed the required passes.
For example, let a 6 sided dice be used and the required number of passes be 9. (That is n=6,p=9) Let us say that A has the ball in the starting.
In the first round, A rolls 4 and B rolls 5. A completes 4 passes within that round.A now has to complete at least 5 more passes in order to get an opportunity to shoot.
In the second round, A rolls 3 and B rolls 3. Possession is overturned and now B controls the ball.B now has 0 passes completed.
In the third round, B rolls 5 and A rolls 2.B has completed 5 passes in this round and needs to complete at least 4 more passes in order to get an opportunity to shoot.
In the fourth round, B rolls 6 and A rolls 4. Although B has overshot the number of passes, B is permitted to flip the coin. B tosses the coin and gets heads. B scores a goal.B leads by one goal to nil. Now A starts off with possession of the ball.
So my question is as follows:If we use an n-sided dice and the number of required passes is p, the what is the probability that the player in possession of the ball(starting off from 0 passes completed) scores a goal without losing possession of the ball. Is it possible to calculate the value for this in a non-computational /formulaic manner?
(BONUS QUESTION) Form a graph with the Y-Axis being the probability that the player in possession scores without losing possession of the ball and the number of required passes be the X-Axis. (Let a 10 sided dice be used for this case.)
If possible, determine how different the graph would look [in a visual sense, no need for deeper calculations] if the number of sides on the dice is made to vary. Is the graph just shifted along the x-axis or do the peaks become more gradual?).
PS This isn't a homework question/problem so please don't close it because an answer is being requested. It is an original question. Consider this to be a challenge of sorts.
Disclaimers
I cannot provide a full explanation of all techniques used in this answer. If there is something I use where searching for the key terms does not provide enough clarity for you to find a good resource explaining the idea, just ask.
I don't have a background in things like numerical approximations, and am not very familiar with theorems about things like Markov models or signal processing, so there may be better approaches to many of the calculations I do below.
At the time of writing, I have the calculations/references for most of the Appendix ready, but have not written up all of the exposition just yet.
Setup
Without loss of generality, call the player in possession "A", and the other player "B". They will use an $n$-sided die. Let $P(p)$ be the probability that A scores before any change in possession, when A needs $p$ more passes to have the opportunity to shoot. The question essentially asks for a nice way to calculate $P(p)$ for arbitrary $n$, as well as graphical insight.
Manual Calculation
Recursive Formula
Note that $P(0)=\dfrac{1}{2}$ since that's the chance of winning the coin flip.
When the $n$-sided die is rolled by the players, there are various possibilities of varying probabilities: One possibility is that of an interception, which happens with probability $\dfrac{1}{n}$ when B rolls the right number. And there are $n$ possibilities are of the form "A makes $i$ passes" for some $i\in\left\{ 1,\ldots,n\right\}$. In each of those cases, there is a $\dfrac1n$ chance that A rolls an $i$, and a $1-\dfrac1n$ chance that B does not roll an $i$. Thus, each one of those occurs with probability $\dfrac{1}{n}\left(1-\dfrac{1}{n}\right)=\dfrac{n-1}{n^{2}}$ by the product rule for probability.
When A makes $i$ passes: if $i\le p$ then there are $p-i$ passes left to go, and if $i\ge p$, there are $0$ passes left to go, only the coin flip. Then the sum rule of probability tells us that, for $p>0$, $P(p)=\dfrac{1}{n}*0+{\displaystyle \sum_{i=1}^{n}}\dfrac{n-1}{n^{2}}P\left(\max\left\{ p-i,0\right\}\right)=\dfrac{n-1}{n^{2}}{\displaystyle \sum_{i=1}^{n}}P\left(\max\left\{ p-i,0\right\}\right)$. The first term of $0$ was because there is no chance of winning before a change in possession when there is an interception.
Computer Code
This formula is relatively straightforward to plug into a computer. For example, if we wanted to calculate $P(7)$ for $n=10$, the following Python code outputs
0.41614499445795. Try it online!And the following translation into Wolfram Language code outputs the exact answer of $\dfrac{8322899889159}{20000000000000}$. Try it online!
However, I don't think this qualifies as a "non-computational /formulaic manner", as requested in the question.
Mathematical Formula
General Information
The "$\max$" part makes this challenging to fit into existing theory about recurrences. So, for convenience, define $P\left(1-i\right)=\dfrac{1}{2}$ for $1\le i\le n$ . With those as initial values (sometimes called "initial conditions"), we can simply write $P(p)=\dfrac{n-1}{n^{2}}{\displaystyle \sum_{i=1}^{n}}P\left(p-i\right)=\dfrac{n-1}{n^{2}}{\displaystyle \sum_{j=p-n}^{p-1}}P(j)$ for $p>0$. This is now a "linear homogeneous recurrence relation of order $n$ with constant coefficients", which is a well-studied situation.
Our relation has "characteristic equation" $x^{n}=\dfrac{n-1}{n^{2}}{\displaystyle \sum_{j=0}^{n-1}}x^{j}$. It turns out that this equation always has $n$ distinct complex roots $r_{1},\ldots,r_{n}$ (see the Appendix below for a proof). That means that $P(p)$ has a particularly nice form: $P(p)={\displaystyle \sum_{i=1}^{n}}C_{i}r_{i}^{p}$ for some constants $C_{i}$. In fact, the symmetry allows us to write the constants $C_{i}$ explicitly in terms of the roots $C_{i}=\dfrac{1}{2}r_{i}^{n-1}{\displaystyle \prod_{j\ne i}}\dfrac{1-r_{j}}{r_{i}-r_{j}}$. And Using the theory of symmetric polynomials, we can express this in terms of just $r_i$ and $n$, if desired. See the Appendix for all details.
Small $n$
For $n=1$, there is always an interception if $p>0$, so $P(p)=0$ for $p>0$, and we might still take $P(0)=\dfrac{1}{2}$. This is actually consistent with the approach and formula outlined above. The characteristic equation becomes $x^{1}=\dfrac{1-1}{1^2}{\displaystyle \sum_{j=0}^{1-1}}x^{j}=\dfrac01*x^{0}=0$, with root $r_{1}=0$. Then we have $P(p)=\dfrac{1}{2}0^{1-1}[\text{empty product}]0^{p}$. Assuming the combinatorics conventions that $0^{0}=1$ in discrete situations and the empty product is $1$, this reduces to $P(p)=\dfrac{1}{2}*\begin{cases}1 & \text{ if }p=0\\0 & \text{ otherwise}\end{cases}$, as expected.
For $n=2$, the characteristic equation becomes $x^{2}=\dfrac{2-1}{2^2}{\displaystyle \sum_{j=0}^{2-1}}x^{j}=\dfrac{1}{4}\left(x+1\right)$. The roots of this quadratic are $r_{1,2}=\dfrac{1\pm\sqrt{17}}{8}$. Then $C_{1,2}=\dfrac{1}{2}\left(\dfrac{1\pm\sqrt{17}}{8}\right)\dfrac{1-\dfrac{1\mp\sqrt{17}}{8}}{\dfrac{1\pm\sqrt{17}}{8}-\dfrac{1\mp\sqrt{17}}{8}}=\dfrac{17\pm3\sqrt{17}}{68}$. So we have $P(p)=\dfrac{17+3\sqrt{17}}{68}\left(\dfrac{1+\sqrt{17}}{8}\right)^{p}+\dfrac{17-3\sqrt{17}}{68}\left(\dfrac{1-\sqrt{17}}{8}\right)^{p}$. For $p$ from $0$ to $9$, this produces expressions which simplify down to $\frac{1}{2},\frac{1}{4},\frac{3}{16},\frac{7}{64},\frac{19}{256},\frac{47}{1024},\frac{123}{4096},\frac{311}{16384},\frac{803}{65536},\frac{2047}{262144}$.
Note that for large $p$, $\left(1+\sqrt{17}\right)^{p}\approx5^{p}$ is much larger in size than $\left(1-\sqrt{17}\right)^{p}\approx\left(-3\right)^{p}$, so we can approximate by dropping the term with the minus sign: $P(p)\approx\dfrac{17+3\sqrt{17}}{68}\left(\dfrac{1+\sqrt{17}}{8}\right)^{p}$ for large $p$.
For $n=3$, the roots $r_{1,2,3}$ of the characteristic equation $x^{3}=\dfrac{8}{81}\left(1+r+r^{2}\right)$ are given (via the cubic formula) by $r_{k}=\dfrac{1}{27}\left(2+\omega^{k}\sqrt[3]{2357+81\sqrt{817}}+\omega^{-k}\sqrt[3]{2357-81\sqrt{817}}\right)$ where $\omega=\dfrac{-1+i\sqrt{3}}{2}$ (so that $\omega^{-1}=\omega^{2}=\overline{\omega}=\dfrac{-1-i\sqrt{3}}{2}$). We have $c_{1}=\dfrac{r_{1}^{2}\left(1-r_{2}\right)\left(1-r_{3}\right)}{2\left(r_{1}-r_{2}\right)\left(r_{1}-r_{3}\right)}$ and similarly for the other constants. Then $P(p)=c_{1}r_{1}^{p}+c_{2}r_{2}^{p}+c_{3}r_{3}^{p}$. For $p$ from $0$ to $7$, this produces expressions which simplify down to $\frac{1}{2},\frac{1}{3},\frac{8}{27},\frac{61}{243},\frac{428}{2187},\frac{3250}{19683},\frac{24086}{177147},\frac{176008}{1594323}$.
Note that $\left|27r_{1}\right|\approx22$ but (using complex absolute value), $\left|27r_{2,3}\right|\approx14$. So for large $p$, we may write $P(p)\approx c_{1}r_{1}^{p}$. This approximate formula can be written using real numbers alone. Let $s_{1,2}=\sqrt[3]{2357\pm81\sqrt{817}}$, so that $r_{1}=\dfrac{1}{27}\left(2+s_{1}+s_{2}\right)$. Then $c_{1}$ simplifies to $\dfrac{74358+2438\left(s_{1}+s_{2}\right)+243\left(s_{1}^{2}+s_{2}^{2}\right)+25\left(s_{1}s_{2}^{2}+s_{1}^{2}s_{2}\right)}{1458\left(58+s_{1}^{2}+s_{2}^{2}\right)}$.
Higher $n$
For $n=4$, we can do something similar to $n=3$, but the expressions become worse since we're dealing with quartics. For $n=5,6$, and probably all higher $n$, we cannot write $r_{i}$ in terms of radicals. (Abel's Theorem gives us reason to expect that.) We can still use approximations of the roots and evaluate $P(p)={\displaystyle \sum_{i=1}^{n}}c_{i}r_{i}^{p}$. These powers of complex numbers are not so bad if we use something like De Moivre's formula.
If we want to avoid taking so many powers, we can always use a single term approximation as we used for low values of $n$. By Descartes' rule of signs, there is always a single positive real root of the characteristic equation. By a result of Cauchy (discussed in the Appendix), this will be the root of largest absolute value of a polynomial like the characteristic polynomials we have here. So we can always approximate $P(p)\approx c_{1}r_{1}^{p}$ by choosing $r_{1}$ to be the positive real root of the characteristic equation $x^{n}=\dfrac{n-1}{n^{2}}\left(1+x+\cdots+x^{n-1}\right)=\dfrac{\left(n-1\right)\left(1-x^{n}\right)}{n^{2}\left(1-x\right)}$. The last expression comes from the sum of a geometric progression.
Graphs
Graph for $n=10$
When $n=10$, we can first look at low values of $p$ ($0\le p\le 15$, say), for which $P(p)$ is a little erratic: $\frac{1}{2},\frac{9}{20},\frac{891}{2000},\frac{88119}{200000},\frac{8704971}{20000000},\frac{858841839}{2000000000},\frac{84613760451}{200000000000},\frac{8322899889159}{20000000000000 },\frac{817196087918331}{2000000000000000},\frac{80074373583098079}{200000000000000000},\frac{7828106720557690611}{20000000000000000000},\frac{763263632540788276599}{2000000000000000000000} ,\frac{75095735946945922149291}{200000000000000000000000},\frac{7383535218217105514272719}{20000000000000000000000000},\frac{725498238785664501055726371}{2000000000000000000000000000},\frac {71244834127637430615074174439}{200000000000000000000000000000}$
But things smooth out for higher values of $p$. For $15\le p\le100$, we have a graph like:
By the analysis worked out earlier, for large $p$, this can be approximated by an exponential function $E(p)$. For $n=10$, $E(p)$ happens to be approximately $0.472463*.981288^p$. This fits extremely well for $p\ge 15$ (the approximation is in dashed orange):
We can understand the fit better by looking at the signed relative error $\dfrac{E(p)-P(p)}{P(p)}$, plotted for $15\le p\le100$ here:
For the record, the approximation is not very good for small $p$:
Graphs for arbitrary $n$
We can plot various $n$ ($1\le n\le 20$) together for the various $p$ ranges.
$0\le p \le 15$:
$0\le p \le 100$:
$0\le p \le 500$:
One other way to look at things collectively is to examine the approximations $E(p)=C(n)*r(n)^p$. Plotting $r(n)$ (approaching $1$) and $C(n)$ (approaching $\frac12$) for $2\le n\le50$ yields:
Appendix
Distinct roots
The characteristic polynomial is $x^{n}-\dfrac{n-1}{n^{2}}\left(x^{n-1}+\cdots+1\right)$. We would like to show that this has no repeated roots. It doesn't take much extra effort to show a slight generalization: For any positive constant $t$, $\chi\left(x\right)=x^{n}-t\left(x^{n-1}+\cdots+1\right)$ has no repeated roots.
Let $r$ be a multiple root of $\left(x-1\right)\chi\left(x\right)=x^{n+1}-\left(t+1\right)x^{n}+t$ (the simplification comes from the telescoping sum caused by the finite geometric series), so that $\left(x-r\right)^{2}$ is a factor of $\left(x-1\right)\chi\left(x\right)$. Then either $r$ is a multiple root of $\chi\left(x\right)$ or $r$ is a single root of $\chi\left(x\right)$ and $r=1$.
Since $\left(x-r\right)^{2}$ is a factor of $\left(x-1\right)\chi\left(x\right)$, the product rule tells us that $r$ must be a zero of the derivative $\dfrac{\mathrm{d}}{\mathrm{d}x}\left(x-1\right)\chi\left(x\right)=\left(n+1\right)x^{n}-\left(t+1\right)nx^{n-1}=x^{n-1}\left(\left(n+1\right)x-\left(t+1\right)\right)$. This forces $r=0$ or $r=\dfrac{n\left(t+1\right)}{n+1}$. Since $t>0$, $0$ cannot be a root of $\left(x-1\right)\chi\left(x\right)=x^{n+1}-\left(t+1\right)x^{n}+t$, so we must have $r=\dfrac{n\left(t+1\right)}{n+1}$.
This value of $r$ must also be a root of the original $\left(x-1\right)\chi\left(x\right)$, so that $\left(\frac{n(t+1)}{n+1}\right)^{n+1}-(t+1)\left(\frac{n(t+1)}{n+1}\right)^{n}+t=0$. Multiplying through by $\left(n+1\right)^{n+1}$, we obtain $0$ $=n^{n+1}\left(t+1\right)^{n+1}-\left(n+1\right)n^{n}\left(t+1\right)^{n+1}+t\left(n+1\right)^{n+1}$ $=\left(n^{n+1}-\left(n+1\right)n^{n}\right)\left(t+1\right)^{n+1}+t\left(n+1\right)^{n+1}$ $=t\left(n+1\right)^{n+1}-n^{n}\left(t+1\right)^{n+1}$. As a polynomial in $t$, the coefficient of $t$ is $\left(n+1\right)\left(\left(n+1\right)^{n}-n^{n}\right)>0$, and all other coefficients (including the constant term) are negative. By Descartes' rule of signs, either two or zero positive values of $t$ are possible, counting multiplicities. Observe that $t=\dfrac{1}{n}$ is a zero of $t\left(n+1\right)^{n+1}-n^{n}\left(t+1\right)^{n+1}$. And it is also a zero of $\dfrac{\mathrm{d}}{\mathrm{d}t}\left(t\left(n+1\right)^{n+1}-n^{n}\left(t+1\right)^{n+1}\right)=\left(n+1\right)^{n+1}-\left(n+1\right)n^{n}\left(t+1\right)^{n}$. Therefore, the possible positive values of $t$ are both $\dfrac{1}{n}$.
In the original characteristic polynomial, $t=\dfrac{n-1}{n^{2}}\ne\dfrac{1}{n}$, so it has no repeated roots. But we can continue to show the more general result. If $t=\dfrac{1}{n}$, then $r=\dfrac{n\left(t+1\right)}{n+1}=1$ and $\chi\left(x\right)=x^{n}-\dfrac{1}{n}\left(x^{n-1}+\cdots+1\right)=\dfrac{1}{n}\left(x-1\right)\left(nx^{n-1}+\left(n-1\right)x^{n-2}+\cdots+2x+1\right)$ by synthetic division, say. But $1$ cannot be a repeated root of $\chi(x)$ since $n*1^{n-1}+\left(n-1\right)*1^{n-2}+\cdots+2*1+1>0$.
Formulas for $C_i$
Formula in terms of all roots
We have $P(p)={\displaystyle \sum_{i=1}^{n}}C_{i}r_{i}^{p}$ and $P\left(1-i\right)=\dfrac{1}{2}$ for $1\le i\le n$. This yields a linear system of $n$ equations in $n$ unknowns $C_{i}$, so linear algebra can help us approach this. The equations are \begin{align*}C_{1}+\cdots+C_{n}&=\dfrac{1}{2}\\\dfrac{1}{r_{1}}C_{1}+\cdots+\dfrac{1}{r_{n}}C_{n}&=\dfrac{1}{2}\\\vdots&\\\dfrac{1}{r_{1}^{n-1}}C_{1}+\cdots+\dfrac{1}{r_{n}^{n-1}}C_{n}&=\dfrac{1}{2}\end{align*} These can be rewritten into the matrix equation $$\begin{bmatrix}r_{1}^{0} & \cdots & r_{n}^{0}\\ \vdots & \ddots & \vdots\\ r_{1}^{1-n} & \cdots & r_{n}^{1-n} \end{bmatrix}\begin{bmatrix}C_{1}\\ C_{2}\\ \vdots\\ C_{n} \end{bmatrix}=\begin{bmatrix}\dfrac{1}{2}\\ \dfrac{1}{2}\\ \vdots\\ \dfrac{1}{2} \end{bmatrix}$$
Incidentally, the square matrix is the transpose of the Vandermonde matrix with $\alpha_{i}=r_{i}^{-1}$. I had tried to form the augmented matrix and to use the beginning steps of Gaussian elimination, but it was too complicated to work through manually. Since the roots $r_{i}$ are distinct, the Vandermonde matrix is invertible and hence so is its transpose. Unfortunately, the inverse of a Vandermonde matrix is quite complicated, so I will not be deriving $C_{i}=\dfrac{1}{2}r_{i}^{n-1}{\displaystyle \prod_{j\ne i}}\dfrac{1-r_{j}}{r_{i}-r_{j}}$ in this answer. It can be confirmed for small $n$ by a CAS like Mathematica. You can see verification for $n\le7$ online.
Formula in terms of $n$ and $r_i$
Taking that formula on faith, we can exploit the symmetry to write $C_{i}$ in terms of the elementary symmetric polynomials in the $r_i$: $s_{0}=1,s_{1}={\displaystyle \sum_{j=1}^{n}}r_{j},\ldots,s_{n}={\displaystyle \prod_{j=1}^{n}}r_{j}$. We have \begin{align*}C_{i}&=\dfrac{1}{2}r_{i}^{n-1}{\displaystyle \prod_{j\ne i}}\dfrac{1-r_{j}}{r_{i}-r_{j}}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{{\displaystyle \prod_{j}\left(1-r_{j}\right)}}{{\displaystyle \prod_{j\ne i}}\left(r_{i}-r_{j}\right)}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{\sum_{j=0}^{n}\left(-1\right)^{j}s_{j}}{\left(-1\right)^{n-1}\dfrac{s_{n}}{r_{i}}+\sum_{j=1}^{n-1}\left(-1\right)^{j-1}\left(n-j\right)r_{i}^{n-j}s_{j-1}}\text{.}\end{align*} (I suspect this is very close to the form one would obtain using the inverse of the Vandermonde matrix.)
Then we can use Vieta's formulas to write the $s_{j}$s in terms of the coefficients of the characteristic polynomial. For convenience, set $t=\dfrac{n-1}{n^{2}}$, so that the characteristic polynomial is $x^{n}-t\left(x^{n-1}+\cdots+1\right)$. This means $s_{0}=1,s_{1}=-\left(-t\right),s_{2}=-t,\ldots$ with $s_{j}=\left(-1\right)^{j+1}t$ for $j>0$. As such, \begin{align*}C_{i}&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{s_{0}+\sum_{j=1}^{n}\left(-1\right)^{j}s_{j}}{\left(-1\right)^{n-1}\dfrac{s_{n}}{r_{i}}+\left(n-1\right)r_{i}^{n-1}s_{0}+\sum_{j=2}^{n-1}\left(-1\right)^{j-1}\left(n-j\right)r_{i}^{n-j}s_{j-1}}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{1+\sum_{j=1}^{n}\left(-1\right)^{j}\left(-1\right)^{j+1}t}{\left(-1\right)^{n-1}\dfrac{\left(-1\right)^{n+1}t}{r_{i}}+\left(n-1\right)r_{i}^{n-1}+\sum_{j=2}^{n-1}\left(-1\right)^{j-1}\left(n-j\right)r_{i}^{n-j}\left(-1\right)^{j}t}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{1-\sum_{j=1}^{n}t}{\dfrac{t}{r_{i}}+\left(n-1\right)r_{i}^{n-1}-\sum_{j=2}^{n-1}\left(n-j\right)r_{i}^{n-j}t}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{1-nt}{\dfrac{t}{r_{i}}+\left(n-1\right)r_{i}^{n-1}-t\sum_{k=1}^{n-2}kr_{i}^{k}}\text{.}\end{align*} The remaining summation resembles a geometric series, and indeed there is trick (often used when dealing with generating functions) involving the derivative of the geometric series formula (see this MathSE question) which allows us to write ${\displaystyle \sum_{k=1}^{n-2}}kr_{i}^{k}=\dfrac{r_{i}+\left(1-n-2r_{i}+nr_{i}\right)r_{i}^{n-1}}{\left(1-r_{i}\right)^{2}}$. Then we have: \begin{align*}C_{i}&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}}{1-r_{i}}*\dfrac{1-nt}{\dfrac{t}{r_{i}}+\left(n-1\right)r_{i}^{n-1}-t\dfrac{r_{i}+\left(1-n-2r_{i}+nr_{i}\right)r_{i}^{n-1}}{\left(1-r_{i}\right)^{2}}}\\&=\dfrac{1}{2}r_{i}^{n-1}\dfrac{1-nt}{\dfrac{t}{r_{i}}\left(1-r_{i}\right)+\left(n-1\right)r_{i}^{n-1}\left(1-r_{i}\right)-t\dfrac{r_{i}+\left(1-n-2r_{i}+nr_{i}\right)r_{i}^{n-1}}{\left(1-r_{i}\right)}}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n-1}\left(1-nt\right)\left(1-r_{i}\right)}{\dfrac{t}{r_{i}}\left(1-r_{i}\right)^{2}+\left(n-1\right)r_{i}^{n-1}\left(1-r_{i}\right)^{2}-t\left(r_{i}+\left(1-n-2r_{i}+nr_{i}\right)r_{i}^{n-1}\right)}\\&=\dfrac{1}{2}\dfrac{r_{i}^{n}\left(1-nt\right)\left(1-r_{i}\right)}{t\left(1-r_{i}\right)^{2}+\left(n-1\right)r_{i}^{n}\left(1-r_{i}\right)^{2}-t\left(r_{i}^{2}+\left(1-n-2r_{i}+nr_{i}\right)r_{i}^{n}\right)}\text{.}\end{align*}
This can be written in many different forms, especially after substituting $\dfrac{n-1}{n^{2}}$ back in for $t$. One form (for $n>1$) is $$\boxed{C_{i}=\dfrac{n}{2\left(n-1\right)}*\dfrac{1-r_{i}}{r_{i}^{-n}-2r_{i}^{1-n}+n^{2}+n-1-\left(2n^{2}+n-2\right)r_{i}+n^{2}r_{i}^{2}}}$$
Large $p$ approximation
As discussed in On finding the largest root of a polynomial by Davenport and Mignotte and New Bounds on the Real Polynomial Roots by Prodanov, Cauchy showed that the absolute values of the zeros of $a_{n}x^{n}+\cdots+a_{0}$ are bounded by the unique positive root of $\left|a_{n}\right|x^{n}-\left(\left|a_{n-1}\right|x^{n-1}+\cdots+\left|a_{0}\right|\right)$. As a consequence, the unique positive root $r\left(n\right)$ of the characteristic polynomial $x^{n}-\dfrac{n-1}{n^{2}}\left(x^{n-1}+\cdots+1\right)$ has the greatest absolute value of any root. However, we need it to be strictly greater than that of the other roots in order to justify a large-$p$ approximation of the form $C\left(n\right)*r\left(n\right)^{p}$.
We can extract more information from a proof of Cauchy's result (see Theorem 8.1.3 in Analytic Theory of Polynomials by Rahman and Schmeisser): If $r$ is a zero of $x^{n}-\dfrac{n-1}{n^{2}}\left(x^{n-1}+\cdots+1\right)$, then $r^{n}=\dfrac{n-1}{n^{2}}\left(r^{n-1}+\cdots+1\right)$, so that $\left|r\right|^{n}=\dfrac{n-1}{n^{2}}\left|r^{n-1}+\cdots+1\right|\le\dfrac{n-1}{n^{2}}\left(\left|r\right|^{n-1}+\cdots+1\right)$ by the triangle inequality, so that $x=\left|r\right|$ makes the characteristic polynomial nonpositive. If we had $\left|r\right|=\left|r(n)\right|$ then we would need equality in the triangle inequality, which would require the powers of $r$ to all point in the same direction. This is impossible (for $n>1$) if $r$ is not a positive real by De Moivre's theorem or similar.
In conclusion, for large $p$ we can indeed approximate $P(p)$ by $C(n)*r(n)^{p}$ where $r(n)$ is the unique positive root of the characteristic polynomial, and $C(n)$ is given by, say, $C(n)=\dfrac{n}{2\left(n-1\right)}*\dfrac{1-r(n)}{r(n)^{-n}-2r(n)^{1-n}+n^{2}+n-1-\left(2n^{2}+n-2\right)r(n)+n^{2}r(n)^{2}}$.
However, it would be helpful to be able to approximate $r(n)$, too. For even moderately large $n$, $r(n)$ is just barely less than $1$. We can use Newton's method starting at $1$ to get approximations. Our first approximation is $1$. Then a better approximation is $1-\dfrac{1^{n}-\dfrac{n-1}{n^{2}}\left(1^{n-1}+\cdots+1\right)}{n1^{n-1}-\dfrac{n-1}{n^{2}}\left(\left(n-1\right)1^{n-2}+\cdots1\right)}$ $=1-\dfrac{1-\dfrac{n-1}{n^{2}}*n}{n-\dfrac{n-1}{n^{2}}*\dfrac{\left(n-1\right)n}{2}}$ $=1-\dfrac{1-\dfrac{n-1}{n}}{n-\dfrac{\left(n-1\right)^{2}}{2n}}$ $=1-\dfrac{2}{n^{2}+2n-1}$.
This approximation $r_{\approx}(n)$ is already quite good (relative error less than $0.5\%$) for $n\ge10$. The corresponding approximation to $C(n)$ simplifies to $C_{\approx}(n)=\dfrac{\left(n^{2}+2n-1\right)^{n+1}\left(n^{2}+2n-3\right)}{2\left(n-1\right)\left(n^{2}+4n-1\right)r_{\approx}(n)^{n}}$.
Now, $C_{\approx}(n)r_{\approx}(n)^{p}$ is not a good approximation to $P(p)$ when $p$ is much greater than $n$. For example, when $n=10$ and $p=100$, we have about $35\%$ relative error. But $n=70$ and $p=100$ brings the relative error below $2\%$, so $C_{\approx}(n)r_{\approx}(n)^{p}$ is a straightforward way to approximate $P(p)$ for large $n$ and $p$ of comparable size.