@sarra   In order that the matrix was square, the number of equations should equal the number of unknowns. You have 16 unknowns but only four equations.

What do you mean by a transposition?  Transposition (accoding to Wiki) a permutation which exchanges two elements and keeps all others fixed. In this sense, for your example is enough 2 transpositions: [w,x,y,z]->[y,x,w,z]->[y,x,z,w]

evalf(Int(abs(exp((2*Pi*I)*x)+exp((2*Pi*I)*y)), [x = 0 .. 1, y = 0 .. 1]));

identify(%);

                                                                             4/Pi

@Carl Love 

A := Matrix([[1, 2, 3], [3, 1, 2], [2, 3, 1]]):
B := Matrix([[1, 2, 3], [3, 1, 2], [2, 3, 1]]):

is(A=B);

                             false

@Carl Love   I do not understand why we need to override the matrices. Maple tests  A  and  A1  as the identical matrices:

restart;

with(LinearAlgebra):

A := Matrix([[1, 2, 3], [3, 1, 2], [2, 3, 1]]);

A1 := Matrix(3, shape = Circulant[[1, 2, 3]]);

Equal(A, A1);

@Muhammad Usman   

applyrule(exp(x::anything)^y::anything=exp(x*y), expand(exp(2*a+b)));

@Markiyan Hirnyk   I do not know the causes of your problems. Here's a screenshot of the playing this animation in Power Point. Below is visible panel to control animation:

@Markiyan Hirnyk   I do not understand what's the problem? Click on the picture  and the bottom panel appears to start the animation.

Here is a sample ppt-file with the animation of the tangent line as the limiting position of the secant. Animation is made in Maple (M 12 classic), saved as a gif-file and inserted to Power Point in the above method.

Animation_tangent_line.pptx

@Markiyan Hirnyk   Everything works, I checked. On the instrument panel of PP click insert  (not paste)  and so on (see the order of actions above). After  All Files  click your  gif-file from documents.

@laporte bernard   I do not see any errors. To better see the elliptic paraboloid, you just need to extend the domain of definition for the graph of the function  f (x, y)  (the red line on the surface shows the steepest increase of the function  f(x,y) ):

restart;

f:=(x,y)->2+x^2 +1/4*(y^2);

Px, Py:=1, 2;

DPx, DPy:=D[1](f)(Px,Py), D[2](f)(Px,Py);  # Partial derivatives

maximize(DPx*cos(alpha)+DPy*sin(alpha), location);  # Maximum of directional derivative

A:=plot3d(f(x,y), x=-3..3, y=-sqrt(36-4*x^2)..sqrt(36-4*x^2), axes=normal):

B:=plots[arrow]([Px,Py,0],[DPx,DPy,0], color=red):

C1:=plots[spacecurve]( [1,2,t], t=0..4, color=black,linestyle=2):

C2:=plots[spacecurve]( [2*t+1,t+2,f(2*t+1,t+2)], t=0..0.8, color=red, thickness=3):

T:=plots[textplot3d]([[1,2.2,0,"P"], [1,2.2,4,"Q"]], align=right, font=[TIMES,BOLD,16], color=blue):

plots[display](A,B,C1,C2,T, orientation=[-15,65], scaling=constrained);

@Markiyan Hirnyk   Thanks for the clarification of the causes of issues.

@acer  Thanks for the detailed answer! 

@Kitonum  If you replace the fourth line of code with the following line, the procedure will run faster:

L:=[seq(ithprime(k) $ min(n, floor(a/ithprime(k))), k=1..pi(a))];

@fairuzalwani 

E := [0, 2, 7, 15, 26, 40]:

map(C, [3, 10, 24], [0, 2, 7, 15, 26, 40]);

                           [3, 4, 5]