@did 

You have no choice. A perfect CAS does not exist and never will.

P.S. Fourier series? Trivial result?

@Carl Love 

gamma is known by the assume facility. So its tables contain the information.

about(gamma);
gamma:
  is assumed to be: gamma

@Markiyan Hirnyk 

The question was: why reject the 0 answer and accept the other one without some math investigations?
I think that Carl's suggestion (+ some evalf testss) is more reliable.

BTW. See this:

restart: 
int( exp((1+a+b*I)*exp(I*x)) , x=-Pi..Pi)  assuming real;  # a  and  b  are assumed to be real
simplify(%)  assuming real;

    0
    0

@Kitonum 

So, we have 2 symbolic results. How can we decide which one is true.  Is it always true? (e.g. a=0 in your answer).
Hard to answer without some maths.

 

@tomleslie 

can be found with a few computations by hand using the inequalities  u-1 < floor(u) <= u.
They can be obtained also like this:

Download bounds.mw

@Markiyan Hirnyk 

I think that the code for the projection and the regions are enough for this problem.
Any other solid people can start from here.

@Markiyan Hirnyk 

Yes, the MMA expression seems to be the same (the 13 regions), but that triangle is a nonsense. I am not sure about the boundary points.
The MMA projection  is:

@Markiyan Hirnyk 

A projection of a triangle is a triangle. As can be seen from the first plot (with code) this is false. So, MMA is totally wrong.

@Markiyan Hirnyk 

The MMA plot does not seem to be correct.
Can you plot the (x1,x2) projection?

The 13 regions obtained with Maple are:

Region1
{442/7 < x2, x1 < 80, x2 < 64, 381/2-(7/4)*x2 < x1}
Region2
{432/7 < x2, x1 < 78, x2 < 62, -(7/4)*x2+186 < x1}
Region3
{414/7 < x2, x1 < 78, x2 < 60, -(7/4)*x2+363/2 < x1}
Region4
{438/7 < x2, x1 < 75, x2 < 442/7, -(7/4)*x2+369/2 < x1}, {442/7 <= x2, x2 <= 444/7, 74 < x1, x1 < 75}, {x1 = 74, 442/7 < x2, x2 < 64}, {74 < x1, 444/7 < x2, x1 < -(7/4)*x2+186, x2 < 64}, {x1 = -(7/4)*x2+186, 444/7 < x2, x2 < 64}
Region5
{426/7 < x2, x1 < 75, x2 < 430/7, -(7/4)*x2+363/2 < x1}, {430/7 <= x2, 74 < x1, x1 < 75, x2 < 62}, {x1 = 74, 430/7 < x2, x2 < 62}
Region6
{410/7 < x2, x1 < 76, x2 < 414/7, -(7/4)*x2+357/2 < x1}, {414/7 <= x2, 75 < x1, x1 < 76, x2 < 60}, {x1 = 75, 414/7 < x2, x2 < 60}
Region7
{404/7 < x2, x1 < 76, x2 < 58, -(7/4)*x2+177 < x1}
Region8
{408/7 < x2, x1 < 75, x2 < 412/7, -(7/4)*x2+177 < x1}, {412/7 <= x2, x2 <= 414/7, 74 < x1, x1 < 75}, {x1 = 74, x2 <= 418/7, 412/7 < x2}, {74 < x1, 414/7 < x2, x1 < -(7/4)*x2+357/2, x2 < 418/7}, {x1 = -(7/4)*x2+357/2, 414/7 < x2, x2 < 418/7}
Region9
{x2 = 60, 147/2 < x1, x1 < 74}, {60 < x2, x1 < 74, x2 < 426/7, -(7/4)*x2+357/2 < x1}, {426/7 <= x2, x2 <= 430/7, 72 < x1, x1 < 74}, {x1 = 72, 426/7 < x2, x2 < 62}, {72 < x1, 430/7 < x2, x1 < -(7/4)*x2+363/2, x2 < 62}, {x1 = -(7/4)*x2+363/2, 430/7 < x2, x2 < 62}
Region10
{58 <= x2, 72 <= x1, x2 <= 412/7, x1 < 74}, {x1 = 72, 412/7 < x2, x2 < 60}, {72 < x1, 412/7 < x2, x1 < -(7/4)*x2+177, x2 < 60}, {x1 = -(7/4)*x2+177, 412/7 < x2, x2 < 60}
Region11
{56 <= x2, x2 <= 402/7, 70 < x1, x1 < 72}, {x1 = 70, 56 < x2, x2 < 58}, {70 < x1, 402/7 < x2, x1 < 345/2-(7/4)*x2, x2 < 58}, {x1 = 345/2-(7/4)*x2, 402/7 < x2, x2 < 58}
Region12
{x2 = 54, 70 <= x1, x1 < 72}, {x1 = 70, 54 < x2, x2 < 55}, {54 < x2, 70 < x1, x1 < 333/2-(7/4)*x2, x2 < 55}, {x1 = 333/2-(7/4)*x2, 54 < x2, x2 < 55}
Region13
{x2 = 56, 68 <= x1, x1 <= 137/2}, {x1 = 68, x2 <= 394/7, 56 < x2}, {56 < x2, 68 < x1, x1 < 333/2-(7/4)*x2, x2 < 394/7}, {x1 = 333/2-(7/4)*x2, 56 < x2, x2 < 394/7}

The correct 3d plot:

@Preben Alsholm 

To obtain the 1D version I have pasted the Int in a 2D execution group and then lprint(%).
BTW, the exact and approx value of the integral are computed below.
 

@Carl Love 

OP's function written in a natural and simplified form is

f:=piecewise(0<-w*x*y+u*x and -w*x*y+u*x<1,1)*piecewise(0<-w*x*y+v*y and-w*x*y+v*y<1,1) ;

Now,

int(f,  [x=0..1,y=0..1,u=0..1,v=0..1,w=0..1]);
       7/12

Converting to Heaviside
g:=convert(f,Heaviside);
        g := Heaviside(-w*x*y+u*x)*Heaviside(-w*x*y+v*y)*Heaviside(w*x*y-u*x+1)*Heaviside(w*x*y-v*y+1)
int(g,  [x=0..1,y=0..1,u=0..1,v=0..1,w=0..1]);
    (interrupted)
but ...

int(g,  [u=0..1,v=0..1,w=0..1,x=0..1,y=0..1]);   # ok!
                               7/12

I wanted to check manually the result.
I was surprised to find that Maple is not able to solve the system

[0 < -w*x*y+u*x, -w*x*y+u*x < 1, 0 < -w*x*y+v*y, -w*x*y+v*y < 1, 0 < x, x < 1, 0 < y, y < 1, 0 < u, u < 1, 0 < v, v < 1, 0 < w, w < 1]

using solve, or SolveTools:-SemiAlgebraic

even if it's easy to see that the solution is:

0<x<1, 0<y<1, 0<w<1, w*y<u<1, w*x<v<1

so the integral is indeed
int((1-w*x)*(1-w*y), x=0..1, y=0..1, w=0..1) = 7/12.

@mmcdara 

It's about the user, not the laptop ... :-)

@mmcdara 

Have you tried the code?

discont was added just for a correct plot.

It seems that you actually do not want an answer. In this case you can simply ignore mine.

kernelopts(version);
   Maple 2018.1, X86 64 WINDOWS, Jun 8 2018, Build ID 1321769

bin.mw

@Klausklabauter 

So, LDLt is actually implemented in

Student:-NumericalAnalysis:-MatrixDecomposition(A, 'method'='LDLt')

Note however that it is too slow (compared with Cholesly for floats) for large matrices.

@Ahmed111 

For p:=1; q:=1; lambda:=2; i:=1  the system has only the trivial solution.
Maple says u1=u2=0  for i = I  too. If you want to check this, try to use the method from the provided worksheet.