MuPAD evaluation of local variables, Sum of row of array

Question

I have discovered a strange behavior in MuPAD version 5.7.0 (MATLAB R2011b), and I would like to know whether this is a bug, and if not so, what I am doing wrong. Ideally, I would also like to know why MuPAD does what it does.

Consider the array C of size 3x3, the elements of which have some example values. I would like to regard this array as an array of arrays and thus use cascaded indexing.

The problem apparently arises when both indices are local variables of different nested scopes, namely when the first index's scope is wider than the second index's one. There is no problem if the first index is a constant.

When I enter:

reset();
C := [[a,b,c],[d,e,f],[g,h,i]];
sum((C[3])[t], t = 1..3);
S := j -> sum((C[j])[t], t = 1..3);
S(3);

I get the following result:

I would expect lines 3 and 5 in the code (2 and 4 in the output) to yield the same result: g+h+i. Instead, line 5 produces a+e+i, which seems to be the diagonal of C.

When I do the same with product instead of sum, the result is even stranger, but might reveal more about the "error's" source, particularly DOM_VAR(0,2):

reset();
C := [[a,b,c],[d,e,f],[g,h,i]];
product((C[3])[t], t = 1..3);
eval(product((C[3])[t], t = 1..3));
S := j -> product((C[j])[t], t = 1..3);
S(3);
eval(S(3));

I get:

I might be on the wrong track here, but I suspect that a closure is created that tried to save the surrounding scope's local variables, which are undetermined at the time of the closure's creation. Also, substitutions seem to stop at some point, which is overridden by eval().

Practical Problem

The practical problem I am trying to solve is the following:

reset();
aVec := Symbol::accentUnderBar(Symbol::alpha)

Problem: Calculate multinomials of the form

hold(sum(x(i), i=i_0..i_k)^n)

On Wikipedia, the following form is defined:

sumf := freeze(sum):
hold(sum(x[i], i=1..m)^n)=sumf(binomial(n,aVec)*product(x[t]^aVec[t], t = 1..m), abs(aVec)=n);

In order to implement this, we need to define the set of vectors alpha, the sum of which equals m. These correspond to the set of possible compositions of n with length m and possible zero elements:

C := (n,m) -> combinat::compositions(n, MinPart = 0, Length = m)

For example, the sum

n := 3:
m := 3:
sumf(x[i], i=1..m)^n = sum(x[i], i=1..m)^n;

would call for these combinations of powers, each of which is one vector alpha:

A := C(n,m)

Additionally, we need the multinomial coefficients. Each such coefficient depends on vector alpha and the power n:

multinomial := (n, aVec) -> fact(n) / product(fact(aVec[k]), k = 1..nops(aVec))

For example, the number of times that the second composition appears, is:

multinomial(n, A[2])

Summation over all compositions yields:

sum(multinomial(n,A[i])*product(x[t]^A[i][t], t = 1..m), i = 1..nops(A))

The powers are correct, but the coefficients are not. This seems related to the boiled-down abstract problem stated first in this question.

Answer 1

I am posting this answer only to provide a working example of the practical problem from the question. horchler provided the decisive solution by proposing to use a matrix instead of a list of lists.

Basically, this is a modified transscript of the practical problem.

Practical Solution

The practical problem I am trying to solve is the following:

reset();
aVec := Symbol::accentUnderBar(Symbol::alpha)

Problem: Calculate multinomials of the form

hold(sum(x(i), i=i_0..i_k)^n)

On Wikipedia, the following form is defined:

sumf := freeze(sum):
hold(sum(x[i], i=1..m)^n)=sumf(binomial(n,aVec)*product(x[t]^aVec[t], t = 1..m), abs(aVec)=n);

In order to implement this, we need to define the set of vectors alpha, the sum of which equals m. These correspond to the set of possible compositions of n with length m and possible zero elements:

C := (n,m) -> combinat::compositions(n, MinPart = 0, Length = m)

For example, the sum

n := 3:
m := 3:
sumf(x[i], i=1..m)^n = sum(x[i], i=1..m)^n;

would call for these combinations of powers, each of which is one vector alpha:

A := matrix(nops(C(n,m)),m,C(n,m));

Additionally, we need the multinomial coefficients. Each such coefficient depends on vector alpha and the power n:

multinomial := (n, aVec) -> fact(n) / product(fact(aVec[k]), k = 1..nops(aVec))

For example, the number of times that the second composition appears, is:

multinomial(n,A[2,1..m])

Summation over all compositions yields:

sum(multinomial(n,A[i,1..m])*product(x[t]^A[i,t], t = 1..m), i = 1..nops(C(n,m)));

Finally, prove that the result transforms back:

simplify(%)

Answer 2

The expression [[a,b,c],[d,e,f],[g,h,i]] is not an array in MuPAD. It is a "list of lists." I'm guessing that's not what you're after. Lists are commonly used to initialize arrays and matrices and other objects (more here). For these examples, either arrays or matrices will work, but note that these two data types have different advantages.

Using array:

reset();
C := array([[a,b,c],[d,e,f],[g,h,i]]);
sum(C[3, t],t=1..3);
S := j -> sum(C[j, t], t = 1..3);
S(3);

which returns

          

Note the different way row/column indexing is represented relative to that in your question. Similarly, modifying your other example

reset();
C := matrix([[a,b,c],[d,e,f],[g,h,i]]);
product(C[3, t], t = 1..3);
S := j -> product(C[j, t], t = 1..3);
S(3);

results in

          

---

If you do happen to want to use lists for this, you can do so like this

reset();
C := [[a,b,c],[d,e,f],[g,h,i]];
_plus(op(C[3]));
S := j -> _plus(op(C[j]));
S(3);

which returns

          

The _plus is the functional form of + and op extracts each element from the list. There are other ways to do this, but this is one of the simplest.