import numpy

# find a, b, c such that y = ax^2 + by + c
# passes through the points (8,2), (15,3), and (6,3):
A = numpy.matrix([[64,8,1],[225,15,1],[36,6,1]])
B = numpy.matrix([[2],[3],[3]])
A**-1 * B

# equivalent to A**-1 * B,
# but faster when the matrices get big:
numpy.linalg.solve(A,B)

# easier way to build our Vandermonde matrix A:
A = numpy.vander([8,15,6])