Make cosmology distance functions act like ufuncs

* Change FLRW.comoving_distance and FLRW.absorption_distance
  to handle input arrays with shape.  This automatically propogates
  through to the other distance functions in FLRW.
* Add unit tests for this functionality
* This starts to address #1816, but not completely.
* Closes #3177
* There seems to be a 1-2% performance penalty for distance calculations
  caused by this change.