ENH: interpolate: piecewise N-D tensor-product polynomials

This PR adds piecewise tensor-product polynomials to scipy.interpolate.

This is partly WIP, the following features are still missing:

- definite integration
- check how many copies it makes --- at least antiderivative computation makes quite many

Constructing interpolating polynomials is also on the TODO list, but out of scope for this PR which just tries to put the data structures in place.

(The next step would then be to implement similar structures and evaluation routines as for PPoly/NdPPoly for the B-splines.)