The procedure

expMod := proc(pol::polynom(integer,x),p::prime)
   map(y -> lcoeff(y)*x^(ldegree(y) mod p),pol):
end proc:

takes a polynomial in x with integer coefficients and prime exponents, and performs the requested operation. A couple of examples:

expMod(2*x^21 - 3*x^7 + 10*x^3 - 18,5);
expMod(2*x^21 - 3*x^7 + 10*x^3 - 18,17);

The task can be performed using

restart:
with(LinearAlgebra):
A := Matrix([
   [1,-1,0],
   [0, 1,0],
   [0, 0,1]
]):
B := Matrix([
   [ 1,0,2],
   [-1,1,0],
   [ 0,2,0]
]):
T := Transpose(LinearSolve(Transpose(A),Transpose(B))):
T . Vector([3,0,1]);

giving the following result:

The columns of the matrix A are the ones multiplied by T; the columns of the matrix B are the results of these multiplications. The use of Transpose in the penultimate line above is due to the following: Solving the equation T . A = B is equivalent to solving A^T T^T = B^T , where T means the transposed, and the latter equation can be solved for T^T using LinearSolve. [PS: Perhaps I am overlooking something in the help pages, but I do wonder why a command for solving X . A = B is not also available.]

Note added: But, coming to think of it, the result may be obtained simply as

3*Vector([1,-1,0]) + Vector([2,0,0]);

because 3*Vector([1,0,0]) + Vector([0,0,1]) is equal to Vector([3,0,1]), and T is a linear map.

Two explicit examples of construction (with three data points for the baseball), neither one using arrays (which does seem less than optimal in regards to this task):

• List of lists:

velocities := [[1,2,3],[2,3,4],[4,3,2]];

• List of Vectors:

velocities := [Vector([1,2,3]),Vector([2,3,4]),Vector([4,3,2])];

The mapping remove(i->L[i]<>3,L) fails for L := [4,5,6,3,7], for instance, because neither L[6] nor L[7] exist.

A working code is:

L := [-3,4,5,-6,3,7]:
remove(x -> evalb(x <> 3),L);

or, equivalently,

L := [-3,4,5,-6,3,7]:
select(x -> evalb(x = 3),L);

Another way is as follows:

m := Matrix(14):
m[1..7 ,1..7 ] := a:
m[1..7 ,8..14] := b:
m[8..14,1..7 ] := c:
m[8..14,8..14] := d:

Maybe you could use multidimensional Arrays. An example: Define a 2 x 3 x 4 3-D Array A by:

A := Array(1..2,1..3,1..4,(i,j,k) -> i*j*k);

Define also the following three Arrays which are different permuted versions of A (the first one, though, being A itself, but being included for expositional reasons):

A1 := Array(1..2,1..3,1..4,(i,j,k) -> A[i,j,k]);
A2 := Array(1..3,1..4,1..2,(i,j,k) -> A[k,i,j]);
A3 := Array(1..4,1..2,1..3,(i,j,k) -> A[j,k,i]);

Then, study the output of the following code which mimic your 'layering' (depth level) along different 'axes' of the 3-D Array:

A1[1],A1[2];
A2[1],A2[2],A2[3];
A3[1],A3[2],A3[3],A3[4];

The association may also be obtained as follows: Define the procedure

assoc := (n::posint) -> table([seq(seq([i,j],j = i+1..n),i = 1..n)]):

Then, for instance,

t := assoc(4):
t[1],t[2],t[3],t[4],t[5],t[6];
         [1, 2], [1, 3], [1, 4], [2, 3], [2, 4], [3, 4]

You can also define them as follows:

e := Vector(4,i -> Vector(4,j -> `if`(i = j,1,0)));

Error discovered after posting - I apologize.

You might do something like the following: Define a procedure expandDiff by

expandDiff := proc(expr,var)
	local a,b;
	if   type(expr,`+`) then map(expandDiff,expr,var)
	elif type(expr,`*`) then
		a := op(1,expr);
		b := `*`(op([2..-1],expr));
		expandDiff(a,var)*b + expand(a*expandDiff(b,var))
	elif type(expr,`^`) then
		a := op(1,expr);
		b := op(2,expr);
		expr*(expandDiff(b,var)*ln(a) + b/a*expandDiff(a,var))
	else Diff(expr,var)
	end if
end proc:

where for completeness functionality corresponding to the type `^` has been included. Then, try the following examples:

expandDiff(u+v,t);
expandDiff(u+v+w,t);
expandDiff(u*v,t);
expandDiff(u*v*w,t);
expandDiff(u^v,t);

Do you mean that the line should pass through two specific points p and q of the plane? If affirmative, then you can do something like the following:

p := [p1,p2]:
q := [q1,q2]:
y := (q[2] - p[2])/(q[1] - p[1])*(x - p[1]) + p[2];

Concerning item 1, you might do something like:

y[10] := x[10]+x[1]^2+x[2]*x[5]+x[3]*x[8];
f := unapply(
	subs({seq(x[i] = x_||i,i = 1..10)},y[10]),
	[seq(x_||i,i = 1..10)]
);
f(1,2,3,4,5,6,7,8,9,10);

Concerning item 2, you might do something like (where I have invented two extra polynomials y[2] and y[3] by just cyclically permuting the indices of the above y[10], now called y[1]):

y[1] := x[10]+x[1]^2+x[2]*x[5]+x[3]*x[8];
y[2] := x[1]+x[2]^2+x[5]*x[3]+x[8]*x[10];
y[3] := x[2]+x[5]^2+x[3]*x[8]+x[10]*x[1];
f := unapply(
	subs({seq(x[i] = x_||i,i = 1..10)},Vector([seq(y[i],i = 1..3)])),
	[seq(x_||i,i = 1..10)]
);
f(1,2,3,4,5,6,7,8,9,10);

Are you sure that it did not read something like arctan(a,b)? See the help page ?invtrig.

You can do something like the following: Define

x := (r,theta,phi) -> r*sin(theta)*cos(phi):
y := (r,theta,phi) -> r*sin(theta)*sin(phi):
z := (r,theta,phi) -> r*cos(theta):

Then (note that r1 and r2, and phi1 and phi2 are adjustable parameters):

r1,r2 := 1,3:
phi1,phi2 := 0,Pi:
spheres := plot3d([
	[x(r1,theta,phi),y(r1,theta,phi),z(r1,theta,phi)],
	[x(r2,theta,phi),y(r2,theta,phi),z(r2,theta,phi)]
],theta = 0..Pi,phi = phi1..phi2,grid = [20,20],color = [blue,red]):
planes := plot3d([
	[x(r,theta,phi1),y(r,theta,phi1),z(r,theta,phi1)],
	[x(r,theta,phi2),y(r,theta,phi2),z(r,theta,phi2)]
],r = r1..r2,theta = 0..Pi,grid = [10,20],color = green):
plots[display]([spheres,planes],
	scaling = constrained,
	orientation = [-30,70]
);

PS: The spheres part may be written more concisely as

spheres := plot3d([r1,r2],
	theta = 0..Pi,phi = phi1..phi2,
	coords = spherical,color = [blue,red]
):

but unless the planes part can be written similarly consisely I prefer the above code because it makes the parametrizations explicit.

You can do something like the following:

y1 := (x,z) -> x^2 + z^2:
y2 := (x,z) -> 1 - z^2:
plot3d(
	[y1(x,z),y2(x,z)],x = -1.2..1.2,z = -1.2..1.2,
	orientation = [40,50],labels = [x,z,y],axes = BOXED
);

The projection to the x-z-plane of the intersection curve of the two surfaces is given by the equation y1(x,z) = y2(x,z) which may be rearranged as x^2 + 2*z^2 = 1 which is the equation of an ellipse with major axis a = 1 along the x-axis, and minor axis b = 1/sqrt(2) along the z-axis, visually consistent with the above plot.