On the distribution of the number of internal equilibria in random evolutionary games

The analysis of equilibrium points is of great importance in evolutionary game theory with numerous practical ramifications in ecology, population genetics, social sciences, economics and computer science. In contrast to previous analytical approaches which primarily focus on computing the expected number of internal equilibria, in this paper we study the distribution of the number of internal equilibria in a multi-player two-strategy random evolutionary game. We derive for the first time a closed formula for the probability that the game has a certain number of internal equilibria, for both normal and uniform distributions of the game payoff entries. In addition, using Descartes’ rule of signs and combinatorial methods, we provide several universal upper and lower bound estimates for this probability, which are independent of the underlying payoff distribution. We also compare our analytical results with those obtained from extensive numerical simulations. Many results of this paper are applicable to a wider class of random polynomials that are not necessarily from evolutionary games.


Motivation
Evolutionary Game Theory (EGT) (Maynard Smith and Price 1973) has become one of the most diverse and far reaching theories in biology finding its applications in a plethora of disciplines such as ecology, population genetics, social sciences, economics and computer science (Maynard Smith 1982;Axelrod 1984;Hofbauer and Sigmund 1998;Nowak 2006;Broom and Rychtář 2013;Perc and Szolnoki 2010;Sandholm 2010;Han et al. 2017), see also recent reviews (Wang et al. 2016;Perc et al. 2017). For example, in economics, EGT has been employed to make predictions in situations where traditional assumptions about agents' rationality and knowledge may not be justified (Friedman 1998;Sandholm 2010). In computer science, EGT has been used extensively to model dynamics and emergent behaviour in multiagent systems (Helbing et al. 2015;Tuyls and Parsons 2007;Han 2013). Furthermore, EGT has provided explanations for the emergence and stability of cooperative behaviours which is one of the most well-studied and challenging interdisciplinary problems in science (Pennisi 2005;Hofbauer and Sigmund 1998;Nowak 2006). A particularly important subclass in EGT is random evolutionary games in which the payoff entries are random variables. They are useful to model social and biological systems in which very limited information is available, or where the environment changes so rapidly and frequently that one cannot describe the payoffs of their inhabitants' interactions (May 2001;Fudenberg and Harris 1992;Han et al. 2012;Gross et al. 2009;Galla and Farmer 2013).
Similar to the foundational concept of Nash equilibrium in classical game theory (Nash 1950), the analysis of equilibrium points is of great importance in EGT. It provides essential understanding of complexity in a dynamical system, such as its behavioural, cultural or biological diversity (Haigh 1988(Haigh , 1990Broom et al. 1997;Broom 2003;Traulsen 2010, 2014;Han et al. 2012;Han 2015, 2016;Broom and Rychtář 2016). A large body of literature has analysed the number of equilibria, their stability and attainability in concrete strategic scenarios such as the public goods game and its variants, see for example Broom et al. (1997), Broom (2000), ), Peña (2012, Peña et al. (2014) and Sasaki et al. (2015). However, despite their importance, equilibrium properties in random games are far less understood with, to the best of our knowledge, only a few recent efforts Traulsen 2010, 2014;Han et al. 2012;Galla and Farmer 2013;Han 2015, 2016;Broom and Rychtář 2016). One of the most challenging problems in the study of equilibrium properties in random games is to characterise the distribution of the number of equilibria (Gokhale and Traulsen 2010;Han et al. 2012): What is the distribution of the number of (internal) equilibria in a d-player random evolutionary game and how can we compute it?
This question has been studied in the literature to some extent. For example, in Traulsen (2010, 2014) and Han et al. (2012), the authors studied this question with a small number of players (d ≤ 4) and only focused on the probability of attaining the maximal number of equilibrium points, i.e. p d−1 , where p m (0 ≤ m ≤ d − 1) is the probability that a d-player game with two strategies has exactly m internal equilibria. These works use a direct approach by analytically solving a polynomial equation, expressing the positivity of its zeros as domains of conditions for the coefficients and then integrating over these domains to obtain the corresponding probabilities. However, it is impossible to extend this approach to games with a large or arbitrary number of players as in general, a polynomial of degree five or higher is not analytically solvable (Abel 1824). In more recent works Han 2015, 2016;Duong et al. 2017), we have established the links between random evolutionary games, random polynomial theory (Edelman and Kostlan 1995) and classical polynomial theory (particularly Legendre polynomials), employing techniques from the latter to study the expected number of internal equilibria, E. More specifically, we provided closed form formulas for E, characterised its asymptotic limits as the number of players in the game tends to infinity and investigated the effect of correlation in the case of correlated payoff entries. On the one hand, E offers useful information regarding the macroscopic, average behaviour of the number of internal equilibria a dynamical system might have. On the other hand, E cannot provide the level of complexity or the number of different states of biodiversity that will occur in the system. In these situations, details about how the number of internal equilibrium points distributed is required. Furthermore, as E can actually be derived from p m using the formula E = d−1 m=0 mp m , a closed form formula for p m would make it possible to compute E for any d, hence filling in the gap in the literature on computing E for large d (d ≥ 5). Therefore, it is necessary to estimate p m .

Summary of main results
In this paper, we address the above question by providing a closed-form formula for the probability p m (0 ≤ m ≤ d − 1). Our approach is based on the links between random polynomial theory and random evolutionary game theory established in our previous work Han 2015, 2016). That is, an internal equilibrium in a d-player game with two strategies can be found by solving the following polynomial equation (detailed derivation in Sect. 2), where β k = A k − B k , with A k and B k being random variables representing the entries of the game payoff matrix. We now summarise the main results of this paper. Detailed derivations and proofs will be given in subsequent sections. The first main result is an explicit formula for the probability distribution of the number of internal equilibria.
This theorem, which is stated in detail in Theorem 4 in Sect. 3, is derived from a more general theorem, Theorem 3, where we provide explicit formulas for the probability p m,2k,n−m−2k that a random polynomial of degree n has m (0 ≤ m ≤ n) positive, 2k (0 ≤ k ≤ n−m 2 ) complex and n − 2m − 2k negative roots. Note that results from Theorem 3 are applicable to a wider class of general random polynomials, i.e. beyond those derived from evolutionary random games considered in this work.
Theorem 1 is theoretically interesting and can be used to compute p m , 0 ≤ m ≤ d − 1 for small d. We use it to compute all the probabilities p m , 0 ≤ m ≤ d − 1, for d up to 5, and compare the results with those obtained through extensive numerical simulations (for validation). However, when d is larger it becomes computationally expensive to compute these probabilities using formula (2) because one needs to calculate all the probabilities p m,2k,d−1−2k , 0 ≤ k ≤ n−m 2 , which are complex multiple integrals. To overcome this issue, in Sect. 5, we develop our second main result, Theorem 2 below, which offers simpler explicit estimates of p m in terms of d and m. The main idea in developing this result is employing the symmetry of the coefficients β k . Specifically, we consider two cases Case 2: P(β k > 0) = α and P(β k < 0) = 1 − α, for all k = 0, . . . , d − 1 and for some 0 ≤ α ≤ 1. Note here that Case 1 is an instance of Case 2 when α = 1 2 and can be satisfied when a k and β k are exchangeable (see Lemma 1 below). Interestingly, the symmetry of β k allows us to obtain a much simpler treatment. The general case allows us to move beyond the exchangeability condition capturing the fact that different strategies might have different payoff properties.  (1) For d = 2: p 0 = α 2 + (1 − α) 2 and p 1 = 2α(1 − α); (2) For d = 3: This theorem is a summary of Theorems 6, 7 and 8 in Sect. 4 that are derived using Descartes' rule of signs and combinatorial methods. We note that results of the aforementioned theorems are applicable to a wider class of random polynomials that are not necessarily from random games.

Organisation of the paper
The rest of the paper is organised as follows. In Sect. 2, we recall and summarise the replicator dynamics for multi-player two-strategy games. The main contributions of this paper and the detailed analysis of the main results described above will be presented in subsequent sections. Section 3 is devoted to the proof of Theorem 1 on the probability distribution. The proof of Theorem 2 will be given in Sect. 4. In Sect. 5 we show some numerical simulations to demonstrate analytical results. In Sect. 6, further discussions are given. Finally, Appendix 1 contains proofs of technical results from previous sections.

Replicator dynamics
A fundamental model of evolutionary game theory is replicator dynamics (Taylor and Jonker 1978;Zeeman 1980;Hofbauer and Sigmund 1998;Schuster and Sigmund 1983;Nowak 2006), describing that whenever a strategy has a fitness larger than the average fitness of the population, it is expected to spread. For the sake of completeness, below we derive the replicator dynamics for multi-player two-strategy games.
Consider an infinitely large population with two strategies, A and B. Let x, 0 ≤ x ≤ 1, be the frequency of strategy A. The frequency of strategy B is thus (1 − x). The interaction of the individuals in the population is in randomly selected groups of d participants, that is, they play and obtain their fitness from d-player games. The game is defined through a (d − 1)-dimensional payoff matrix (Gokhale and Traulsen 2010), as follows. Let A k (respectively, B k ) be the payoff that an A-strategist (respectively, a B-strategist) obtained when playing with a group of d − 1 players that consists of k A-strategists. In this paper, we consider symmetric games where the payoffs do not depend on the ordering of the players. Asymmetric games will be studied in a forthcoming paper. In the symmetric case, the probability that an A strategist interacts with k other A strategists in a group of size d − 1 is Thus, the average payoffs of A and B are, respectively The replicator equation of a d-player two-strategy game is given by (Hofbauer and Sigmund 1998;Sigmund 2010; Gokhale and Traulsen 2010) Since x = 0 and x = 1 are two trivial equilibrium points, we focus only on internal ones, i.e. 0 < x < 1. They satisfy the condition that the fitnesses of both strategies are the same, i.e. π A = π B , which gives rise to Using the transformation y = x 1−x , with 0 < y < +∞, dividing the left hand side of the above equation by (1− x) d−1 we obtain the following polynomial equation for y Note that this equation can also be derived from the definition of an evolutionarily stable strategy (ESS), an important concept in EGT (Maynard Smith 1982), see e.g., Broom et al. (1997). Note however that, when moving to random evolutionary games with more than two strategies, the conditions for ESS are not the same as for those of stable equilibrium points of replicator dynamics. As in Gokhale and Traulsen (2010), Han (2015, 2016), we are interested in random games where A k and B k (thus β k ), for 0 ≤ k ≤ d − 1, are random variables. In Sect. 3 where we provide estimates for the number of internal equilibria in a d-player two-strategy game, we will use the information on the symmetry of β k . The following lemma gives a necessary condition to determine when the difference of two random variables is symmetrically distributed.
Lemma 1 (Duong et al. 2017, Lemma 3.5) Let X and Y be two exchangeable random variables, i.e. their joint probability distribution f X , . In addition, if X and Y are i.i.d then they are exchangeable.
For the sake of completeness, the proof of this Lemma is provided in Sect. 1.

The distribution of the number of positive zeros of random polynomials and applications to EGT
This section focuses on deriving the distribution of the number of internal equilibria of a d-player two-strategy random evolutionary game. We recall that an internal equilibria is a real and positive zero of the polynomial P(y) in (5). We denote by κ the number of positive zeros of this polynomial. For a given m, 0 ≤ m ≤ d − 1, we need to compute the probability p m that κ = m. To this end, we first adapt a method introduced in Zaporozhets (2006) (see also Butez and Zeitouni 2017;Götze et al. 2017 for its applications to other problems) to establish a formula to compute the probability that a general random polynomial has a given number of real and positive zeros. Then we apply the general theory to the polynomial P.

The distribution of the number of positive zeros of a random polynomial
Consider a general random polynomial P(t) = ξ 0 t n + ξ 1 t n−1 + · · · + ξ n−1 t + ξ n .
As consequences, (1) The probability that P has m positive zeros is (2) In particular, the probability that P has the maximal number of positive zeros is Proof The reference (Zaporozhets 2006, Theorem 1) provides a formula to compute the probability that the polynomial P has n − 2k real and 2k complex roots. In the present paper, we need to distinguish between positive and negative real zeros. We now sketch and adapt the proof of Theorem 1 of Zaporozhets (2006) to obtain the formula (9) for the probability that the polynomial P has m positive, 2k complex and n −m −2k negative roots. Consider a (n + 1)-dimensional vector space V of polynomials of the form Q(t) = a 0 t n + a 1 t n−1 + · · · + a n−1 t + a n , and a measure μ on this space defined as the integral of the differential form d Q = p(a 0 , . . . , a n ) da 0 ∧ · · · ∧ da n .
Our goal is to find μ(V m,2k ) where V m,2k is the set of polynomials having m positive, 2k complex and n − m − 2k negative roots. Let Q ∈ V m,2k . Denote all zeros of Q as To find μ(V m,2k ) we need to integrate the differential form (12) over the set V m,2k .
The key idea in the proof of Theorem 1 of Zaporozhets (2006) is to make a change of coordinates (a 0 , . . . , a n ) → (a, x 1 , . . . , x n−2k , r 1 , . . . , r k , α 1 , . . . , α k ), with a = a 0 , and find d Q in the new coordinates. The derivation of the following formula is carried out in detail in Zaporozhets (2006): hence the assertion (9) follows.

The distribution of the number of internal equilibria
Next we apply Theorem 3 to compute the probability that a random evolutionary game has m, 0 ≤ m ≤ d − 1, internal equilibria. We derive formulas for the three most common cases (Han et al. 2012): . The main result of this section is the following theorem (cf. Theorem 2).

Theorem 4 The probability that a d-player two-strategy random evolutionary game
where p m,2k,d−1−m−2k is given below for each of the cases above: where σ i , for i = 0, . . . , d − 1, and Δ are given in (10)-(11) and -For the case (C3) In particular, the probability that a d-player two-strategy random evolutionary game has the maximal number of internal equilibria is: (1) for the case (C1) (2) for the case (C2) (3) for the case (C3) Note that in formulas (16)-(18) above Therefore, where Using the following formula for moments of a normal distribution, Applying Theorem 3 to the polynomial P given in (5) and using the above identity we obtain which is the desired equality (13) by definition of C and σ . ( we have (for simplicity of notation, in the subsequent computations we shorten min i∈{0,...,d−1} by min) Similarly as in the normal case, using this identity and applying Theorem 3 we obtain (3) Now we assume that A j and B j are i.i.d. uniformly distributed with the common distribution γ (x) = 1 2 1 [−1,1] (x). Since β j = A j − B j , its probability density is given by The probability density of δ j β j is

Corollary 1 The expected numbers of internal equilibria and stable internal equilibria, E(d) and S E(d), respectively, of a d-player two-strategy game, are given by
Note that this formula for E(d) is applicable for non-normal distributions, which is in contrast to the method used in previous works Han 2015, 2016) that can only be used for normal distributions. The second part, i.e. the formula for the expected number of stable equilibrium points, was obtained based on the following property of stable equilibria in multi-player two-strategy evolutionary games, as shown in Han et al. (2012, Theorem 3): Remark 1 In Theorem 4 for the case (C1), the assumption that β k 's are standard normal distributions, i.e. having variance 1, is just for simplicity. Suppose that β k 's are normal distributions with mean 0 and variance η 2 . We show that the probability p m , for 0 ≤ m ≤ d − 1, does not depend on η. In this case, the formula for p is given by (19) but with C being replaced by η 2 C. To indicate its dependence on η, we write p η . We use a change of variable a = ηã. Then
(2) Two internal equilibria: The probability that a four-player two-strategy evolutionary game has 2 internal equilibria is (3) Three internal equilibria: p 3 = p 3,0,0 The probability that a four-player two-strategy evolutionary game has 3 internal equilibria is (4) No internal equilibria: the probability that a four-player two-strategy evolutionary game has no internal equilibria is:

Universal estimates for p m
In Sect. 3, we have derived closed-form formulas for the probability distributions p m (0 ≤ m ≤ d − 1) of the number of internal equilibria. However, it is computationally expensive to compute these probabilities since it involves complex multiple-dimensional integrals. In this section, using Descartes' rule of signs and combinatorial techniques, we provide universal estimates for p m . Descartes' rule of signs is a technique for determining an upper bound on the number of positive real roots of a polynomial in terms of the number of sign changes in the sequence formed by its coefficients. This rule has been applied to random polynomials before in the literature (Bloch and Pólya 1932); however this paper only obtained estimates for the expected number of zeros of a random polynomial.
Theorem 5 (Descartes' rule of signs, see e.g., Curtiss 1918) Consider a polynomial of degree n, p(x) = a n x n + · · · + a 0 with a n = 0. Let v be the number of variations in the sign of the coefficients a n , a n−1 , . . . , a 0 and n p be the number of real positive zeros. Then (v − n p ) is an even non-negative integer.
We recall that an internal equilibrium of a d-player two-strategy game is a positive root of the polynomial P given in (5). We will apply Descartes' rule of signs to find an upper bound for the probability that a random polynomial has a certain number of positive roots. This is a problem that is of interest in its own right and may have applications elsewhere; therefore we will first study this problem for a general random polynomial of the form p(y) := n k=0 a k y k , and then apply it to the polynomial P. It turns out that the symmetry of {a k } will be the key: the asymmetric case requires completely different treatment from the symmetric one.

Estimates of p m : symmetric case
We first consider the case where the coefficients {a k } in (22) are symmetrically distributed. The main result of this section will be Theorem 6 that provides several upper and lower bounds for the probability that a d-player two strategy game has m internal equilibria. Before stating Theorem 6, we need the following auxiliary lemmas.
Proposition 1 Suppose that the coefficients a k , 0 ≤ k ≤ n in the polynomial (22) are i.i.d. and symmetrically distributed. Let p k,n , 0 ≤ k ≤ n, be the probability that the sequence of coefficients (a 0 , . . . , a n ) has k changes of signs. Then p k,n = 1 2 n n k .
Proof See Appendix 2.
The next two lemmas on the sum of binomial coefficients will be used later on.
Lemma 2 Let 0 ≤ k ≤ n be positive integers. Then it holds that n j=k j:even where it is understood that n j = 0 if j < 0. In particular, for k = 0, we get n j=0 j:even n j = n j=0 j:odd n j = 2 n−1 .
Proof See Appendix 3.
The following lemma provides estimates on the sum of the first k binomial coefficients.

Lemma 3 Let n and
where δ = 0.98 and H is the binary entropy function where 0 log 2 0 is taken to be 0. In addition, if n = 2n is even and 0 ≤ k ≤ n , we also have the following estimate (Lovász et al. 2003, Lemma 3.8.2) We now apply Proposition 1 and Lemmas 2 and 3 to derive estimates for the probability that a d-player two-strategy evolutionary game has a certain number of internal equilibria. The main theorem of this section is the following. As In addition, if d − 1 = 2d is even and 0 ≤ m ≤ d then (b) Lower-bound for p 0 and p 1 : Proof (a) This part is a combination of Decartes' rule of signs, Proposition 1 and Lemmas 2 and 3. In fact, as a consequence of this rule and by Proposition 1, we have p m ≤ j: j≥m j−m: even which is the inequality part in (29). Next, applying Lemma 2 for k = m and n = d − 1 and then Lemma 3, we obtain This proves the equality part in (29) and (30). As a result, the estimate p m ≤ 1 2 for all 0 ≤ m ≤ d − 1 is followed from (29) and (24); the estimates p d−1 ≤ 1 2 d−1 and p d−2 ≤ d−1 2 d−1 are special cases of (29) for m = d − 1 and m = d − 2, respectively.
(b) It follows from Decartes' rule of signs and Proposition 1 that (c) For d = 2: from parts (a) and (b) we have which implies that p 0 = p 1 = 1 2 as claimed.
(d) Finally, for d = 3: also from parts (a) and (b) we get so p 1 = 1 2 . This finishes the proof of Theorem 6. Remark 2 Note that in Theorem 6 we only assume that β k are symmetrically distributed but do not require that they are normal distributions. When {β k } are normal distributions, we have derived Han 2015, 2016) Fig. 2 shows that the new ones do better for m closer to 0 or d − 1 but worse for intermediate m (i.e. closer to (d − 1)/2).

Estimates of p m : general case
In the proof of Proposition 1 the assumption that {a k } are symmetrically distributed is crucial. In that case, all the 2 n binary sequences constructed are equally distributed, resulting in a compact formula for p k,n . However, when {a k } are not symmetrically distributed, those binary sequences are no longer equally distributed. Thus computing p k,n becomes much more intricate. We now consider the general case where Note that the general case allows us to move beyond the usual assumption in the analysis of random evolutionary games that all payoff entries a k 's and b k 's have the same probability distribution resulting in α = 1/2 (see Lemma 1). In the general case it only requires that all a k 's have the same distribution and all b k 's have the same distribution, capturing the fact that different strategies, i.e. A and B in Sect. 2, might have different payoff properties (e.g., defectors always have a larger payoff than cooperators in a public goods game).
The main results of this section will be Theorem 7 and Theorem 8. The former provides explicit formulas for p k,n while the latter consists of several upper and lower bounds for p m . We will need several technically auxiliary lemmas whose proofs will be given in Appendix 1. We start with the following proposition that provides explicit formulas for p k,n for k ∈ {0, 1, n − 1, n}.

Proposition 2
The following formulas hold: if n is even, if n is odd.
The computations of p k,n for other k are more involved. We will employ combinatorial techniques and derive recursive formulas for p k,n . We define u k,n = P(there are k variations of signs in {a 0 , . . . , a n } a n > 0), v k,n = P(there are k variations of signs in {a 0 , . . . , a n } a n < 0).
We have the following lemma.

Lemma 4
The following recursive relations hold: We can decouple the recursive relations in Lemma 4 to obtain recursive relations for {u k,n } and v k,n separately as follows: Lemma 5 The following recursive relations hold Proof See Appendix 6.
Using the recursive equations for u k,n and v k,n we can also derive a recursive relation for p k,n .
Proposition 3 { p k,n } satisfies the following recursive relation.
Proof See Appendix 7.
Using the transformation a k,n := 1 2 n p k,n as in the proof of Theorem 7, then a k,n = a k−2,n−2 − a k,n−2 + 2a k,n−1 , which is exactly the Chu-Vandermonde identity for m = 2 above. Then it is no surprise that in Theorem 7 we obtain that a k,n is exactly the same as the binomial coefficient In the next main theorem we will find explicit formulas for { p k,n } from the recursive formula in the previous lemma using the method of generating functions. The case α = 1 2 will be a special one.
Theorem 7 p k,n is given explicitly by: for α = 1 2 , p k,n = 1 2 n n k .

Numerical simulations
In this section, we perform several numerical (sampling) simulations and calculations to illustrate the analytical results obtained in previous sections. Figure 1 shows the values of { p m } for d ∈ {3, 4, 5}, for the three cases studied in Theorem 4, i.e., when β k are i.i.d. standard normally distributed (GD), uniformly distributed (UD1) and when β k = a k − b k with a k and β k being uniformly distributed (UD2). We compare results obtained from analytical formulas in Theorem 4 and from samplings. The figure shows that they are in accordance with each other agreeing to at least 2 digits after the decimal points. Figure

Further discussions and future research
In this paper, we have provided closed-form formulas and universal estimates for the probability distribution of the number of internal equilibria in a d-player two-strategy random evolutionary game. We have explored further connections between evolutionary game theory and random polynomial theory as discovered in our previous works Han 2015, 2016;Duong et al. 2017). We believe that the results reported in the present work open up a new exciting avenue of research in the study of equilibrium properties of random evolutionary games. We now provide further discussions on these issues and possible directions for future research.
Computations of probabilities { p m }. Although we have found analytical formulas for p m it is computationally challenging to deal with them because of their complexity. Obtaining an effective computational method for { p m } would be an interesting problem for future investigation. The payoff entries a k and b k were drawn from a normal distribution with variance 1 and mean 0 (GD) and from a standard uniform distribution (UD2). We also study the case where β k = a k − b k itself is drawn from a standard uniform distribution (UD1). Results are obtained from analytical formulas (Theorem 2) (a) and are based on sampling 10 6 payoff matrices (b) where payoff entries are drawn from the corresponding distributions. Analytical and simulations results are in accordance with each other. All results are obtained using Mathematica Quantification of errors in the mean-field approximation theory (Schehr and Majumdar 2008). Consider a general polynomial P as given in (6) with dependent coefficients, and let P m ([a, b], n) be the probability that P has m real roots in the interval [a, b] (recall that n is the degree of the polynomial, which is equal to d − 1 in Equation (1)). The mean-field theory (Schehr and Majumdar 2008) neglects the correlations between the real roots and simply considers that these roots are randomly and independently distributed on the real axis with some local density f (t) at point t, with f (t) being the density that can be computed from the Edelman-Kostlan theorem (Edelman and Kostlan 1995). Within this approximation in the large n limit, the probability P m ([a, b], n) is given by a non-homogeneous Poisson distribution, see Schehr and Majumdar (2008, Section 3.2.2 and Equation (70)). By applying the mean-field theory one can approx-  imate the probability p m that a random d-player two-strategy evolutionary game has m internal equilibria by a simpler and computationally feasible formula. However, it is unclear to us how to quantify the errors of approximation. We leave this topic for future research.
Extensions to multi-strategy games. We have focused in this paper on random games with two strategies (with an arbitrary number of players). The analysis of games with more than two strategies is much more intricate since in this case one needs to deal with systems of multi-variate random polynomials. We have provided Han 2015, 2016) a closed formula for the expected number of internal equilibria for a multi-player multi-strategy games for the case of normal payoff entries. We aim to extend the present work to the general case in future publications. In particular, Decartes' rule of signs for multi-variate polynomials (Itenberg and Roy 1996) might be used to obtain universal estimates, regardless of the underlying payoff distribution.
According to Duong and Tran (2018, Lemma 5.4) or equivalently: This finishes the proof of this lemma.
Therefore we obtain the first relationship in (33). The second one is proved similarly.

Proof of Proposition 3
From Lemmas 4 and 5 we have This finishes the proof.
According to Proposition 2 a 0,n = A n p 0,n = A n α n+1 + (1 − α) n+1 = α α 1 − α n 2 a 1,n = A n p 1,n = n if α = 1 2 , Also a k,n = 0 for k > n. Let F(x, y) be the generating function of a k,n , that is Note that in the above computations we have the following identities And finally for the last sum ∞ n=2 a k,n−1 x k y n = y(F(x, y) − g(x, y)).
For α = 1 2 , we get (1 + x) n y n = ∞ n=0 n k=0 n k x k y n , which implies that α k,n = n k . Hence for the case α = 1 2 , we obtain p k,n = 1 2 n n k .
Finding the series expansion for this case is much more involved than the previous one. Using the multinomial theorem we have Similarly, if k is odd, k = 2k + 1, then to obtain the coefficient of x k y n on the right-hand side of (40), we select (i, m, l) such that and obtain a k,n = 2 n m= n−1 2 m k , n − k − m − 1, 2m − n + 1 (−1) n−k −m−1 A 2m−n .
From a k,n we compute p k,n using the relations p k,n = a k,n A n and A 2 = 1 α(1−α) and obtain the claimed formulas. This finishes the proof of this theorem. Remark 4 We can find a k,n by establishing a recursive relation. We have It is not trivial to obtain an explicit formula from this recursive formula. However, it is easily implemented using a computational software such as Mathematica or Mathlab.