The easiest way - translation of all constants to symbols and use the concatenation operator.

Examples:

`3`+`7`+`1`=`11`;

`3+7+1`=`11`;

`(-2)`||` (-1)`/`(-5)`||` (-1)`;

add(a[i]*mul(s-b[j], j={$1..n} minus {i}), i= 1..n);

Example:

n:=5:
add(a[i]*mul(s-b[j], j={$1..n} minus {i}), i= 1..n);

PS:  In Standard instead of  minus  you can use  \  :

n := 5:  add(a[i]*mul(s-b[j], j={$1..n}\{i}), i=1..n);

There are the functions  max  and  min  but not  Max  and  Min.  

The equation  max(min(x, 1), min(x, 2))=a  is equivalent to  

min(x,1)<=a  and  min(x,2)<=a  and  (min(x,1)=a  or  min(x,2)=a)

so

solve(min(x,1)<=a  and  min(x,2)<=a  and  (min(x,1)=a  or  min(x,2)=a), x);

eq is your last equation in which all terms are moved to one side.

eq := 2*sin(alpha(t))*cos(a[1])*h[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(beta(t))*cos(gamma[1](t))*e[1]*rBTP[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*cos(beta(t))*cos(a[1])^2*rBTP[1]*rF[1]+2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*cos(gamma[1](t))*e[1]*rBTP[1]-2*z(t)*h[1]+2*sin(gamma[1](t))*e[1]*h[1]-2*z(t)*sin(gamma[1](t))*e[1]-2*cos(beta(t))*rBTP[1]*rF[1]-2*cos(gamma[1](t))*e[1]*rF[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]-2*sin(beta(t))*sin(a[1])*cos(gamma[1](t))*h*e[1]+2*cos(alpha(t))*cos(a[1])^2*cos(gamma[1](t))*e[1]*rBTP[1]+z(t)^2+h[1]^2+rBTP[1]^2+e[1]^2+rF[1]^2+h^2-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*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);  # 2 series of solutions
gamma2 := unapply(eval(sol[2], {a = A, b = B, c = C}), t);

It would be better to use  solve  instead of  isolate, because  isolate does not provide all solutions . See the example:

eq:=a*sin(t)+b*cos(t)+c=0:

isolate(eq, t);

solve(eq, t, AllSolutions);

This is a cylindrical body, in the base of which the circle  x^2+y^2=2*x  or  (x-1)^2+y^2=1, and on the top  it is limited by a paraboloid of revolution x^2+y^2=z .

plot3d(x^2+y^2, y = -sqrt(2*x-x^2) .. sqrt(2*x-x^2), x = 0 .. 2, scaling = constrained, filled = true, axes = normal, numpoints = 10000);
Volume=int(x^2+y^2, [y = -sqrt(2*x-x^2) .. sqrt(2*x-x^2), x = 0 .. 2]);

restart;

X:=[3,4,5,6,7]:

Y:=[7.42494922444550, 3.67768248674133, 2.52235142453921, 1.95610223891559, 1.61770309810016]:

A:=plot(X,Y, style=point, symbolsize=20):

CurveFitting[LeastSquares](X,Y, x, curve=a/(x-2)^0.94);

B:=plot(%, x=0..8, -1..9, discont=true,color=blue):

plots[display](A,B, scaling=constrained);  

PS.  Your attached file is blank. 

I was wrong about "a lot of time." Still, the capabilities of modern computers boggle my mind ! Less than 7 minutes the program went through  21^6 = 85766121 variants and  found all solutions. The first 10 and last 10 solutions displayed.

restart;

ts := time(): N := 0:

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 a >= -10 and a <= 10 and b >= -10 and b <= 10 then N := N+1:
Sol[N] := [a, b, c, d, e, f, g, h] end if:
end do: end do: end do: end do: end do: end do:

N; seq(Sol[i], i = 1 .. 10); seq(Sol[i], i = N-9 .. N);
time()-ts;

Of a given system  a  and  b  can be easily expressed  through the other variables, and then in the loop we find one solution:

restart;

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 a>=-10 and a<=10 and b>=-10 and b<=10 then Sol:=[a,b,c,d,e,f,g,h]: break: fi:

od: od: od: od: od: od:

Sol;

                                        [-10, -10, 10, 10, 10, 1, -7, 0]

Finding all solutions will require a lot of time.

Use  mtaylor  command. In your example,  a  can not be 0, because it is in the denominator:

mtaylor((1+1/(2*a^2)+(1/2)*b^2-1/(2*a))^a, [a=1,b=0], 5);

From the plot  

plot([ln(x+1),x^2-1], x=-2..3, -3..3);  

we can see that there are exactly two points of intersection of the graphs of functions. Solving numerically the corresponding equations, we finally obtain

RealRange(Open(-1),fsolve(ln(x+1)=x^2-1, x=-infinity..0)),RealRange(fsolve(ln(x+1)=x^2-1, x=0..infinity), infinity) ;

This can be done as follows (example):

A := plot(sqrt(x), x = 0 .. 4, tickmarks = [0, default]):
B := seq(plottools[line]([k, 0], [k, 0.4e-1]), k = 1 .. 4):
C := seq(plottools[line]([k, 0], [k, 0.2e-1]), k = 0 .. 4, 0.1):
T := plots[textplot]([seq([k, 0.3e-1, k], k = 0.5 .. 4, 0.5)], align = above): plots[display](A, B, C, T, scaling = constrained);

It is easy to derive a formula for the  n-th  member in this sequence:

f:=n->ceil((-1+sqrt(1+8*n))/2);

Examples:

f(2013);   f(10^10);

                 63
             141421

PS.  About  ceil  function see  ?trunc  and  http://en.wikipedia.org/wiki/Floor_and_ceiling_functions

Basic concepts, motivation, examples and references, see here

http://en.wikipedia.org/wiki/Group_representation

Use  curve=a*x  option.

Example:

with(CurveFitting):

LeastSquares([1, 2, 3], [1, 5/2, 3], x);  #  without the option

LeastSquares([1, 2, 3], [1, 5/2, 3], x, curve=a*x);  #  with the option

                                         1/6+x

                                       15/14*x