First define some auxiliary function to test for whether or not some positive integer has a prime factor greater than 11:

primeFactorGT11 := (n::posint) ->
   evalb(select(`>`,map(x -> x[1],ifactors(n)[2]),11) <> []):

Then the following can test the conjecture up to some limit N:

N := 10000:
# Initialization
p1 := primeFactorGT11(12):
p2 := primeFactorGT11(13):
p3 := primeFactorGT11(14):
p4 := primeFactorGT11(15):
if not `or`(p1,p2,p3,p4) then error "Conjecture is false" end if;
# Looping
for n from 16 to N do
   p1 := p2;
   p2 := p3;
   p3 := p4;
   p4 := primeFactorGT11(n);
   if not `or`(p1,p2,p3,p4) then error "Conjecture is false" end if;
end do:

The cyclic reuse of p1,p2,p3,p4 inside the loop is to avoid having to factorize the same number four times.

A variation on the theme of tomleslie:

A[..nops(remove(`>`,ListTools:-PartialSums(A),50))];

Concerning the first error, you cannot use a summation variable (a dummy variable) that has previously been assigned a specific value, in this case n. If you remove the assignment n := 1 above the summation, then the error vanishes.

Concerning the second error, replace assignment := with equality =.

Let f(z) be the integrand:

f := (z) -> (2*z+1)/(z*(z-1)-5)^2;

It has two real-valued poles, each of second order, only one of which, though, lies within the contour:

poles := {solve(denom(f(z)) = 0)};                              # The two second-order poles
pole  := select(x -> evalf(x) > -2 and evalf(x) < 2,poles)[];   # The pole lying within the contour

The contour integral thus evaluates to

simplify(2*Pi*I*residue(f(z),z = pole));

This result may actually also be obtained by direct integration, without recourse to the residue theorem: Let

gamma1 := (t) -> 2*(+I-1)*t + 2;
gamma2 := (t) -> 2*(-1-I)*t + 2*I;
gamma3 := (t) -> 2*(-I+1)*t - 2*1;
gamma4 := (t) -> 2*(+1+I)*t - 2*I;

be the four sub-contours of the full contour, each with parameter t running from 0 to 1. Then the contour integral is given by

simplify(
   +int(f(gamma1(t))*diff(gamma1(t),t),t = 0..1)
   +int(f(gamma2(t))*diff(gamma2(t),t),t = 0..1)
   +int(f(gamma3(t))*diff(gamma3(t),t),t = 0..1)
   +int(f(gamma4(t))*diff(gamma4(t),t),t = 0..1)
);

as before.

map(c -> select(x -> coeffs(x) = c,p),{coeffs(p)});

PS: I assume that p is always in expanded form.

nops(sort~({seq(seq([a,b,36-a-b],b = floor((35-a)/2)+1..17),a = 1..17)}));

And generally for any perimeter P:

N := proc(P::posint)
   local p := floor((P-1)/2);
   nops(sort~({seq(seq([a,b,P-a-b],b = floor((P-a-1)/2)+1..p),a = 1..p)}));
end proc:

The output of the procedure is consistent with the formula by Alcuin:

seq(N(2*n)                 ,n = 1..20);
seq(round((2*n)^2/48)      ,n = 1..20);   # The formula by Alcuin for even P
seq(N(2*n+1)               ,n = 0..20);
seq(round(((2*n+1)+3)^2/48),n = 0..20);   # The formula by Alcuin for odd P

I would define the tensors myself using the Define command:

Define(S[mu,nu]     = 1/(d-2)*(Ricci[mu,nu] - 1/(2*d-2)*Ricci[~alpha,alpha]*g_[mu,nu])):
Define(C[mu,nu,rho] = D_[nu](S[rho,mu]) - D_[rho](S[nu,mu])):

S and C are then the Schouten tensor and Cotton tensor, respectively, following the formulas for them given at the help pages ?schouten and ?cotton. I call the dimension of spacetime d.

The Physics package calculates all the required quantities 'on the fly':

Setting up the geometry:

with(Physics):
Setup(
   Signature   = `-+++`,
   Coordinates = (X = (t,r,theta,phi)),
   metric = -exp(alpha(r))*(dt^2)+exp(beta(r))*(dr^2)+r^2*(dtheta^2)+r^2*(sin(theta)^2)*(dphi^2)
,quiet):
g_[matrix];   # Check the metric

Now, all quantities, like the Christoffel coefficients and the components of the Riemann tensor, are readily available, an examples being:

Riemann[ 1,2,1,2];
Riemann[~1,2,1,2];

You can access all nonzero components of, say, the Riemann tensor as

Riemann[nonzero];

Concerning contractions: The contraction of the Riemann tensor, say, being by definition the Ricci tensor, can be found either via the predefined components of Ricci, as Ricci[1,1], for instance, or by explicit contraction like follows [this example is of course a bit artificial as Ricci is already available, but it is just to illustrate the point]:

Riemann[~mu,1,mu,1];   # For a proper contraction, remember to have upper first index using ~
Ricci[1,1];   # Comparison

Here is another suggestion:

find := (m) -> map(x -> `if`(lhs(x) = m,rhs(x),NULL),convert(zip(`=`,M,N),list)):
find(1);
find(2);
find(3);

Note added: As Christopher correctly points out below, my code above does not work if matrix (lower case) is used. But it works if Matrix (upper case) is used, see my reply to his comment.

This is probably neither the smartest nor the fastest idea, but it seems to work:

ranges := [1..2,1..2,1..3]:   # An example of ranges provided by the user
T := combinat:-cartprod(map(x -> [$x],ranges)):
S := NULL:
while not T[finished] do
   S := S,F(T[nextvalue]()[])
end do:
S;

Update: The following much simpler construction suddenly occurred to me:

ranges := 1..3,1..3,1..4:
map(x -> F(x[]),[indices(Array(ranges))])[];

Concerning the bold part of your code, you may consider writing

not member(i,{k,l,m,o,q,r,s})

for the conditions on i, etc.

Is something along the following lines what you have in mind?

f := (m::posint,n::posint) -> Matrix(m,n,(i,j) ->
   expand(binomial(m,j)*binomial(n,i)*x^i*(1-x)^(m-i)*y^j*(1-y)^(n-j))
):

Example of use:

f(2,2);

Edit: It suddenly occurs to me that there is a slight problem with the above construction: it excludes the possibility of having i = 0 and/or j = 0. To include these as well, the construction could perhaps be changed to the following:

f := (m::posint,n::posint) -> Matrix(m+1,n+1,(i,j) ->
   expand(binomial(m,j-1)*binomial(n,i-1)*x^(i-1)*(1-x)^(m-i-1)*y^(j-1)*(1-y)^(n-j-1))
):

Note that the matrix has been enlarged with one row and one column, and that the row and column indices, i and j, respectively, have both been shifted by -1 in the formula. But perhaps that is not pretty?

I don't think that you can do much about these square roots. If I call your expressions for expr1 and expr2, respectively, and with the following code lines pick out the expressions underneath the square roots and try to factor each one of them, nothing happens:

factor(op([1,5,3,1],expr1));
factor(op([1,3,2,1],expr2));

If they factored into squares, say, then something could be done.

Your answer to Preben makes me think that what you are asking for is the number of nonzero elements in each row of the array. If that is indeed the case, then your could do as follows:

map(x -> nops(remove(`=`,x,0)),convert(A,listlist));

This gives you as a list the number of nonzero elements of the individual rows.

I was thinking something along the following lines:

with(Physics):
Setup(coordinatesystems = (X = [t,r,x,y]),quiet):
Lambda := (F) -> D_[mu](D_[~mu](F)):
# Unperturbed metric
U := Matrix(4,4,{(1,1) = -f(r),(2,2) = 1/f(r),(3,3) = r^2,(4,4) = r^2}):
Setup(metric = U):
SumOverRepeatedIndices(Lambda(Phi(x,y)));
# Perturbed metric
P := Matrix(4,4,{(1,3) = h(r),(2,3) = r^2*g(r),(3,1) = h(r),(3,2) = r^2*g(r)}):
Setup(metric = U + P):
expr := SumOverRepeatedIndices(Lambda(Phi(x,y)));

To perform (multi-dimensional) Taylor expansions of expr to various orders, substitutions back and forth between the perturbation functions g(r) and h(r) and their derivatives, and some symbols, s1,s2,s3,..., say, will be performed below:

INDS := remove(has,indets(expr,function),{Phi(x,y)} union indets(U,function));
SUBS := {seq(INDS[i] = s||i,i = 1..nops(INDS))};
SUBS_REVERSE := map(rhs = lhs,SUBS);
expr := eval(expr,SUBS):

First order Taylor expansion (it equals the unperturbed sum):

eval(mtaylor(expr,rhs~(SUBS),1),SUBS_REVERSE);

Second order Taylor expansion:

eval(mtaylor(expr,rhs~(SUBS),2),SUBS_REVERSE);

Third order Taylor expansion:

eval(mtaylor(expr,rhs~(SUBS),3),SUBS_REVERSE);