@Carl Love 

Nice and professional (as always :-)).

@digerdiga 

asympt cannot switch branches. You can obtain the same result without puiseux.

restart;
sol:=RootOf((8*n-8)*_Z^6+(n^4+36*n^2-68*n+56)*_Z^5+(n^5+10*n^4+80*n^3-200*n^2+224*n-152)*_Z^4+(n^6+28*n^5+69*n^4-268*n^3+468*n^2-356*n+200)*_Z^3+(3*n^7+32*n^6+7*n^5-204*n^4+380*n^3-544*n^2+272*n-128)*_Z^2+(3*n^8+14*n^7-20*n^6-32*n^5+252*n^4-240*n^3+304*n^2-80*n+32)*_Z-n^9-12*n^8-44*n^7-40*n^6-4*n^5-128*n^4+48*n^3-64*n^2):
f:=eval(op(sol), [_Z=z,n=1/x]):
g:=numer(eval(f, z=a/x+b+c*x)):
S:=solve({seq(coeff(g,x,k),k=0..2)}):
subs(S[],x=1/n, a/x+b+c*x);

@Kitonum 

I know. SymFun just calls solve/identity which could (and should) be improved a lot (you are free to send a SCR for it, and if successful, SymFun will work for your examples). Actually, if we are interested only in an odd/even result, this is easy to obtain.
BTW, the "hack" part was included just to answer (c). SymFun can be of course improved, but I prefer a capable solve/identity.

@Carl Love 

To find the initial point of a p-cycle simply do:

f:=y-> 4*y*(1-y):
p:=2;
solve((f@@p)(x)=x); # or fsolve

@tomleslie 

No need for DirectSearch. The min of expr^2 is of course 0 (attained for all variables 0).

@mamehrashi 

Write your problem in the text area, not in the title.
Maybe a smaller example will help.
If   P = 2*x*y + 3*x*z + 4*y*t + 5*t^2   is the polynomial in the variables x,y,z,t,
what result do you expect?

You just posted a polynomial in 16 variables, homogeneous, total degree 4,  with float coefficients.
What do you want to do with it?

@acer 

You have used the parametric version

plot3d([1,u,v],u=0..2*Pi-2*v,v=0..Pi, coords=spherical)

instead of the explicite

plot3d(1, u=0..2*Pi-2*v,v=0..Pi, coords=spherical)

It seems that you knew that the latter does not work correctly. Is it documented, or just a bug?

@Carl Love 

If we want separated colors, without interpolation, I don't think that an alternative to grid exists.

c:=piecewise(u+v<Pi, 1/6, u+v<2*Pi, 2/3, 1):
plot3d(1,u=0..2*Pi,v=0..Pi, coords=spherical, color=c, grid=[300,300], style=surface);

Compare with the default grid:

@acer Consider for simplicity  plot3d(f(x,y), x=0..1, y=0..1,  color = [c, colorscheme = [c1,c2,c3,c4]])

c = c(x,y) being the color and c1,c2,... fixed colors.

I can imagine a (useful) meaning of this for c = x (i.e. a kind of xgradient), even for c = c(x), but not for c = c(x,y).
What would be your definition? To me, it's nonsense for the moment.

plot3d( 1,u=0..2*Pi,v=0..Pi, coords=spherical, 
colorscheme= ["zgradient",[black, red, yellow, white],  colorspace="HSV"] , 
grid=[300,300]);

rr:=(x,y,z) -> piecewise(x<-1/3,0, x<1/3,1, x<=1,1):
gg:=(x,y,z) -> piecewise(x<-1/3,0, x<1/3,1, x<=1,0):
bb:=(x,y,z) -> piecewise(x<-1/3,1, x<1/3,0, x<=1,0):
plot3d(1,u=0..2*Pi,v=0..Pi, coords=spherical,
colorscheme=["xyzcoloring",[rr,gg,bb]], grid=[300,300]);

@dharr 

Yes, noticing that a "1" must exist in each row and column allows to reduce the number of the candidates.
This number can be halved (so, also the timing). And it is actually necessary to do so because otherwise each A appears twice (there are only 40320).
It is enough to change
for m to n^2 do
to

for m from k+1 to n^2 do

And preferably also k to n^2 - 1.

Edit. It's possible to reduce further the candidates by identifying two matrices A,B with A=B^+. (Maybe other symmetries too.)
BTW, for the second polynomial f:=x^8+x^4+x^3+x^2+1, A does not exist.

@Magma 

There are 2 possibilities.
1. A solution does not exist.
2. There are very few solutions and more time is needed. Actually in this case a random approach may fail, but unfortunatelly the brute force in not an option because there are > 10^11 candidates for A. [A rough estimation shows that on a regular computer almost a month could be needed].

The only solution I see is to try to reduce the number of candidates using a theoretical examination.
Or, try the "reversed" approach: take a matrix having f as characteristic polynomial (e.g. the matrix companion) and try to describe  somehow the form of its conjugates in order to find a conjugate with 10 ones.

1. Please post your code (in order to correct it). An example and the expected result would be useful too.

2. What kind of decomposition? How many terms?
What about using the elementary matrices obtained during a Gauss-Jordan elimination?

@mmcdara 

simp := eq -> `if`( eq=(0=0), undefined, eq/gcd(op(eq)));
or maybe
simp:=eq -> `if`( eq=(0=0), undefined, `if`(`*`(op(eq))=0,eq,eq/gcd(op(eq))) );