I’ll let the code speak for itself:

In [81]: m = Matrix(array([[1.0,1.0],[0.0,1.0]])) In [82]: def my_sqrt(x, num_iters): ....: r = 0.5*x ....: for i in xrange(num_iters): ....: r = 0.5*(r+x/r) ....: return r ....: In [83]: m*m Out[83]: array([[ 1., 2.], [ 0., 1.]]) In [84]: my_sqrt(m*m, 10) Out[84]: array([[ 1., 1.], [ 0., 1.]]) In [85]: m Out[85]: array([[ 1., 1.], [ 0., 1.]])

It’s always fun to see the math work out. At first when I learned that e^A for a matrix A may also be defined using the Taylor series in the usual way, I was really amazed. It still amazes me that this stuff works out so well. This is another kind of beauty.