Solving linear systems in Maxima

Question

I'm trying to write a general solver for a linear system in Maxima using linsolve(eqlist, varlist), but without having to explicitly specify the dimension of the problem.

This works, but fixes the dimension to 3:

linsolve( [ eq[0],eq[1],eq[2] ], [ a[0],a[1],a[2] ])

This doesn't:

solution(p):=(  
  array(eq,p+1), /* creating arrays of length p+1 */
  array(a,p+1),

  for i:0 thru p do (
    eq[i]: sum(binom(j+1,i)*a[j],j,i,p) = binom(p,i)
  ),

  linsolve(eq,a)
)

Any insight on how to get this to work?

---

Background behind the problem: this linear system arises when solving the finite summation of integer powers, i.e. the sum of finitely many squares, cubes, or general powers p. Although the finite sum of squares is straightforward, the general solution is surprisingly complicated: a discussion can be found here: Finite Summation by Recurrence Relations, Part 2.

Answer 1

Apparently in Maxima, lists and arrays are NOT the same underlying object.

Arrays are more complex and a bit of mess to get working (as suggested in this posting to the Maxima mailing list).

The problem goes away if we stay away from arrays and work with lists instead:

solution(p):= block([a, eq],        /* give subroutine variables local scope */
    v : makelist(a[i], i, 0, p),    /* create list of unknowns (0-indexed) */
   eq : makelist(sum(binom(j+1,i)*a[j],j,i,p) = binom(p,i), i, 0, p),  
                                    /* create list of equations (0-indexed) */
   linsolve(eq, v)
)