More precise mean_var for non-Dask dense arrays?

In #127, I changed how mean_var works just for Dask:

https://github.com/scverse/fast-array-utils/blob/940727fe307d7c6e0a3d0a2b8907860a285a11c0/src/fast_array_utils/stats/_mean_var.py#L40

We already have a custom numba implementation for sparse-in-core which operates on float64s, but doesn’t convert everything at once.

We could make the change also affect dense-in-core. Should we

• leave things as they are and keep power for float32 dense arrays as float32? It takes less space and is faster, but less precise.
• do exponentiation in 64 bit everywhere? this would allocate an exponentiated array of the input size in 64 bit, but I guess if people’s machines can handle a temporary 2n (prod(shape)×32×2), they can handle a temporary 3n (prod(shape)×32 + prod(shape)×64)?
• implement a custom numba implementation for dense-in-core as well?

Start by reading src/fast_array_utils/stats/_mean_var.py around line 40 and the change described in #127. Compare the dense-in-core, Dask, and sparse-in-core paths, focusing on float32 versus float64 exponentiation and temporary allocations. Done means selecting and documenting one precision strategy for dense arrays, with its memory and performance trade-offs clear.