Solution (source code)

= Solution

Multiplying the system by $Q^T$ turns it into $Rx=Q^Tb$. For nonsingular $A$, compute the transformed right side and solve the upper triangular system by <back substitution>. <Orthogonal> transformations preserve Euclidean lengths and do not amplify errors in that transformation.

If $A$ is singular,
$$
0=\det A=\det Q\det R=\det Q\prod_{i=1}^nr_{ii}.
$$
Since $\det Q=\pm1$, at least one diagonal entry of $R$ is zero. The resulting triangular system must then be checked for consistency and any free variables, rather than applying divisions by zero.