Adds support for sparse A,B,C,D matrices in AffineSystem.

[RussTedrake]

Towards properly supporting MatrixGain for sparse matrix multiplicatoin of the actuation matrix in #20971.

This PR relocates sparse_and_dense_matrix from solvers into math with no changes (apart from the namespace), and uses it in AffineSystem.

---

This change is 

[jwnimmer-tri]

Do you have any (even manual, unreliable) benchmarks that using sparse instead of dense is faster in cases we care about (e.g., inverse dynamics)? There is a possibility that sparse matmul is a serious performance degradation for a meaningful number of user of this class, so at least we should be sure that there is any upside before we dig into that question. In mathematical programs with hundreds or thousands of variables, the sparsity is usually a win (especially as it helps the solver know that some variables are irrelevant). Here though, for a num_robot_dofs-sized matmul, it's not clear to me that the branchy sparse index chasing will be a net win. It would be a more convincing PR if you had a performance measurement you could post alongside it.

I also could use your help with the python bindings ...

Is this maybe the first time we have sparse matrix arguments with default values? That might be the problem. I think the default arguments get eagerly converted to Python objects during loading, even if the user never calls that overload. If so, you might need to change e.g. py::arg("D") = Eigen::SparseMatrix<double>(), to be py::arg("D") = Eigen::MatrixXd(), instead. Or maybe default them to None and then in the body of the binding convert nullptr to an empty C++ matrix.

[jwnimmer-tri]

One way around the perf question would be to remember whether the user did the dense constructor or sparse constructor, and then likewise the Calc function would then use either dense or sparse matrix operations. That way, most existing users would be unaffected.