@dharr  But if we use the  is  command to check for identity, we don't get the desired result:

is(sqrt(x+2*sqrt(x-1))=sqrt(x-1)+1);
is(sqrt(x+2*sqrt(x-1))=sqrt(x-1)+1) assuming x>=1 ;

When trying to calculate this sum symbolically, we get an incorrect result (returns  0 ):

restart;
sum(sum(`if`(igcd(m,n)=1, 1/(m^2*n^2), 0), m=1..infinity), n=1..infinity);

                                   0

@nm  Maple probably thinks that  B  is no simpler than  A . 
It's also worth noting that  expand  is a standard command in Maple for breaking an algebraic fraction (with a sum in the numerator) into the sum of several fractions:

expand((a+b+c)/d);
expand((a+b+c)/(d+e));

@nm  I think OP himself should clarify which option he had in mind.

You can use the procedure  Partition  from my post  https://mapleprimes.com/posts/200677-Partition-Of-An-Integer-With-Restrictions . This procedure works directly with the branch and bound method. It finds all partitions of a positive integer  n  with various restrictions. Below is an example of how the procedure works. We find all partitions of the number 100, where the number of partition members will be from  1  to  5 , and the members themselves are greater than 1.

L:=Partition(100, 1..5, 2):
nops(L);
L[1..100][];
       

@Ronan Thanks for the explanation, now I understand your solution. Vote up. Just note that this solution contains a logical gap. In fact, you have proven that if such a configuration as in the original figure exists, then the desired area is 23.

In Maple, the number pi should be encoded as Pi, not pi.

exp(I*Pi);

   #   -1

In my opinion, this whole check is so obvious and is done manually in a few seconds (using periodicity of sinus)  that it seems a bit strange to involve Maple here. 

@jalal  Use  Matrix(2,2, ... )  instead of  Vector[row]

@SSMB I understood your original question as simply replacing the symbols in the f[3] expression with some expressions you specified. I am completely unfamiliar with the specific issues involved in your task. I can't help you any further.

@mmcdara  The variant with

n->Matrix(n,(i,j)->a[i,j])

is more universal, I often use it and therefore remembered it well. It is suitable when instead of  a[i, j]  stands any expression , depending on  i  and  j .

@jalal 

restart;
F:=proc(t)
local P1, P2;
uses plots;
P1:=plot3d([y^2*cos(alpha),y,y^2*sin(alpha)],y=1..2,alpha=Pi/2..t, color="LightGreen");
P2:=spacecurve([y^2*cos(t),y,y^2*sin(t)],y=1..2, color=red, thickness=3);
display(P1,P2);
end proc: 
P:=plots:-spacecurve([0,t,t^2], t=-3..3, color=red, thickness=3):
P1:=plot([[0,0]], x=-4..4,y=-3..3);
f:=plottools:-transform((x, y)->[x, y, 0]):
plots:-animate(F,[t], t=Pi/2..2*Pi+Pi/2, frames=90, labels=[x,y,z], background=P, paraminfo=false, axes=normal, axis[3]=[color=white],view=[-4..4,-3..3,-4..9]);

@jalal  The following works :

restart;
F := t->plot3d([(y-1)*(y+2)*y*cos(alpha), y, (y-1)*(y+2)*y*sin(alpha)], y = -2.5 .. 0.5, alpha = (1/2)*Pi .. t, color = "LightBlue"):
P := plots:-spacecurve([0, t, (t-1)*(t+2)*t], t = -3 .. 2, color = red, thickness = 2):
plots:-animate(F, [t], t = (1/2)*Pi .. 2*Pi+(1/2)*Pi, frames = 90, labels = [x, y, z], background = P, paraminfo = false, axes = normal, view = [-4 .. 4, -3 .. 3, -4 .. 9], orientation=[-20,65]);

@nm  For problems with a parameter in Mathematica there is a very powerful  Reduce  command that does a full analysis of an equation or inequality:

I think developers should pay attention to this post, because none of the natural approaches give the desired result (code in Maple 2018.2):

solve( x^2+(y-2)^2=0, {x,y}) assuming real;  # Incorrect output
RealDomain:-solve( x^2+(y-2)^2=0, {x,y});  # NULL