We describe an algorithm which computes the sign function of a complex matrix by using the continued fraction expansion of the inverse of the principal square root function at each step of the iteration. We show that the algorithm iteratively computes globally convergent main diagonal Pade approximants. The proposed algorithm avoids computing large matrix powers and performs fewer matrix inversions than Newton’s method. The algorithm is multiplication-rich and particularly suitable for implementation on vector and parallel computers. The stability analysis of the algorithm suggests that the errors introduced during a step are either suppressed or have limited effect on the next step. Finally, we summarize the results of our experiments on computing the sign function of certain matrices.