The command roots was specifically designed for finding rational roots:
 

P:=`*`('(x*rand(-99..99)()+rand(-99..99)())^op(rand(1..4)(),[1$3,2])'$8)*
   `*`('(x^2*rand(-99..99)()+rand(-99..99)())^op(rand(1..4)(),[1$3,2])'$8);
expand(P);
roots(expand(P));

Thumb if You like.

Note: the use of has can be weak.

seq(i = nops(ListTools[Classify](has, Set, i)[true]), i = (`union`@op)(Set));

seq(i = nops(ListTools[Classify]((s,i)->member(i,s), Set, i)[true]), i = (`union`@op)(Set));

This one looks cleanest:

map((i->i) = ListTools[Occurrences], ((`union`@op) , [op])(Set), member);

Without ListTools:

F:=(Set::{set,list}({set,list})) -> map((i->i) = (i-> (nops@select)(has, Set, i)), ({op}@map)(op, Set));
F(Set);

F:=(Set::{set,list}({set,list})) -> map((i->i) = (i-> (nops@select)((s,i)->member(i,s), Set, i)), ({op}@map)(op, Set));
F(Set);

Thumb if You like.

This is the way to these substitutions, although I do not approve of using eval:

expression_ABC:=eval(expression,[x^(-(3*c)/2 + a/2 + 1/2)*(c + a + 1)*b^(-(3*c + a + 1)/(2*c))=A,c*x^(a/2 + 1/2 - c/2)*b^(-(c + a + 1)/(2*c))=B]);
simplify(%, [((a + 1)*(c + a + 1)*(2*c + a + 1))=C]);

Thumb if You like.

There was a typo in your code.

MapleQuestions_227394.mw

Thumb if You like.

I would expect LCState is a variable that belongs the module containing the GetState proc. That would be one of the modules:  RandomTools:-LinearCongruence and RandomTools.

MapleQuestions_227381.mw

I do not have the shoot library. Where is it stored?

Perhaps this gives the correct values:

F := proc(f) Re(evalf(convert(f, Sum), 20)); end proc;
seq(sum(((binomial(2*k, k))/(4^k))^d, k = 1 .. infinity), d = [3, 4, 5, 6]);
F~([%]);

Thumb if You like.

[seq(convert(X, set), X = n__a)];
map(convert, [n__a], set);

or

[seq(convert(X, list), X = n__a)];
map(convert, [n__a], list);

The anser appears to be:
E__fv = -48.46001884+(6.283185307*I)*n, E__fv = -46.13690406+3.141592654*I+(6.283185307*I)*n;
So the first one gives
E__fv = -48.46001884;
as the one real solution.


Thumb if you like.

 frontend(collect, [eqq, [d+H, H, A, a], factor], [({Non})(identical(d+H)), {}]);

 subs(`d+H` = d+H, collect(simplify(eqq, [d+H-`d+H`], [d, `d+H`]), [`d+H`, H, A, a], factor));

Thumb if You like.

Edit:
I also managed to obtain this form:

https://www.mapleprimes.com/DocumentFiles/260082_Answer/Annotation_2019-06-19_230004.png

#1:
L := [[1,2,3],[7,8,9],[13,12,11]] ;
map2(op, 3, L);

F := proc (R0, M, K) local X, Y, beta; unapply([(M*X + R0)*cos(beta*K), (M*X + R0)*sin(beta*K), M*Y], X, Y, beta) end proc;
f := F(1, 1, (Pi/2)/0.5);
A := U[1,6](x, theta);
cylinderplot(A - 0.1, theta = 0..2*Pi, x = 0..0.5, grid = [50,50], transparency = 0.3, scaling = constrained);
plottools[transform](f)(%);

Thumb if you like.

I believe there is a Windows setting that controls the presentation of numbers. Could it be this control affects Maple too?

The reason for this is the sqrt command. Instead use ^(1/2).

Thumb if you like.

and you can obtain the following form of your system:

{p^3*c[2]+p^2*q*c[1]+p*q^2*c[2]+q^3*c[1]-p*c[4], -p^3*c[1]+p^2*q*c[2]-p*q^2*c[1]+q^3*c[2]-q*c[4]};

Quite straightforward from here.