ENH: Add Tck support for tensor products and multiple dimension values.

The setup is that tensored b-splines map a polytope in E^m into the
vector space R^n. Note that the vector space may be an array with
multiple dimensions, i.e., an arrays with shape (2, 3) are taken to
be a vector space of dimension 6.

To support this functionality, the degrees of the b-splines that will be
tensored are kept in a list and the corresponding knot points in the
form of ndarrays are kept in another list. The coeffients are kept in a
single ndarray. The first m dimensions of the coefficient array
correspond to the m b-splines that will be tensored, and the remaining
dimensions are for the range space. As mentioned above, the range space
can be multi-dimensional and as a vector space its dimension is the
product of the remaining dimensions.

This can be illustrated in two tensored dimensions as

    f(x,y) = \sum_{j,m} c[j,m,...]*N_{i,j}(x)*N_{l,m}(y)

where f(x,y) is the b-spline value at (x, y), c is the coefficient array
and N_{i,j} is the j'th b-spline of degree i.

degrees into a list of degrees and knot points into a list of knot
points. The coefficient array has the tensor dimensions in the first
shape elements, and the vector
ENH: Add addition and subtraction operators to Tck object.

There are two cases here:

    1) Both objects are Tck objects
    2) One of the object is array_like

In the first case, both Tck objects must have the same type, knot
points. If n is the number of leading indices in the coefficient shapes
that index over the domain, then ones are inserted into the shape of the
smaller coefficient array starting at the n'th position until ther
arrays have the same number of dimensions. The effect of this is that
the vector values of the splines broadcast in the same way as standard
numpy ndarrays.

In the second case, the array_like object is first converted to an array
of the same dtype as the existing coefficient array and, if it has equal
or fewer dimensions than the range dimensions of the Tck object and
broadcasts to those coefficients, then it is added/subtracted. The
reason is isn't simply added without the dimension check is that some of
the Tck domain dimensions can be legitimately 1 and we do not want the
domain dimensions to change. Effectively, the array being added is a
scalar.
ENH: Add multiplication/division operators to Tck objects.

The valid types for these operations are 'array scalars', Tck objects do
not support multiplication/division with other Tck object. An array
scalar in this context is an array_like object that broadcasts over
coefficient array indices devoted to the domain dimensions. For example,
if two splines are tensored to give two spacial dimensions, then the
first two indexes of the coefficient array are the domain indexes and
the broadcast array must have at least two fewer indexes then the
coefficient array. We are essentially using the same array (vector) at
every knot point of the domain, it is like a constant function.
ENH: Make functions and tests tensor/multidimensional compatible.

The changes to Tck to support tensor product splines with vector
values made Tck incompatible with the previous functions in bsplines.py
that made use of the Tck type. The fixes those functions the test that go
with them.

The biggest change that may be notices by the user is that when a bspline
is evaluated it always returns a Tck object. When the tensored b-splines
are completely reduced the domain dimensionality becomes zero and the
results are in the coefficient array. It might be desireable to simply
return the results in that case, but that is for another day.

Tests pass, but this part of the work is incomplete and needs much more
testing.