@ Sorry. What does "previous and more robust" mean? If I understand right, do you mean that the current version is less robust instead? 

Besides, it seems that neither the built-in Optimization package nor the additional GlobalOptimization package is well maintained; I cannot find any relevant news in most recent "Release Notes". (Actually, in some cases, MatLab's Optimization Toolbox and Global Optimization Toolbox are more powerful, but they are not easy-to-use for advanced users that are just familiar with Maple syntax.)

@Carl Love Amazing! Your code works well. Many thanks!
In addition, in that link, the (Mathematica) example uses 

but it appears that there is not a similar built-in functionality (Am I wrong?) in Maple. So, to verify the result, I have to copy them to Maple (and then execute `convert/FromMma`), which is more or less inconvenient. Since I never want to make an ugly procedure again, can this be implemented in a single line of input elegantly?
Unfortunately, the following implementation is quite inelegant, and does not work in general cases: 

(*  list all possible trees formed by certain binary combiners inefficiently, where leaves are distributed keeping the order at the deepest level (if possible, as the result below goes wrong) *)
(z -> ListTools:-FlattenOnce(((y -> ListTools:-MakeUnique(ListTools:-FlattenOnce(`~`[x -> ['[x[() .. k - 1][], F(x[k], x[k + 1]), x[k + 2 .. ()][]]' $ ('k' = 1 .. numelems(x) - 1)]](y))))@@(numelems(z[]) - 1))(z)))([[a, b, c, d, e]]);
 = 
    [F(F(F(F(a, b), c), d), e), F(F(F(a, b), c), F(d, e)), 

      F(F(F(a, b), F(c, d)), e), F(F(a, b), F(F(c, d), e)), 

      F(F(a, b), F(c, F(d, e))), F(F(F(a, F(b, c)), d), e), 

      F(F(a, F(b, c)), F(d, e)), F(F(a, F(F(b, c), d)), e), 

      F(a, F(F(F(b, c), d), e)), F(a, F(F(b, c), F(d, e))), 

      F(F(a, F(b, F(c, d))), e), F(a, F(F(b, F(c, d)), e)), 

      F(a, F(b, F(F(c, d), e))), F(a, F(b, F(c, F(d, e))))]

Any ideas?

@nm Actually, Mma's LeafCount has an alternate (open-source) official version (though this may not give the total number of indivisible subexpressions): The automatic complexity function.
Besides, Maple's dismantle should be used when the "size" is measured (since the `rtable`, especially the `Matrix`, is not equivalent to nested lists). However, it's hard to sort a list of expressions by Depth and then LeafCount likewise in Maple.

@C_R I find a way to "simplify" these, though it is not guaranteed in general. Fortunately, it works here: 

alias(seq(alpha[i] = RootOf(_Z^3 + _Z - 1, index = i), i = 1 .. 3)):
f__5 := collect(evala(collect(f__0, {y, x}, 'distributed', n -> convert(RootOf(PolynomialTools[MinimalPolynomial](evalf(n), _X, 3), _X, evalf(n)), RootOf))), x, evala);
 = 
[ 3   1 /        2                 \          
[x  + - \alpha[1]  + 2 alpha[1] + 2/ (y + 2) x
[     3                                       

     1 /          2               \ / 3                     \   3
   - - \2 alpha[1]  + alpha[1] + 4/ \y  - 2 alpha[1] + y + 1/, x 
     3                                                           

     1 /        2                 \          
   + - \alpha[3]  + 2 alpha[3] + 2/ (y + 2) x
     3                                       

     1 /          2               \ / 3                     \   3
   - - \2 alpha[3]  + alpha[3] + 4/ \y  - 2 alpha[3] + y + 1/, x 
     3                                                           

     1 /        2                 \          
   + - \alpha[2]  + 2 alpha[2] + 2/ (y + 2) x
     3                                       

     1 /          2               \ / 3                     \]
   - - \2 alpha[2]  + alpha[2] + 4/ \y  - 2 alpha[2] + y + 1/]
     3                                                       ]


verify(p, -3*`?()`(`*`, f__5), equal);
 = 
                              true

Compare:  

> length[shortname](f__5);
                                                          771

> MmaTranslator[Mma][LeafCount](f__5);
                                                          247

@acer You are right. However, if one can obtain an equivalent but still readable (exact) absolute factorization of p (I mean, in another way), this is also welcome. 
Edit. I find a Magma program to do so: Programmes Magma, Maple.

@Carl Love The reason is: 

It appears that the former one does more things (e.g., prevents premature evaluation, as described in this recent question) internally, hence in my view, `$` should be a low-level level. 

@Derek Wright Thanks for your concise comparisons. Unfortunately, this link appears dead. Are such comparisons still available? By the way, it seems that Maple's numerical solvers have certain limitations on solving differential equations: A Comparison Between Differential Equation Solver Suites In MatLab, R, Julia, Python, C, Mathematica, Maple, and Fortran - Stochastic Lifestyle. Will this degrade simulation accuracy and even cause major errors in some model?

@Harry Garst In my opinion, seq should also work here, but I found some issues before, so I have to use the lower-level command `$` at last. By the way, how about the so-called Chió's condensation (compared to Dodgson condensation)?

@ Well, I'm not sure if Maple includes (optimized) implementations of those algorithms from recent research which routinely outperform the classic and standard C/Fortran methods. Any benchmarks (liike Benchmarks for SciML and Equation Solvers)? 

I don't know what "tolerance = tolerance" means. If you use "res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars})):" instead, your code will work well.

Not an answer. But if you use "[x]" (instead of "x" or "{x}"), the output does agree with each other: 

with(RealDomain):
sol1 := solve(x - 1 > 0, [x]);
sol2 := solve(1/x^2 > 0, [x]);
sol3 := solve(1/(x^2 - 1) > 0, [x]);
sol4 := solve(x^2 + 1 > 0, [x]);
sol5 := solve(x^2 + 1 > 0, x); = 
                       sol1 := [[1 < x]]

                       sol2 := [[x <> 0]]

                  sol3 := [[x < -1], [1 < x]]

                       sol4 := [[x = x]]

                           sol5 := []

@acer Thank you. Though this might produce useful solutions (especially in physics and engineering), I think that for the original problem, this does not guarantee the completeness of results.

@rcorless Thanks. I doubt that this solver will be stuck or hanged as well, but if complete exact solutions can be obtained at long last, so much the better.

@tomleslie Thank you. But what about [a = 3, b = 9, c = 3]? It seems that the "AllSolutions" opinion creates a false impression. Besides, I think that it should be "a::positive, b::positive, c::positive".
By the way, it appears that these are not enough to show that a "sixth" solution does not exist. (Right or not?)