Last semester, I completed my first Linear Algebra course and I am currently trying to make a 3d rendering engine from scratch. For this program, I decided to create a matrix math library. I've implemented a few methods that can solve the linear system of equations Ax=b where A is an n x n matrix. One of them uses Cramer's Rule and the other uses the following method: $x=A^{-1}b$ where $A^{-1}$ is calculated using the adjoint (or adjugate) of A. However, I would like to implement a method that can calculate x where A is an m x n matrix where $m ≠ n$. I also need to consider the case that the solution may be a span.
Note #1: I found a similar question titled "solve linear system Ax=b, when A is m.n matrix where m<n". One of the answers was (written by Robert Israel):
"If A has rank <m, there may be no solutions, depending on b. Assuming A has rank m, $AA^T$ is invertible and you can take $G=A^T(AA^T)^{-1}$" where $x=Gb$.
However, what if $\operatorname{rank}(A) ≠ m$? How do I find the solution in this case (if there is one)? Also, this method would be very slow in my opinion, especially for larger matrices. There is another solution written by Daryl, but I do not understand it.
Note #2: I know that there are iterative methods to solve linear systems of equations, but they are inaccurate (and I like to make calculations as accurate as possible). Also, they don't consider the fact that the solution may be a span.