I also just see that there is a function lagrange_polynomial in sage

e.g.

# using the definition of Lagrange interpolation polynomial

sage: R = PolynomialRing(QQ, 'x')

sage: R.lagrange_polynomial([(0,1),(2,2),(3,-2),(-4,9)])

I mean this should be super easy now. Just plug in your argument-value-pairs and you have the interpolating polynomial (no matrix fuzz).

Then you can apply this interpolating polynomial to non-real values.

Or extract the coefficients as you like.

However I didnt check how long it takes

For variants see

http://wiki.sagemath.org/sage-4.0.1

e.g.

# using the definition of Lagrange interpolation polynomial

sage: R = PolynomialRing(QQ, 'x')

sage: R.lagrange_polynomial([(0,1),(2,2),(3,-2),(-4,9)])

I mean this should be super easy now. Just plug in your argument-value-pairs and you have the interpolating polynomial (no matrix fuzz).

Then you can apply this interpolating polynomial to non-real values.

Or extract the coefficients as you like.

However I didnt check how long it takes

For variants see

http://wiki.sagemath.org/sage-4.0.1