The randomized Kaczmarz algorithm has received considerable attention recently because of its simplicity, speed, and the ability to approximately solve large-scale linear systems of equations. In this paper we propose randomized double and triple Kaczmarz algorithms to solve extended normal equations of the form \(\mathbf{A}^\top \mathbf{Ax}=\mathbf{A}^\top \mathbf{b}-\mathbf{c}\). The proposed algorithms avoid forming \(\mathbf{A}^\top \mathbf{A}\) explicitly and work for arbitrary \(\mathbf{A}\in \mathbb {R}^{m\times n}\) (full rank or rank-deficient, \(m\ge n\) or \(m<n\)). Tight upper bounds showing exponential convergence in the mean square sense of the proposed algorithms are presented and numerical experiments are given to illustrate the theoretical results.