New implementation of cov/cor with extended API

This PR is made based on #5249.

* Here is the extended API:
```julia
cov(x[, corrected=true, vardim=1, zeromean=false])
cov(x, y[, corrected=true, vardim=1, zeromean=false])
cor(x[, vardim=1, zeromean=false])
cor(x, y[, vardim=1, zeromean=false])

covm(x, xmean[, corrected=true, vardim=1])
covm(x, xmean, y, ymean[, corrected=true, vardim=1])
corm(x, xmean[, vardim=1])
corm(x, xmean, y, ymean[, vardim=1])
```
Here, the arguments ``x`` and ``y`` can be either vectors or matrices.

Below is the explanation of keyword arguments:
     - ``corrected``: scale the matrix by ``1 / (n - 1)`` if set to true (by default), or ``1 / n`` otherwise. Note that this option does not apply to ``cor`` and ``corm``, as it does not affect the resultant correlation values.
     - ``vardim``: consider variables in columns and observations in rows if set to 1 (by default), or the other way round, if set to 2.
     - ``zeromean``: ``x`` and ``y`` are considered to have zero means, if set to true. Otherwise, the means will be computed from ``x`` and ``y`` and subtracted therefrom, if set to false (default).

* The computation is internally delegated to ``covzm`` and ``corzm``, which take centered data. These functions further delegate to ``unscaled_covzm``, which is the core implementation. All these internal functions are not exported.

* The functions are thoroughly tested -- I expand the test cases to ensure correctness under all combinations of keyword arguments.