We derive here an algorithm for a complete square root implementation of the modified Bryson-Frazier (MBF) smoother. The MBF algorithm computes the smoothed covariance as the difference of two symmetric matrices. Numerical errors in this differencing can result in the covariance matrix not being positive semi-definite. Earlier algorithms implemented the computation of intermediate quantities in square root form but still computed the smoothed covariance as the difference of two matrices. We show how to compute the square root of the smoothed covariance by solving an equation in the form CCT = AA(T) - BBT using QR decomposition with hyperbolic Householder transformations.