# Solve[polynomial, x, Reals] doesn't get all real roots or correct ones?

My main confusion is about the difference between the two code blocks at the end of this long spiel, but the spiel contains the code to create the polynomial if its coefficients are helpful for explaining the difference. Thanks in advance!

I have this polynomial, eu1-eu2, of degree n (where n is a multiple of 3 between 3 and 93, inclusive).

solutions = Simplify[Solve[eu1 == eu2]];
solutions = p /. solutions;
Select[solutions, Element[#, Reals] &]


and

p /. Simplify[Solve[eu1 == eu2, p, Reals]]


should return the same roots, from what I can tell. And they do, for lower values of n. But as n increases, at some point they start returning slightly different roots, and later on, the first code block returns a bunch of different roots while the second one usually returns one root repeated many times. As n increases, the second one's repeated root seems to be farther away from any of the roots gotten from the first one.

The first one probably gets the true roots, since the size of "solutions" (which includes complex roots) equals n. But then what number is the second one getting?

Edit: 30 is the smallest n for which the number of roots is the same but the roots themselves are slightly different. Here the first method produces {0.333794} whereas the second produces {0.333788}.

33 is the smallest n for which the number of roots is different for each method. Here the first method produces {0.251609, 0.251683, 0.273157, 0.274114, 0.3336} whereas the second produces {0.33363}.

Edit: The code that sets up the polynomial is as follows:

n = 93;
q1 = (n - Mod[n, 3])/3;
q2 = q1;
upperbound = 0.5;

(mat1 = Array[0, {n, 3}];
prob1 = Array[0, {n, 1}];
prob2 = Array[0, {n, 1}];
prob3 = Array[0, {n, 1}];
For[m = 0, m <= n - 1, m++,
If[m + 1 - q1 != 0,
prob1[[m + 1, 1]] = q1/(m + 1);
prob2[[m + 1,
1]] = ((q2 - (n - m - 1))/(m + 1 - q1))*(1 - q1/(m + 1));
prob3[[m + 1, 1]] = 1 - prob1[[m + 1, 1]] - prob2[[m + 1, 1]];
mat1[[m + 1]] = {prob1[[m + 1, 1]], prob2[[m + 1, 1]],
prob3[[m + 1, 1]]},
mat1[[m + 1]] = {0, 0, 0}]])

(mat = Array[0, {n, 3}];
For[m = 0, m <= n - 1, m++,
If[m + 1 <= q1,
mat[[m + 1]] = {1, 0, 0},
mat[[m + 1]] = mat1[[m + 1]]]])

(mat3 = Array[0, {n, 3}];
prob4 = Array[0, {n, 1}];
prob5 = Array[0, {n, 1}];
prob6 = Array[0, {n, 1}];
For[m = 0, m <= n - 1, m++,
If[m + 1 - q2 != 0,
prob4[[m + 1,
1]] = ((q1 - (n - m - 1))/(m + 1 - q2))*(1 - q2/(m + 1));
prob5[[m + 1, 1]] = q2/(m + 1);
prob6[[m + 1, 1]] = 1 - prob4[[m + 1, 1]] - prob5[[m + 1, 1]];
mat3[[m + 1]] = {prob4[[m + 1, 1]], prob5[[m + 1, 1]],
prob6[[m + 1, 1]]},
mat3[[m + 1]] = {0, 0, 0}]])

(mat2 = Array[0, {n, 3}];
For[m = 0, m <= n - 1, m++,
If[m + 1 <= q2,
mat2[[m + 1]] = {0, 1, 0},
mat2[[m + 1]] = mat3[[m + 1]]]])

u = {{1 - p}, {p}, {0}};

(eu1 = 0;
For[m = 0, m <= n - 1, m++,
eu1 = eu1 + (mat[[m + 1]].u)*(p/upperbound)^
m*(1 - p/upperbound)^(n - m - 1)])
eu1 = Simplify[eu1];

(eu2 = 0;
For[m = 0, m <= n - 1, m++,
eu2 = eu2 + (mat2[[m + 1]].u)*(p/upperbound)^(n - m - 1)*(1 -
p/upperbound)^m])
eu2 = Simplify[eu2];


Also note that the following code results in the same solutions as the first code block:

Simplify[Solve[eu1 - eu2 == 0 && p < 15000, p, Complexes]]


Whereas this code results in the same solutions as the second code block:

Simplify[Solve[eu1 - eu2 == 0 && p < 15000, p, Reals]]


So it's something about allowing the function values to be complex vs. not allowing them to be complex. But I don't understand why my function would ever be complex when p is restricted to be real, since all of the coefficients are real.

Also, how exactly does it restrict function values to be real while it's solving for the roots? And why are the roots different from the case where complex function values are allowed?

Sorry that this is kind of a mess, but I would really appreciate any help you can provide!

-
It would help if you could post the smallest example where this discrepancy takes place (a list of coefficents would do as well). BTW, eu1==eu2 is an equation (eu1-eu2 is the polynomial - tonight I feel nitpicking :-) ) –  Peltio Oct 23 at 21:12
–  Sjoerd C. de Vries Oct 23 at 22:36
Impossible to figure out without a concrete example. –  Daniel Lichtblau Oct 24 at 17:03
@Janet Lu - you must give us eu1-eu2 for the n=30 case. Maybe working in machine precision the two ways to compute the approximate roots accumulate errors in a different way and some real roots get 'lost' because they end up complex with very small imaginary parts or they appear to be multiple roots when they are not? When I said "list of coefficients" I meant all the coefficient of the polynomial, not its degree only. –  Peltio Oct 24 at 17:34
Try Select[x /. Solve[-1 - 6 x + 8 x^3 == 0, x], Element[RootReduce[#], Reals] &] instead of Select[x /. Solve[-1 - 6 x + 8 x^3 == 0, x], Element[#, Reals] &]. Element[#, Reals]& may not return True or False if it can't decide. –  RiemannZeta Oct 24 at 19:49