Efficient algorithms for the matrix cosine and sine

Several improvements are made to an algorithm of Higham and Smith for computing the matrix cosine. The original algorithm scales the matrix by a power of 2 to bring the ∞-norm to 1 or less, evaluates the [8/8] Padé approximant, then uses the double-angle formula cos (2A)=2cos 2A−I to recover the cosine of the original matrix. The first improvement is to phrase truncation error bounds in terms of ‖A2‖1/2 instead of the (no smaller and potentially much larger quantity) ‖A‖. The second is to choose the degree of the Padé approximant to minimize the computational cost subject to achieving a desired truncation error. A third improvement is to use an absolute, rather than relative, error criterion in the choice of Padé approximant; this allows the use of higher degree approximants without worsening an a priori error bound. Our theory and experiments show that each of these modifications brings a reduction in computational cost. Moreover, because the modifications tend to reduce the number of double-angle steps they usually result in a more accurate computed cosine in floating point arithmetic. We also derive an algorithm for computing both cos (A) and sin (A), by adapting the ideas developed for the cosine and intertwining the cosine and sine double angle recurrences.