Angle.to_string() is slow for large arrays (per-element string formatting)

While working on the `Time` string formatting speedup (#19929 / #19942) I noticed the coordinates side has basically the same problem: `Angle.to_string()` (and so `Longitude`/`Latitude`/`SkyCoord.to_string()`) gets really slow on large arrays, and it scales linearly with the number of elements.

Here's a quick benchmark on current `main`:

```python
import numpy as np, time
from astropy.coordinates import Angle, SkyCoord
import astropy.units as u

for n in [10_000, 100_000, 1_000_000]:
a = Angle(np.random.uniform(0, 360, n), unit=u.deg)
t0 = time.perf_counter(); a.to_string(sep="dms"); t1 = time.perf_counter()
print(f"Angle.to_string n={n:>9,}: {t1 - t0:.3f} s")

ra = np.random.uniform(0, 360, 1_000_000)
dec = np.random.uniform(-90, 90, 1_000_000)
c = SkyCoord(ra, dec, unit="deg")
t0 = time.perf_counter(); c.to_string("hmsdms"); t1 = time.perf_counter()
print(f"SkyCoord.to_string hmsdms n=1,000,000: {t1 - t0:.3f} s")
```

```
Angle.to_string n= 10,000: 0.043 s
Angle.to_string n= 100,000: 0.434 s
Angle.to_string n=1,000,000: 4.252 s
SkyCoord.to_string hmsdms n=1,000,000: 23.194 s
```

So a million-row catalogue takes ~4 s just to format one angle column, and ~23 s for an `hmsdms` SkyCoord. That's painful when you're saving coordinates to a text table or printing them.

The interesting part is that the actual math is already fast. `_decimal_to_sexagesimal` runs on the whole array with numpy and does a million elements in about 0.015 s. All the time is going into building the strings one element at a time. In `astropy/coordinates/angles/core.py`, `to_string` does:

```python
format_ufunc = np.vectorize(do_format, otypes=["U"])
result = format_ufunc(self.to_value(unit))
```

`np.vectorize` is just a Python loop, and `do_format` calls `_decimal_to_sexagesimal_string` per element, which assembles the string with scalar f-strings (the rounding, the `02d` minutes, the padding, etc.). So every element pays the full Python overhead even though the numbers are already computed in bulk.

This is the same shape of problem we just fixed for `Time` in #19942, so I'd like to do the equivalent here: build the output with array operations so the string assembly runs at numpy speed instead of in a per-element loop. It needs to handle the same options the current code does (`sep`, `pad`, `precision`, `fields`, `alwayssign`, and the `latex`/`unicode` separators).

One thing I wanted to check before starting: would you prefer a pure-numpy vectorized version, or a small C formatter along the lines of the existing fast C parser in `astropy.time`? The numpy route is simpler to maintain and is what I have a prototype in mind for, but if you'd rather have the C version for raw speed I'm happy to go that way instead.

@taldcroft @mhvk — since you looked at the `Time` one, does this seem worth doing, and any preference on the approach?

Start in astropy/coordinates/angles/core.py at Angle.to_string, especially the np.vectorize(do_format) call and _decimal_to_sexagesimal_string. Review the existing options—sep, pad, precision, fields, alwayssign, latex, and unicode—and compare the benchmark in the issue on large arrays. Done means preserving current formatting behavior while avoiding per-element Python string assembly for Angle and the related coordinate string methods.

Written by the indexing model from the issue text.