In a recent paper [1], Moler and Morrison have described an iterative algorithm for the computation of the Pythagorean sum

of two real numbers

and

without computing their squares or taking a square root. The subroutine is robust, short, portable, has a cubic rate of convergence, and is immune to floating-point overflows. In this note the method is extended to the efficient computation of the Pythagorean sum

of two real commuting matrices.