Faster implementation of var across dimensions

Previously, we computed `var` across dimensions using two temporaries of the same size as the input array. This implementation (using Cartesian) avoids those temporaries. (One temporary is still necessary to hold the mean unless it is provided.)

For the code, after warmup:

```julia
a = randn(10000, 10000);
@time var(a, 2);
@time Base.varzm(a, 2);
```

Before:

```julia
elapsed time: 0.514875111 seconds (1600241672 bytes allocated)
elapsed time: 0.257889083 seconds (800177940 bytes allocated)
```

After:

```julia
elapsed time: 0.172945735 seconds (160960 bytes allocated)
elapsed time: 0.111020754 seconds (160512 bytes allocated)
```

A caveat: The old implementation used broadcasting to remove the mean, so it could work even if the size of the mean did not match the size of the output of the reduction. I didn't implement that here because it seemed unlikely to be used, but I could if there is reason to.

A question: Now that we have a `mean` keyword argument for `var`/`std`, can we deprecate `stdm`/`varm`?

cc @lindahua, @timholy