BUG: return one covariance matrix per `rhs` in `polyfit`
Addresses #24073
Taking the first residual from `resids` means that only the first
set of coefficients would get a covariance matrix.
This moves the line to make a `(1,)` shape array into an `int`
to the branch where there is only one `rhs`.
Fixes typos in `polyfit` docstring
fix int/ shape (1,) array test failures
add failing test with multiple `rhs`
Make `weights` of shape `(M,)` in `polyfit` test
Co-authored-by: Dan Foreman-Mackey <foreman.mackey@gmail.com>
Fix rank promotion problems for `polyfit`
Move `fac` expansion out to any size `y_arr`
Fix last `fac` broadcasting test case