In the previous embodiment, the diagonal entries are always  0 , and should be random integers from  -9  to  9 . New version is free from this drawback:

roll:=rand(-9 .. 9):
Matrix(4, {seq((i, i)=roll(), i=1..4), seq(seq((i, j)=roll()+I*roll(), j=i+1..4), i=1..3)}, shape=hermitian);

Factoring  procedure takes 1/2 sec for this calculation, and the improvements made by Joe Riel and Carl Love even less time. See  http://www.mapleprimes.com/posts/141668-Partitions-Of-A-Natural-Number-Into-Factors

ts:=time():

Factoring(9!):

nops(%);

%%[1..20];

time()-ts;

I have already answered a similar question about the matrix  A . About the matrix  B  all is still easier. There is an inaccuracy in determining  B. Apparently instead of  -1/(2*M*(M-2))   should be -1/(M*(M-2)) , otherwise a contradiction with special cases.

@Bendesarts   Do all things similarly, only replace  eq  by  lhs(eq_liaison[1, 1])  and etc.

@Carl Love  Visually  f[`3x`,`4y`]  is identical to  f[`3x, 4y`]

@yihezhi 

f[`{6x}`]; f[`{3x,4y}`];

@Bendesarts  For some reason your last file is not downloaded. In my M 16 all works. See the attached file.

Solving_equation.mw

Here is the copy of the code:

eq := 2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*cos(gamma[1](t))*e[1]*rBTP[1]- 2*cos(alpha(t))*cos(beta(t))*h*h[1]+2*sin(beta(t))*sin(a[1])*h*rF[1]- 2*cos(alpha(t))*cos(a[1])^2*rBTP[1]*rF[1]-2*z(t)*sin(alpha(t))*cos(a[1])*rBTP[1]+ 2*z(t)*cos(alpha(t))*cos(beta(t))*h+2*sin(alpha(t))*cos(a[1])*h[1]*rBTP[1]+ 2*cos(beta(t))*cos(gamma[1](t))*e[1]*rBTP[1]+2*cos(beta(t))*cos(a[1])^2*rBTP[1]*rF[1]+ 2*sin(gamma[1](t))*e[1]*h[1]-2*z(t)*sin(gamma[1](t))*e[1]-2*cos(gamma[1](t))*e[1]*rF[1]- 2*cos(beta(t))*rBTP[1]*rF[1]-2*z(t)*h[1]-l[1]^2+ 2*sin(alpha(t))*cos(beta(t))*cos(a[1])*cos(gamma[1](t))*h*e[1]- 2*cos(alpha(t))*sin(beta(t))*sin(a[1])*sin(gamma[1](t))*e[1]*rBTP[1]- 2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*rBTP[1]*rF[1]+z(t)^2+ h[1]^2+rBTP[1]^2+e[1]^2+rF[1]^2+h^2+2*cos(alpha(t))*cos(a[1])^2* cos(gamma[1](t))*e[1]*rBTP[1]-2*cos(alpha(t))*sin(beta(t))*sin(a[1])*h[1]*rBTP[1]- 2*cos(alpha(t))*cos(beta(t))*sin(gamma[1](t))*h*e[1]-2*sin(beta(t))*sin(a[1])* cos(gamma[1](t))*h*e[1]-2*cos(beta(t))*cos(a[1])^2*cos(gamma[1](t))*e[1]*rBTP[1]+ 2*z(t)*cos(alpha(t))*sin(beta(t))*sin(a[1])*rBTP[1]-2*sin(alpha(t))*cos(beta(t))* cos(a[1])*h*rF[1]+2*sin(alpha(t))*cos(a[1])*sin(gamma[1](t))*e[1]*rBTP[1];

eq1 := collect(subs(sin(gamma[1](t)) = x, cos(gamma[1](t)) = y, eq), [x, y]);

A := coeff(eq1, x);

B := coeff(eq1, y);

C := tcoeff(eq1, [x, y]);

sol := solve(a*sin(t)+b*cos(t)+c, t);

gamma1 := unapply(eval(sol[1], {a = A, b = B, c = C}), t);

gamma2 := unapply(eval(sol[2], {a = A, b = B, c = C}), t);

diff(f(x,y),x$6);

diff(f(x,y),[x$3],[y$4]);

@ctnaz   No, the numbers  2  and  0.94  I got just by selection, banishing my program several times for different values​​.

@Preben Alsholm  Thanks for the solution! I've tried similar variants, but my mistake was that I took too little eps=10^(-6)

@Kitonum  

My first program was unsuccessful. I forgot that break keyword only interrupts the inner loop, but the program continues to run. Therefore returned not the first but the last solution.
One solution (the first one) can be found immediately when the program as a procedure to issue and uses  return keyword:

restart;
P := proc ()
local a, b, c, d, e, f, g, h;
for c from -10 to 10 do
for d from -10 to 10 do
for e from -10 to 10 do
for f from -10 to 10 do
for g from -10 to 10 do
for h from -10 to 10 do
a := c+2*d+4*e+6*f+10*g+22*h-16; b := -2*c-3*d-5*e-7*f-11*g-23*h+20;
if -10 <= a and a <= 10 and -10 <= b and b <= 10 then return [a, b, c, d, e, f, g, h]
end if end do end do end do end do end do end do
end proc:

P();

                                               [10, 10, -10, -10, -10, 8, 7, -1]

@Mikhail Drugov  Yes,  I think there is no automatic way.

@brian bovril  Of course I can. This is really easy.

Let the member of the sequence, standing on the  n-th position is  m. Then it is obvious that

1+2+...+(m-1)< n <=1+2+...+m   or   (m-1)*m/2< n <=m*(m+1)/2

Then, solving this double inequality for m> = 1, we obtain

-(1/2) + 1/2*sqrt(1 + 8*n) <= m < 1/2 + 1/2*sqrt(1 + 8*n)

Because left and right parts of the solution differ by 1, we immediately obtain the formula.

The question arises - what of it relates to Maple?

your question more specifically?