Your problem is easily solved by exhaustive search. Denoted by  a, b, c, d, e, f, g, h  vertices of the cube:

Code:

L:=combinat[permute]({seq(i^2, i=1..8)}):

n:=infinity: S:={}:

for i in L do

a, b, c, d, e, f, g, h := i[1], i[2], i[3], i[4], i[5], i[6], i[7], i[8]:

m := a*b+b*c+c*d+a*d+a*e+b*f+c*g+d*h+e*f+f*g+g*h+e*h: S:=S union {m}:

if m<n then n:=m:

M:=[m, ["a"=a, "b"=b,"c"= c, "d"=d, "e"=e, "f"=f, "g"=g, "h"=h]]:  fi:

od:

S:  # The set of all possible values of this sum 

nops(S);

max(S);

min(S);

M; 

622

9420

3586 

[3586, ["a" = 1, "b" = 36, "c" = 9, "d" = 64, "e" = 49, "f" = 16, "g" = 25, "h" = 4]]

For details see  ?rtable_indexing

f := 1/3 + x/5:

numer(f)/convert(denom(f), symbol);

There are a few extra parentheses and if..end if also not necessary:

aprcos:=proc(x,n::nonnegint)

local resp, i;

resp:=0;

for i from 0 to n do

resp:=resp+(-1)^i*x^(2*i)/(2*i)!;

end do;

resp;

end proc;

c:=Array([1, 3, 5, 7, 8, 2, 1]):

ArrayNumElems(c);

             7

We have a Diophantine equation  4^27+4^1016+4^n=m^2. Rewrite it as follows

4^1016+4^27=(m-2^n)*(m+2^n)

 Hence we find that  2^n<4^1016+4^27 .  Should be n <2033, as

is(2^2033>4^1016+4^27);

               true

So it is enough to make exhaustive search for  n  from 1 to 2032:

L:=[]:

for n from 1 to 2032 do

m:=sqrt(4^27+4^1016+4^n):

if type(m, integer) then L:=[op(L), n]:  fi:

od:

L; 

         [522, 2004]

plot3d(5, x = -3*Pi .. 3*Pi, y = (1/2)*sin(x) .. (1/2)*sin(x)+4, color = white, filled = true, scaling = constrained, grid = [50, 11], orientation = [93, -14, -23], axes = none, lightmodel = light1);

plot3d(5, x = -3*Pi .. 3*Pi, y = (1/2)*sin(x) .. (1/2)*sin(x)+4, color = white, filled = true, scaling = constrained, numpoints = 2500, orientation = [93, -14, -23], axes = none, lightmodel = light1);

Unfortunately, Maple cannot cope with this task:

assume(m, posint, n, posint):

is((2*m-1)^(2^n)-1 mod 2^(n+2) =0);

                             FAIL

But the problem is easily solved manually by induction, if we use the identity

(2*m-1)^(2^(n+1))-1=( (2*m-1)^(2^n)-1 )*((2*m-1)^(2^n)+1)

Change in the normal curve   phi(x)=1/(sigma*sqrt(2*Pi))*exp(-(x-a)^2/(2*sigma^2))  as the parameter  a  is changing from  -2  to  3  (sigma=1) :

plots[animate](plot,[1/sqrt(2*Pi)*exp(-(x-a)^2/2),x=-5..6,thickness=2],a=[seq(-2+0.02*k,k=0..250)]);

Change in the normal curve   phi(x)=1/(sigma*sqrt(2*Pi))*exp(-(x-a)^2/(2*sigma^2))  as the parameter  sigma  is changing from  0.6  to  2  (a=1) :

plots[animate](plot,[1/sigma/sqrt(2*Pi)*exp(-(x-2)^2/2/sigma^2),x=-2..6,thickness=3],sigma=[ seq(0.6+k*0.005,k=0..280)]);

See   PathInt , LineInt and SurfaceInt commands in  Student[VectorCalculus]  Subpackage.

Found an error in my previous solution. Offer exact solution found by the exhaustive search, pre considerably narrowed the search area.

Rewrite the system as follows {a+b=7.11-c-d, a*b=7.11/c/d} . Therefore, for  the specified values ​​of  c  and  d , the values  a  and  b  satisfy the quadratic equation  z^2+(c+d-7,11)*z+7,11/c/d=0 . Discriminant is equal to  (c+d-7.11)^2-4*7.11/c/d >=0 .

Build the domain:

plots[implicitplot]([c*d*(c+d-7.11)^2-4*7.11, c+d=7.11], color=[red, blue],c=0..7.11, d=0..7.11, numpoints=10000);

It is clearly seen that each variable must be in the range  0.7..3.2

restart;

R:=[infinity, []]:

for a from 0.7 by 0.01 to 7.11/4 do

for b from a by 0.01 to (7.11-a)/3 do

for c from b by 0.01 to (7.11-a-b)/2 do

for d from c by 0.01 to min(3.2,7.11-a-b-c) do

delta:=(7.11-(a+b+c+d))^2+(7.11-a*b*c*d)^2;

if delta

od: od: od: od:

                         [0., [1.20, 1.25, 1.50, 3.16]]

R;

It is better to write the procedure as follows. Then the number of factors can be arbitrary.

myproc1:=proc(L)   #  L is a list

convert(L, `*`);

end proc;

Solution of your problem:

map(myproc1, [[3,4], [2,3,4,10], [5,2,1]]);

                      [12, 240, 10]

fa:=proc(x)

local  i, j;

for i while x[i]=0 do; od;

for j from -1 by -1 while x[j]=0 do; od;

[i, ArrayNumElems(x)+j+1];

end proc; 

Because your system contains two redundant unknowns is easy to check that it has infinitely many solutions in four-dimensional space. But you do not look for solutions in the whole space, and on a discrete finite set  a=0.01 ... 7.11 by 0.01  and so on. I think that there does not exist exact solutions on this set.

 Natural formulation of the problem: on a given set of numbers to find a quartet  [a, b, c, d] , for which the left-hand sides of the system have the least differencies from the right-hand sides. As a measure of the difference consider  (7.11-(a+b+c+d))^2+(7.11-a*b*c*d)^2 . The procedure consists of two stages. To speed up work on the first stage we consider the wider step=0.1, and step=0.01 on the second stage near the optimum point of the  first stage.

restart:

R:=[infinity, []]:  # First stage

for a from 0.1 by 0.1 to 7.11 do

for b from a by 0.1 to 7.11 do

for c from b by 0.1 to 7.11 do

for d from c by 0.1 to 7.11 do

delta:=(7.11-(a+b+c+d))^2+(7.11-a*b*c*d)^2;

if delta<0.01 and delta<R[1] then R:=[delta, [a, b, c, d]]: fi:

od: od: od: od:

R;

a0:=R[2,1]: b0:=R[2,2]: c0:=R[2,3]: d0:=R[2,4]:  # Second stage

for a from a0-0.09 by 0.01 to a0+0.09 do

for b from b0-0.09 by 0.01 to b0+0.09 do

for c from c0-0.09 by 0.01 to c0+0.09 do

for d from d0-0.09 by 0.01 to d0+0.09 do

delta:=(7.11-(a+b+c+d))^2+(7.11-a*b*c*d)^2;

if delta<R[1] then R:=[delta, [a, b, c, d]]: fi:

od: od: od: od:

R;

         [.11296e-3, [.8, 1.8, 1.9, 2.6]]  # Result of the first stage

  [.266342400e-7, [.80, 1.76, 1.92, 2.63]]  # Result of the second stage

Verification:

7.11-convert([.80, 1.76, 1.92, 2.63],`+`);

7.11-convert([.80, 1.76, 1.92, 2.63],`*`);

                                 0.

                         0.00016320