This is in the realm of Numerical Approximation, an important mathematical branch.
In Maple see ?numapprox
It is usually needed a careful study of the function to be approximated and an estimation of the error (absolute/relative).

I see two solutions:
1. Revise the automatic simplication for INTERVALs
2. Overload the `*` operator by the user (for the second issue).

To compute evalr(x*x) for
x:=INTERVAL(-1 .. 1);
one may use the ad-hoc
evalr('evalr'(x)*x);

BooleanSimplify returns a sum of products [sum=%or, product=%and].
Your expression is far from having a single term i.e. a product only. It is huge.
In order to obtain something you will have to fix some variables.

Download boolean.mw

Download projectile.mw

After studying Maple such things will be possible. E.g.

funny := module() option package;  export `<`;
   `<`:= proc(a::list(realcons), b::realcons) option overload;
              select(t -> t<b, a);
         end proc;
end module:

with(funny):
a:=[1,2,3,4,5]:
b:=a<4;
                         b := [1, 2, 3]

Suppose you want F[i, j,...](x, y,...)  to return  x^i + y^j + ...

F:=proc()
  local ind:=op(procname); 
  add( [args] ^~ [ind] )
end:

F[i,j,99](x,y,z);
    x^i+y^j+z^99

Of course some parameter checks in the procedure are recommended.

b:=A[..,6]:
C:=A[..,[1..5,7..11]]:
(C^(-1).b)[6];

(simplify@expand@convert)( abs(u+v)^2+abs(u-v)^2, conjugate );

1.   As documented.

2. Because is tries to be smart.
Also:
is(0.5 = 1/2);
     true

3. Because is is not perfect.
Also:
is(exp(1)+Pi, rational);
    false
Should be FAIL, the result is not known (yet)!

f:=exp(t^(1/3)-2*t^(7/4)):
s:=series(f,t,6);

sort(s,t,ascending);

The timing seems to be similar.

t:=time[real]():
L:=combinat:-permute([1,1,1,1,0,0]):
C:=Iterator:-CartesianProduct([seq(1..15)]$6):
n:=0:
T:=table():
for c in C do
if add(L[c[i]][1],i=1..6)<>4 or
   add(L[c[i]][2],i=1..6)<>4 or
   add(L[c[i]][3],i=1..6)<>4 or
   add(L[c[i]][4],i=1..6)<>4 or
   add(L[c[i]][5],i=1..6)<>4 or
   add(L[c[i]][6],i=1..6)<>4 then next fi;
n:=n+1;
T[n]:=Matrix([seq(L[c[i]],i=1..6)]);
od:
'num'=n, 'time'=time[real]()-t;

        num = 67950, time = 65.329

Edit. We gain about 2 seconds if one of the 6 tests inside if is removed.

@shimaa sadk 

The inner sum can be computed symbolically. Do it.
Fot the outer one use evalf(Sum(...,T=0..infinity);

Unfortunately evalf for GAMMA and also Zeta is VERY slow for large arguments. E.g.

evalf( GAMMA(1+10^20*I));

The evalf operations are known to be hard/impossible to interrupt.