diff --git a/base/array.jl b/base/array.jl
index 7435314caac50..70f721e6b7be7 100644
--- a/base/array.jl
+++ b/base/array.jl
@@ -2013,7 +2013,9 @@ julia> findmax([1,7,7,NaN])
 (NaN, 4)
 ```
 """
-function findmax(a)
+findmax(a) = _findmax(a, :)
+
+function _findmax(a, ::Colon)
     if isempty(a)
         throw(ArgumentError("collection must be non-empty"))
     end
@@ -2054,7 +2056,9 @@ julia> findmin([7,1,1,NaN])
 (NaN, 4)
 ```
 """
-function findmin(a)
+findmin(a) = _findmin(a, :)
+
+function _findmin(a, ::Colon)
     if isempty(a)
         throw(ArgumentError("collection must be non-empty"))
     end
diff --git a/base/bitarray.jl b/base/bitarray.jl
index e8b2aaa3b0008..9e6d65b76f908 100644
--- a/base/bitarray.jl
+++ b/base/bitarray.jl
@@ -1639,8 +1639,8 @@ findall(::typeof(!iszero), B::BitArray) = findall(B)
 
 ## Reductions ##
 
-sum(A::BitArray, region) = reducedim(+, A, region)
-sum(B::BitArray) = count(B)
+_sum(A::BitArray, dims)    = reduce(+, A, dims=dims)
+_sum(B::BitArray, ::Colon) = count(B)
 
 function all(B::BitArray)
     isempty(B) && return true
diff --git a/base/deprecated.jl b/base/deprecated.jl
index aa577de658596..744881d3588eb 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -296,9 +296,9 @@ export conv, conv2, deconv, filt, filt!, xcorr
 
 # PR #21709
 @deprecate cov(x::AbstractVector, corrected::Bool) cov(x, corrected=corrected)
-@deprecate cov(x::AbstractMatrix, vardim::Int, corrected::Bool) cov(x, vardim, corrected=corrected)
+@deprecate cov(x::AbstractMatrix, vardim::Int, corrected::Bool) cov(x, dims=vardim, corrected=corrected)
 @deprecate cov(X::AbstractVector, Y::AbstractVector, corrected::Bool) cov(X, Y, corrected=corrected)
-@deprecate cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, vardim::Int, corrected::Bool) cov(X, Y, vardim, corrected=corrected)
+@deprecate cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, vardim::Int, corrected::Bool) cov(X, Y, dims=vardim, corrected=corrected)
 
 # PR #22325
 # TODO: when this replace is removed from deprecated.jl:
@@ -781,8 +781,8 @@ findprev(pred::Function, A, i::Integer) = invoke(findprev, Tuple{Function, Any,
 # issue #20899
 # TODO: delete JULIA_HOME deprecation in src/init.c
 
-@deprecate cumsum(A::AbstractArray)     cumsum(A, 1)
-@deprecate cumprod(A::AbstractArray)    cumprod(A, 1)
+# cumsum and cumprod have deprecations in multidimensional.jl
+# when the message is removed, the `dims` keyword argument should become required.
 
 # issue #16307
 @deprecate finalizer(o, f::Function) finalizer(f, o)
@@ -1305,6 +1305,47 @@ export readandwrite
 @deprecate datatype_name(t::DataType) nameof(t) false
 @deprecate datatype_name(t::UnionAll) nameof(t) false
 
+# issue #25501
+@deprecate sum(a::AbstractArray, dims)        sum(a, dims=dims)
+@deprecate sum(f, a::AbstractArray, dims)     sum(f, a, dims=dims)
+@deprecate prod(a::AbstractArray, dims)       prod(a, dims=dims)
+@deprecate prod(f, a::AbstractArray, dims)    prod(f, a, dims=dims)
+@deprecate maximum(a::AbstractArray, dims)    maximum(a, dims=dims)
+@deprecate maximum(f, a::AbstractArray, dims) maximum(f, a, dims=dims)
+@deprecate minimum(a::AbstractArray, dims)    minimum(a, dims=dims)
+@deprecate minimum(f, a::AbstractArray, dims) minimum(f, a, dims=dims)
+@deprecate all(a::AbstractArray, dims)        all(a, dims=dims)
+@deprecate all(f, a::AbstractArray, dims)     all(f, a, dims=dims)
+@deprecate any(a::AbstractArray, dims)        any(a, dims=dims)
+@deprecate any(f, a::AbstractArray, dims)     any(f, a, dims=dims)
+@deprecate findmax(A::AbstractArray, dims)    findmax(A, dims=dims)
+@deprecate findmin(A::AbstractArray, dims)    findmin(A, dims=dims)
+
+@deprecate mean(A::AbstractArray, dims)                              mean(A, dims=dims)
+@deprecate varm(A::AbstractArray, m::AbstractArray, dims; kwargs...) varm(A, m; kwargs..., dims=dims)
+@deprecate var(A::AbstractArray, dims; kwargs...)                    var(A; kwargs..., dims=dims)
+@deprecate std(A::AbstractArray, dims; kwargs...)                    std(A; kwargs..., dims=dims)
+@deprecate cov(X::AbstractMatrix, dim::Int; kwargs...)               cov(X; kwargs..., dims=dim)
+@deprecate cov(x::AbstractVecOrMat, y::AbstractVecOrMat, dim::Int; kwargs...) cov(x, y; kwargs..., dims=dim)
+@deprecate cor(X::AbstractMatrix, dim::Int)                          cor(X, dims=dim)
+@deprecate cor(x::AbstractVecOrMat, y::AbstractVecOrMat, dim::Int)   cor(x, y, dims=dim)
+@deprecate median(A::AbstractArray, dims; kwargs...)                 median(A; kwargs..., dims=dims)
+
+@deprecate mapreducedim(f, op, A::AbstractArray, dims)     mapreduce(f, op, A, dims=dims)
+@deprecate mapreducedim(f, op, A::AbstractArray, dims, v0) mapreduce(f, op, v0, A, dims=dims)
+@deprecate reducedim(op, A::AbstractArray, dims)           reduce(op, A, dims=dims)
+@deprecate reducedim(op, A::AbstractArray, dims, v0)       reduce(op, v0, A, dims=dims)
+
+@deprecate sort(A::AbstractArray, dim::Integer; kwargs...) sort(A; kwargs..., dims=dim)
+
+@deprecate accumulate(op, A, dim::Integer)               accumulate(op, A, dims=dim)
+@deprecate accumulate!(op, B, A, dim::Integer)           accumulate!(op, B, A, dims=dim)
+@deprecate cumsum(A::AbstractArray, dim::Integer)        cumsum(A, dims=dim)
+@deprecate cumsum!(B, A, dim::Integer)                   cumsum!(B, A, dims=dim)
+@deprecate cumsum!(out, A::AbstractVector, dim::Integer) cumsum!(out, A, dims=dim)
+@deprecate cumprod(A::AbstractArray, dim::Integer)       cumprod(A, dims=dim)
+@deprecate cumprod!(B, A, dim::Integer)                  cumprod!(B, A, dims=dim)
+
 # PR #25196
 @deprecate_binding ObjectIdDict IdDict{Any,Any}
 
diff --git a/base/exports.jl b/base/exports.jl
index 4f9527ef2955f..e8e0bd1cdb40e 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -417,7 +417,6 @@ export
     prod,
     promote_shape,
     range,
-    reducedim,
     reshape,
     reverse!,
     reverse,
@@ -526,7 +525,6 @@ export
     mapfoldl,
     mapfoldr,
     mapreduce,
-    mapreducedim,
     merge!,
     merge,
     pairs,
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 5517f449c0e5c..5147f7f9f058d 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -737,9 +737,9 @@ function accumulate_pairwise(op, v::AbstractVector{T}) where T
     return accumulate_pairwise!(op, out, v)
 end
 
-function cumsum!(out, v::AbstractVector, dim::Integer)
+function cumsum!(out, v::AbstractVector; dims::Integer=1)
     # we dispatch on the possibility of numerical stability issues
-    _cumsum!(out, v, dim, ArithmeticStyle(eltype(out)))
+    _cumsum!(out, v, dims, ArithmeticStyle(eltype(out)))
 end
 
 function _cumsum!(out, v, dim, ::ArithmeticRounds)
@@ -753,9 +753,9 @@ function _cumsum!(out, v, dim, ::ArithmeticStyle)
 end
 
 """
-    cumsum(A, dim::Integer)
+    cumsum(A; dims::Integer)
 
-Cumulative sum along the dimension `dim`. See also [`cumsum!`](@ref)
+Cumulative sum along the dimension `dims`. See also [`cumsum!`](@ref)
 to use a preallocated output array, both for performance and to control the precision of the
 output (e.g. to avoid overflow).
 
@@ -766,20 +766,24 @@ julia> a = [1 2 3; 4 5 6]
  1  2  3
  4  5  6
 
-julia> cumsum(a,1)
+julia> cumsum(a, dims=1)
 2×3 Array{Int64,2}:
  1  2  3
  5  7  9
 
-julia> cumsum(a,2)
+julia> cumsum(a, dims=2)
 2×3 Array{Int64,2}:
  1  3   6
  4  9  15
 ```
 """
-function cumsum(A::AbstractArray{T}, dim::Integer) where T
+function cumsum(A::AbstractArray{T}; dims::Union{Nothing,Integer}=nothing) where T
+    if dims === nothing
+        depwarn("`cumsum(A::AbstractArray)` is deprecated, use `cumsum(A, dims=1)` instead.", :cumsum)
+        dims = 1
+    end
     out = similar(A, rcum_promote_type(+, T))
-    cumsum!(out, A, dim)
+    cumsum!(out, A, dims=dims)
 end
 
 """
@@ -804,24 +808,17 @@ julia> cumsum([fill(1, 2) for i in 1:3])
  [3, 3]
 ```
 """
-cumsum(x::AbstractVector) = cumsum(x, 1)
-
-"""
-    cumsum!(B, A, dim::Integer)
+cumsum(x::AbstractVector) = cumsum(x, dims=1)
 
-Cumulative sum of `A` along the dimension `dim`, storing the result in `B`. See also [`cumsum`](@ref).
 """
-cumsum!(B, A, dim::Integer) = accumulate!(+, B, A, dim)
+    cumsum!(B, A; dims::Integer)
 
+Cumulative sum of `A` along the dimension `dims`, storing the result in `B`. See also [`cumsum`](@ref).
 """
-    cumsum!(y::AbstractVector, x::AbstractVector)
+cumsum!(B, A; dims::Integer) = accumulate!(+, B, A, dims=dims)
 
-Cumulative sum of a vector `x`, storing the result in `y`. See also [`cumsum`](@ref).
 """
-cumsum!(y::AbstractVector, x::AbstractVector) = cumsum!(y, x, 1)
-
-"""
-    cumprod(A, dim::Integer)
+    cumprod(A; dims::Integer)
 
 Cumulative product along the dimension `dim`. See also
 [`cumprod!`](@ref) to use a preallocated output array, both for performance and
@@ -834,18 +831,24 @@ julia> a = [1 2 3; 4 5 6]
  1  2  3
  4  5  6
 
-julia> cumprod(a,1)
+julia> cumprod(a, dims=1)
 2×3 Array{Int64,2}:
  1   2   3
  4  10  18
 
-julia> cumprod(a,2)
+julia> cumprod(a, dims=2)
 2×3 Array{Int64,2}:
  1   2    6
  4  20  120
 ```
 """
-cumprod(A::AbstractArray, dim::Integer) = accumulate(*, A, dim)
+function cumprod(A::AbstractArray; dims::Union{Nothing,Integer}=nothing)
+    if dims === nothing
+        depwarn("`cumprod(A::AbstractArray)` is deprecated, use `cumprod(A, dims=1)` instead.", :cumprod)
+        dims = 1
+    end
+    return accumulate(*, A, dims=dims)
+end
 
 """
     cumprod(x::AbstractVector)
@@ -869,15 +872,15 @@ julia> cumprod([fill(1//3, 2, 2) for i in 1:3])
  Rational{Int64}[4//27 4//27; 4//27 4//27]
 ```
 """
-cumprod(x::AbstractVector) = cumprod(x, 1)
+cumprod(x::AbstractVector) = cumprod(x, dims=1)
 
 """
-    cumprod!(B, A, dim::Integer)
+    cumprod!(B, A; dims::Integer)
 
-Cumulative product of `A` along the dimension `dim`, storing the result in `B`.
+Cumulative product of `A` along the dimension `dims`, storing the result in `B`.
 See also [`cumprod`](@ref).
 """
-cumprod!(B, A, dim::Integer) = accumulate!(*, B, A, dim)
+cumprod!(B, A; dims::Integer) = accumulate!(*, B, A, dims=dims)
 
 """
     cumprod!(y::AbstractVector, x::AbstractVector)
@@ -885,12 +888,12 @@ cumprod!(B, A, dim::Integer) = accumulate!(*, B, A, dim)
 Cumulative product of a vector `x`, storing the result in `y`.
 See also [`cumprod`](@ref).
 """
-cumprod!(y::AbstractVector, x::AbstractVector) = cumprod!(y, x, 1)
+cumprod!(y::AbstractVector, x::AbstractVector) = cumprod!(y, x, dims=1)
 
 """
-    accumulate(op, A, dim::Integer)
+    accumulate(op, A; dims::Integer)
 
-Cumulative operation `op` along the dimension `dim`. See also
+Cumulative operation `op` along the dimension `dims`. See also
 [`accumulate!`](@ref) to use a preallocated output array, both for performance and
 to control the precision of the output (e.g. to avoid overflow). For common operations
 there are specialized variants of `accumulate`, see:
@@ -898,22 +901,22 @@ there are specialized variants of `accumulate`, see:
 
 # Examples
 ```jldoctest
-julia> accumulate(+, fill(1, 3, 3), 1)
+julia> accumulate(+, fill(1, 3, 3), dims=1)
 3×3 Array{Int64,2}:
  1  1  1
  2  2  2
  3  3  3
 
-julia> accumulate(+, fill(1, 3, 3), 2)
+julia> accumulate(+, fill(1, 3, 3), dims=2)
 3×3 Array{Int64,2}:
  1  2  3
  1  2  3
  1  2  3
 ```
 """
-function accumulate(op, A, dim::Integer)
+function accumulate(op, A; dims::Integer)
     out = similar(A, rcum_promote_type(op, eltype(A)))
-    accumulate!(op, out, A, dim)
+    accumulate!(op, out, A, dims=dims)
 end
 
 """
@@ -940,12 +943,12 @@ julia> accumulate(*, [1,2,3])
  6
 ```
 """
-accumulate(op, x::AbstractVector) = accumulate(op, x, 1)
+accumulate(op, x::AbstractVector) = accumulate(op, x, dims=1)
 
 """
-    accumulate!(op, B, A, dim::Integer)
+    accumulate!(op, B, A; dims::Integer)
 
-Cumulative operation `op` on `A` along the dimension `dim`, storing the result in `B`.
+Cumulative operation `op` on `A` along the dimension `dims`, storing the result in `B`.
 See also [`accumulate`](@ref).
 
 # Examples
@@ -954,14 +957,14 @@ julia> A = [1 2; 3 4];
 
 julia> B = [0 0; 0 0];
 
-julia> accumulate!(-, B, A, 1);
+julia> accumulate!(-, B, A, dims=1);
 
 julia> B
 2×2 Array{Int64,2}:
   1   2
  -2  -2
 
-julia> accumulate!(-, B, A, 2);
+julia> accumulate!(-, B, A, dims=2);
 
 julia> B
 2×2 Array{Int64,2}:
@@ -969,7 +972,8 @@ julia> B
  3  -1
 ```
 """
-function accumulate!(op, B, A, dim::Integer)
+function accumulate!(op, B, A; dims::Integer)
+    dim = dims
     dim > 0 || throw(ArgumentError("dim must be a positive integer"))
     inds_t = axes(A)
     axes(B) == inds_t || throw(DimensionMismatch("shape of B must match A"))
diff --git a/base/reduce.jl b/base/reduce.jl
index b60f82d817ae8..12b231bda62e5 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -353,11 +353,10 @@ function _mapreduce(f, op, ::IndexLinear, A::AbstractArray{T}) where T
     end
 end
 
-_mapreduce(f, op, ::IndexCartesian, A::AbstractArray) = mapfoldl(f, op, A)
-
-mapreduce(f, op, A::AbstractArray) = _mapreduce(f, op, IndexStyle(A), A)
 mapreduce(f, op, a::Number) = mapreduce_first(f, op, a)
 
+_mapreduce(f, op, ::IndexCartesian, A::AbstractArray) = mapfoldl(f, op, A)
+
 """
     reduce(op, v0, itr)
 
@@ -639,7 +638,9 @@ julia> any(i -> (println(i); i > 3), 1:10)
 true
 ```
 """
-function any(f, itr)
+any(f, itr) = _any(f, itr, :)
+
+function _any(f, itr, ::Colon)
     anymissing = false
     for x in itr
         v = f(x)
@@ -673,7 +674,9 @@ julia> all(i -> (println(i); i < 3), 1:10)
 false
 ```
 """
-function all(f, itr)
+all(f, itr) = _all(f, itr, :)
+
+function _all(f, itr, ::Colon)
     anymissing = false
     for x in itr
         v = f(x)
diff --git a/base/reducedim.jl b/base/reducedim.jl
index dda3057a8764f..df3e249cee54a 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -82,7 +82,7 @@ reducedim_initarray(A::AbstractArray, region, v0::T) where {T} = reducedim_inita
 function reducedim_initarray0(A::AbstractArray{T}, region, f, ops) where T
     ri = reduced_indices0(A, region)
     if isempty(A)
-        if prod(map(length, reduced_indices(A, region))) != 0
+        if prod(length, reduced_indices(A, region)) != 0
             reducedim_initarray0_empty(A, region, f, ops) # ops over empty slice of A
         else
             R = f == identity ? T : Core.Compiler.return_type(f, (T,))
@@ -155,10 +155,10 @@ function check_reducedims(R, A)
     # Check whether R has compatible dimensions w.r.t. A for reduction
     #
     # It returns an integer value (useful for choosing implementation)
-    # - If it reduces only along leading dimensions, e.g. sum(A, 1) or sum(A, (1, 2)),
+    # - If it reduces only along leading dimensions, e.g. sum(A, dims=1) or sum(A, dims=(1,2)),
     #   it returns the length of the leading slice. For the two examples above,
     #   it will be size(A, 1) or size(A, 1) * size(A, 2).
-    # - Otherwise, e.g. sum(A, 2) or sum(A, (1, 3)), it returns 0.
+    # - Otherwise, e.g. sum(A, dims=2) or sum(A, dims=(1,3)), it returns 0.
     #
     ndims(R) <= ndims(A) || throw(DimensionMismatch("cannot reduce $(ndims(A))-dimensional array to $(ndims(R)) dimensions"))
     lsiz = 1
@@ -244,9 +244,9 @@ reducedim!(op, R::AbstractArray{RT}, A::AbstractArray) where {RT} =
     mapreducedim!(identity, op, R, A)
 
 """
-    mapreducedim(f, op, A, region[, v0])
+    mapreduce(f, op[, v0], A::AbstractArray; dims)
 
-Evaluates to the same as `reducedim(op, map(f, A), region, f(v0))`, but is generally
+Evaluates to the same as `reduce(op, f(v0), map(f, A); dims)`, but is generally
 faster because the intermediate array is avoided.
 
 # Examples
@@ -258,24 +258,33 @@ julia> a = reshape(Vector(1:16), (4,4))
  3  7  11  15
  4  8  12  16
 
-julia> mapreducedim(isodd, *, a, 1)
+julia> mapreduce(isodd, *, a, dims=1)
 1×4 Array{Bool,2}:
  false  false  false  false
 
-julia> mapreducedim(isodd, |, a, 1, true)
+julia> mapreduce(isodd, |, true, a, dims=1)
 1×4 Array{Bool,2}:
  true  true  true  true
 ```
 """
-mapreducedim(f, op, A::AbstractArray, region, v0) =
-    mapreducedim!(f, op, reducedim_initarray(A, region, v0), A)
-mapreducedim(f, op, A::AbstractArray, region) =
-    mapreducedim!(f, op, reducedim_init(f, op, A, region), A)
+mapreduce(f, op, v0, A::AbstractArray; dims=:) = _mapreduce_dim(f, op, v0, A, dims)
+
+mapreduce(f, op, A::AbstractArray; dims=:) = _mapreduce_dim(f, op, A, dims)
+
+_mapreduce_dim(f, op, v0, A::AbstractArray, ::Colon) = mapfoldl(f, op, v0, A)
+
+_mapreduce_dim(f, op, A::AbstractArray, ::Colon) = _mapreduce(f, op, IndexStyle(A), A)
+
+_mapreduce_dim(f, op, v0, A::AbstractArray, dims) =
+    mapreducedim!(f, op, reducedim_initarray(A, dims, v0), A)
+
+_mapreduce_dim(f, op, A::AbstractArray, dims) =
+    mapreducedim!(f, op, reducedim_init(f, op, A, dims), A)
 
 """
-    reducedim(f, A, region[, v0])
+    reduce(f[, v0], A; dims)
 
-Reduce 2-argument function `f` along dimensions of `A`. `region` is a vector specifying the
+Reduce 2-argument function `f` along dimensions of `A`. `dims` is a vector specifying the
 dimensions to reduce, and `v0` is the initial value to use in the reductions. For `+`, `*`,
 `max` and `min` the `v0` argument is optional.
 
@@ -292,24 +301,31 @@ julia> a = reshape(Vector(1:16), (4,4))
  3  7  11  15
  4  8  12  16
 
-julia> reducedim(max, a, 2)
+julia> reduce(max, a, dims=2)
 4×1 Array{Int64,2}:
  13
  14
  15
  16
 
-julia> reducedim(max, a, 1)
+julia> reduce(max, a, dims=1)
 1×4 Array{Int64,2}:
  4  8  12  16
 ```
 """
-reducedim(op, A::AbstractArray, region, v0) = mapreducedim(identity, op, A, region, v0)
-reducedim(op, A::AbstractArray, region) = mapreducedim(identity, op, A, region)
+reduce(op, v0, A::AbstractArray; dims=:) = _reduce_dim(op, v0, A, dims)
+
+_reduce_dim(op, v0, A, dims) = mapreduce(identity, op, v0, A, dims=dims)
+_reduce_dim(op, v0, A, ::Colon) = mapreduce(identity, op, v0, A)
+
+reduce(op, A::AbstractArray; dims=:) = _reduce_dim(op, A, dims)
+
+_reduce_dim(op, A, dims) = mapreduce(identity, op, A, dims=dims)
+_reduce_dim(op, A, ::Colon) = mapreduce(identity, op, A)
 
 ##### Specific reduction functions #####
 """
-    sum(A::AbstractArray, dims)
+    sum(A::AbstractArray; dims)
 
 Sum elements of an array over the given dimensions.
 
@@ -320,17 +336,17 @@ julia> A = [1 2; 3 4]
  1  2
  3  4
 
-julia> sum(A, 1)
+julia> sum(A, dims=1)
 1×2 Array{Int64,2}:
  4  6
 
-julia> sum(A, 2)
+julia> sum(A, dims=2)
 2×1 Array{Int64,2}:
  3
  7
 ```
 """
-sum(A::AbstractArray, dims)
+sum(A::AbstractArray; dims)
 
 """
     sum!(r, A)
@@ -357,7 +373,7 @@ julia> sum!([1 1], A)
 sum!(r, A)
 
 """
-    prod(A::AbstractArray, dims)
+    prod(A::AbstractArray; dims)
 
 Multiply elements of an array over the given dimensions.
 
@@ -368,17 +384,17 @@ julia> A = [1 2; 3 4]
  1  2
  3  4
 
-julia> prod(A, 1)
+julia> prod(A, dims=1)
 1×2 Array{Int64,2}:
  3  8
 
-julia> prod(A, 2)
+julia> prod(A, dims=2)
 2×1 Array{Int64,2}:
   2
  12
 ```
 """
-prod(A::AbstractArray, dims)
+prod(A::AbstractArray; dims)
 
 """
     prod!(r, A)
@@ -405,7 +421,7 @@ julia> prod!([1 1], A)
 prod!(r, A)
 
 """
-    maximum(A, dims)
+    maximum(A::AbstractArray; dims)
 
 Compute the maximum value of an array over the given dimensions. See also the
 [`max(a,b)`](@ref) function to take the maximum of two or more arguments,
@@ -417,17 +433,17 @@ julia> A = [1 2; 3 4]
  1  2
  3  4
 
-julia> maximum(A, 1)
+julia> maximum(A, dims=1)
 1×2 Array{Int64,2}:
  3  4
 
-julia> maximum(A, 2)
+julia> maximum(A, dims=2)
 2×1 Array{Int64,2}:
  2
  4
 ```
 """
-maximum(A, dims)
+maximum(A::AbstractArray; dims)
 
 """
     maximum!(r, A)
@@ -454,7 +470,7 @@ julia> maximum!([1 1], A)
 maximum!(r, A)
 
 """
-    minimum(A, dims)
+    minimum(A::AbstractArray; dims)
 
 Compute the minimum value of an array over the given dimensions. See also the
 [`min(a,b)`](@ref) function to take the minimum of two or more arguments,
@@ -467,17 +483,17 @@ julia> A = [1 2; 3 4]
  1  2
  3  4
 
-julia> minimum(A, 1)
+julia> minimum(A, dims=1)
 1×2 Array{Int64,2}:
  1  2
 
-julia> minimum(A, 2)
+julia> minimum(A, dims=2)
 2×1 Array{Int64,2}:
  1
  3
 ```
 """
-minimum(A, dims)
+minimum(A::AbstractArray; dims)
 
 """
     minimum!(r, A)
@@ -504,7 +520,7 @@ julia> minimum!([1 1], A)
 minimum!(r, A)
 
 """
-    all(A, dims)
+    all(A; dims)
 
 Test whether all values along the given dimensions of an array are `true`.
 
@@ -515,17 +531,17 @@ julia> A = [true false; true true]
  true  false
  true   true
 
-julia> all(A, 1)
+julia> all(A, dims=1)
 1×2 Array{Bool,2}:
  true  false
 
-julia> all(A, 2)
+julia> all(A, dims=2)
 2×1 Array{Bool,2}:
  false
   true
 ```
 """
-all(A::AbstractArray, dims)
+all(A::AbstractArray; dims)
 
 """
     all!(r, A)
@@ -552,7 +568,7 @@ julia> all!([1 1], A)
 all!(r, A)
 
 """
-    any(A, dims)
+    any(A; dims)
 
 Test whether any values along the given dimensions of an array are `true`.
 
@@ -563,17 +579,17 @@ julia> A = [true false; true false]
  true  false
  true  false
 
-julia> any(A, 1)
+julia> any(A, dims=1)
 1×2 Array{Bool,2}:
  true  false
 
-julia> any(A, 2)
+julia> any(A, dims=2)
 2×1 Array{Bool,2}:
  true
  true
 ```
 """
-any(::AbstractArray,dims)
+any(::AbstractArray; dims)
 
 """
     any!(r, A)
@@ -600,18 +616,38 @@ julia> any!([1 1], A)
 """
 any!(r, A)
 
+for (fname, _fname, op) in [(:sum,     :_sum,     :add_sum), (:prod,    :_prod,    :mul_prod),
+                            (:maximum, :_maximum, :max),     (:minimum, :_minimum, :min)]
+    @eval begin
+        # User-facing methods with keyword arguments
+        @inline ($fname)(a::AbstractArray; dims=:) = ($_fname)(a, dims)
+        @inline ($fname)(f::Callable, a::AbstractArray; dims=:) = ($_fname)(f, a, dims)
+
+        # Underlying implementations using dispatch
+        ($_fname)(a, ::Colon) = ($_fname)(identity, a, :)
+        ($_fname)(f, a, ::Colon) = mapreduce(f, $op, a)
+    end
+end
+
+any(a::AbstractArray; dims=:)              = _any(a, dims)
+any(f::Function, a::AbstractArray; dims=:) = _any(f, a, dims)
+_any(a, ::Colon)                           = _any(identity, a, :)
+all(a::AbstractArray; dims=:)              = _all(a, dims)
+all(f::Function, a::AbstractArray; dims=:) = _all(f, a, dims)
+_all(a, ::Colon)                           = _all(identity, a, :)
+
 for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod),
                     (:maximum, :max), (:minimum, :min),
-                    (:all, :&), (:any, :|)]
+                    (:all, :&),       (:any, :|)]
     fname! = Symbol(fname, '!')
+    _fname = Symbol('_', fname)
     @eval begin
         $(fname!)(f::Function, r::AbstractArray, A::AbstractArray; init::Bool=true) =
             mapreducedim!(f, $(op), initarray!(r, $(op), init, A), A)
         $(fname!)(r::AbstractArray, A::AbstractArray; init::Bool=true) = $(fname!)(identity, r, A; init=init)
 
-        $(fname)(f::Function, A::AbstractArray, region) =
-            mapreducedim(f, $(op), A, region)
-        $(fname)(A::AbstractArray, region) = $(fname)(identity, A, region)
+        $(_fname)(A, dims)    = $(_fname)(identity, A, dims)
+        $(_fname)(f, A, dims) = mapreduce(f, $(op), A, dims=dims)
     end
 end
 
@@ -679,9 +715,9 @@ function findmin!(rval::AbstractArray, rind::AbstractArray, A::AbstractArray;
 end
 
 """
-    findmin(A, region) -> (minval, index)
+    findmin(A; dims) -> (minval, index)
 
-For an array input, returns the value and index of the minimum over the given region.
+For an array input, returns the value and index of the minimum over the given dimensions.
 `NaN` is treated as less than all other values.
 
 # Examples
@@ -691,14 +727,16 @@ julia> A = [1.0 2; 3 4]
  1.0  2.0
  3.0  4.0
 
-julia> findmin(A, 1)
+julia> findmin(A, dims=1)
 ([1.0 2.0], CartesianIndex{2}[CartesianIndex(1, 1) CartesianIndex(1, 2)])
 
-julia> findmin(A, 2)
+julia> findmin(A, dims=2)
 ([1.0; 3.0], CartesianIndex{2}[CartesianIndex(1, 1); CartesianIndex(2, 1)])
 ```
 """
-function findmin(A::AbstractArray{T}, region) where T
+findmin(A::AbstractArray; dims=:) = _findmin(A, dims)
+
+function _findmin(A, region)
     ri = reduced_indices0(A, region)
     if isempty(A)
         if prod(map(length, reduced_indices(A, region))) != 0
@@ -726,9 +764,9 @@ function findmax!(rval::AbstractArray, rind::AbstractArray, A::AbstractArray;
 end
 
 """
-    findmax(A, region) -> (maxval, index)
+    findmax(A; dims) -> (maxval, index)
 
-For an array input, returns the value and index of the maximum over the given region.
+For an array input, returns the value and index of the maximum over the given dimensions.
 `NaN` is treated as greater than all other values.
 
 # Examples
@@ -738,14 +776,16 @@ julia> A = [1.0 2; 3 4]
  1.0  2.0
  3.0  4.0
 
-julia> findmax(A,1)
+julia> findmax(A, dims=1)
 ([3.0 4.0], CartesianIndex{2}[CartesianIndex(2, 1) CartesianIndex(2, 2)])
 
-julia> findmax(A,2)
+julia> findmax(A, dims=2)
 ([2.0; 4.0], CartesianIndex{2}[CartesianIndex(1, 2); CartesianIndex(2, 2)])
 ```
 """
-function findmax(A::AbstractArray{T}, region) where T
+findmax(A::AbstractArray; dims=:) = _findmax(A, dims)
+
+function _findmax(A, region)
     ri = reduced_indices0(A, region)
     if isempty(A)
         if prod(map(length, reduced_indices(A, region))) != 0
diff --git a/base/sort.jl b/base/sort.jl
index f56081b982179..5713bd729b89f 100644
--- a/base/sort.jl
+++ b/base/sort.jl
@@ -858,13 +858,15 @@ julia> sort(A, 2)
  1  2
 ```
 """
-function sort(A::AbstractArray, dim::Integer;
+function sort(A::AbstractArray;
+              dims::Integer,
               alg::Algorithm=DEFAULT_UNSTABLE,
               lt=isless,
               by=identity,
               rev::Union{Bool,Nothing}=nothing,
               order::Ordering=Forward,
               initialized::Union{Bool,Nothing}=nothing)
+    dim = dims
     if initialized !== nothing
         Base.depwarn("`initialized` keyword argument is deprecated", :sort)
     end
diff --git a/base/statistics.jl b/base/statistics.jl
index 9a008191b9905..e358e197d6cad 100644
--- a/base/statistics.jl
+++ b/base/statistics.jl
@@ -33,7 +33,6 @@ function mean(f::Callable, iterable)
 end
 mean(iterable) = mean(identity, iterable)
 mean(f::Callable, A::AbstractArray) = sum(f, A) / _length(A)
-mean(A::AbstractArray) = sum(A) / _length(A)
 
 """
     mean!(r, v)
@@ -65,16 +64,18 @@ function mean!(R::AbstractArray, A::AbstractArray)
 end
 
 """
-    mean(v[, region])
+    mean(v; dims)
 
-Compute the mean of whole array `v`, or optionally along the dimensions in `region`.
+Compute the mean of whole array `v`, or optionally along the given dimensions.
 
 !!! note
     Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
     to represent missing values, and the [`skipmissing`](@ref) function to omit them.
 """
-mean(A::AbstractArray{T}, region) where {T} =
-    mean!(reducedim_init(t -> t/2, +, A, region), A)
+mean(A::AbstractArray; dims=:) = _mean(A, dims)
+
+_mean(A::AbstractArray{T}, region) where {T} = mean!(reducedim_init(t -> t/2, +, A, region), A)
+_mean(A::AbstractArray, ::Colon) = sum(A) / _length(A)
 
 ##### variances #####
 
@@ -82,7 +83,9 @@ mean(A::AbstractArray{T}, region) where {T} =
 realXcY(x::Real, y::Real) = x*y
 realXcY(x::Complex, y::Complex) = real(x)*real(y) + imag(x)*imag(y)
 
-function var(iterable; corrected::Bool=true, mean=nothing)
+var(iterable; corrected::Bool=true, mean=nothing) = _var(iterable, corrected, mean)
+
+function _var(iterable, corrected::Bool, mean)
     state = start(iterable)
     if done(iterable, state)
         throw(ArgumentError("variance of empty collection undefined: $(repr(iterable))"))
@@ -165,12 +168,6 @@ function centralize_sumabs2!(R::AbstractArray{S}, A::AbstractArray, means::Abstr
     return R
 end
 
-function varm(A::AbstractArray{T}, m; corrected::Bool=true) where T
-    n = _length(A)
-    n == 0 && return typeof((abs2(zero(T)) + abs2(zero(T)))/2)(NaN)
-    return centralize_sumabs2(A, m) / (n - Int(corrected))
-end
-
 function varm!(R::AbstractArray{S}, A::AbstractArray, m::AbstractArray; corrected::Bool=true) where S
     if isempty(A)
         fill!(R, convert(S, NaN))
@@ -183,10 +180,10 @@ function varm!(R::AbstractArray{S}, A::AbstractArray, m::AbstractArray; correcte
 end
 
 """
-    varm(v, m[, region]; corrected::Bool=true)
+    varm(v, m; dims, corrected::Bool=true)
 
 Compute the sample variance of a collection `v` with known mean(s) `m`,
-optionally over `region`. `m` may contain means for each dimension of
+optionally over the given dimensions. `m` may contain means for each dimension of
 `v`. If `corrected` is `true`, then the sum is scaled with `n-1`,
 whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`.
 
@@ -194,18 +191,25 @@ whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x
     Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
     to represent missing values, and the [`skipmissing`](@ref) function to omit them.
 """
-varm(A::AbstractArray{T}, m::AbstractArray, region; corrected::Bool=true) where {T} =
+varm(A::AbstractArray, m::AbstractArray; corrected::Bool=true, dims=:) = _varm(A, m, corrected, dims)
+
+_varm(A::AbstractArray{T}, m, corrected::Bool, region) where {T} =
     varm!(reducedim_init(t -> abs2(t)/2, +, A, region), A, m; corrected=corrected)
 
+varm(A::AbstractArray, m; corrected::Bool=true) = _varm(A, m, corrected, :)
+
+function _varm(A::AbstractArray{T}, m, corrected::Bool, ::Colon) where T
+    n = _length(A)
+    n == 0 && return typeof((abs2(zero(T)) + abs2(zero(T)))/2)(NaN)
+    return centralize_sumabs2(A, m) / (n - Int(corrected))
+end
 
-var(A::AbstractArray{T}; corrected::Bool=true, mean=nothing) where {T} =
-    real(varm(A, coalesce(mean, Base.mean(A)); corrected=corrected))
 
 """
-    var(v[, region]; corrected::Bool=true, mean=nothing)
+    var(v; dims, corrected::Bool=true, mean=nothing)
 
-Compute the sample variance of a vector or array `v`, optionally along dimensions in
-`region`. The algorithm will return an estimator of the generative distribution's variance
+Compute the sample variance of a vector or array `v`, optionally along the given dimensions.
+The algorithm will return an estimator of the generative distribution's variance
 under the assumption that each entry of `v` is an IID drawn from that generative
 distribution. This computation is equivalent to calculating `sum(abs2, v - mean(v)) /
 (length(v) - 1)`. If `corrected` is `true`, then the sum is scaled with `n-1`,
@@ -216,15 +220,22 @@ The mean `mean` over the region may be provided.
     Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
     to represent missing values, and the [`skipmissing`](@ref) function to omit them.
 """
-var(A::AbstractArray, region; corrected::Bool=true, mean=nothing) =
-    varm(A, coalesce(mean, Base.mean(A, region)), region; corrected=corrected)
+var(A::AbstractArray; corrected::Bool=true, mean=nothing, dims=:) = _var(A, corrected, mean, dims)
+
+_var(A::AbstractArray, corrected::Bool, mean, dims) =
+    varm(A, coalesce(mean, Base.mean(A, dims=dims)); corrected=corrected, dims=dims)
+
+_var(A::AbstractArray, corrected::Bool, mean, ::Colon) =
+    real(varm(A, coalesce(mean, Base.mean(A)); corrected=corrected))
 
-varm(iterable, m; corrected::Bool=true) =
-    var(iterable, corrected=corrected, mean=m)
+varm(iterable, m; corrected::Bool=true) = _var(iterable, corrected, m)
 
 ## variances over ranges
 
-function varm(v::AbstractRange, m)
+varm(v::AbstractRange, m::AbstractArray) = range_varm(v, m)
+varm(v::AbstractRange, m) = range_varm(v, m)
+
+function range_varm(v::AbstractRange, m)
     f  = first(v) - m
     s  = step(v)
     l  = length(v)
@@ -258,14 +269,11 @@ end
 stdm(A::AbstractArray, m; corrected::Bool=true) =
     sqrt.(varm(A, m; corrected=corrected))
 
-std(A::AbstractArray; corrected::Bool=true, mean=nothing) =
-    sqrt.(var(A; corrected=corrected, mean=mean))
-
 """
-    std(v[, region]; corrected::Bool=true, mean=nothing)
+    std(v; corrected::Bool=true, mean=nothing, dims)
 
-Compute the sample standard deviation of a vector or array `v`, optionally along dimensions
-in `region`. The algorithm returns an estimator of the generative distribution's standard
+Compute the sample standard deviation of a vector or array `v`, optionally along the given
+dimensions. The algorithm returns an estimator of the generative distribution's standard
 deviation under the assumption that each entry of `v` is an IID drawn from that generative
 distribution. This computation is equivalent to calculating `sqrt(sum((v - mean(v)).^2) /
 (length(v) - 1))`. A pre-computed `mean` may be provided. If `corrected` is `true`,
@@ -276,11 +284,19 @@ then the sum is scaled with `n-1`, whereas the sum is scaled with `n` if `correc
     Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
     to represent missing values, and the [`skipmissing`](@ref) function to omit them.
 """
-std(A::AbstractArray, region; corrected::Bool=true, mean=nothing) =
-    sqrt.(var(A, region; corrected=corrected, mean=mean))
+std(A::AbstractArray; corrected::Bool=true, mean=nothing, dims=:) = _std(A, corrected, mean, dims)
 
-std(A::AbstractArray{<:AbstractFloat}, region; corrected::Bool=true, mean=nothing) =
-    sqrt!(var(A, region; corrected=corrected, mean=mean))
+_std(A::AbstractArray, corrected::Bool, mean, dims) =
+    sqrt.(var(A; corrected=corrected, mean=mean, dims=dims))
+
+_std(A::AbstractArray, corrected::Bool, mean, ::Colon) =
+    sqrt.(var(A; corrected=corrected, mean=mean))
+
+_std(A::AbstractArray{<:AbstractFloat}, corrected::Bool, mean, dims) =
+    sqrt!(var(A; corrected=corrected, mean=mean, dims=dims))
+
+_std(A::AbstractArray{<:AbstractFloat}, corrected::Bool, mean, ::Colon) =
+    sqrt!(var(A; corrected=corrected, mean=mean))
 
 std(iterable; corrected::Bool=true, mean=nothing) =
     sqrt(var(iterable, corrected=corrected, mean=mean))
@@ -318,7 +334,7 @@ function _getnobs(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int)
 end
 
 _vmean(x::AbstractVector, vardim::Int) = mean(x)
-_vmean(x::AbstractMatrix, vardim::Int) = mean(x, vardim)
+_vmean(x::AbstractMatrix, vardim::Int) = mean(x, dims=vardim)
 
 # core functions
 
@@ -378,14 +394,14 @@ is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is `fals
 cov(x::AbstractVector; corrected::Bool=true) = covm(x, Base.mean(x); corrected=corrected)
 
 """
-    cov(X::AbstractMatrix[, vardim::Int=1]; corrected::Bool=true)
+    cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true)
 
-Compute the covariance matrix of the matrix `X` along the dimension `vardim`. If `corrected`
+Compute the covariance matrix of the matrix `X` along the dimension `dims`. If `corrected`
 is `true` (the default) then the sum is scaled with `n-1`, whereas the sum is scaled with `n`
-if `corrected` is `false` where `n = size(X, vardim)`.
+if `corrected` is `false` where `n = size(X, dims)`.
 """
-cov(X::AbstractMatrix, vardim::Int=1; corrected::Bool=true) =
-    covm(X, _vmean(X, vardim), vardim; corrected=corrected)
+cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true) =
+    covm(X, _vmean(X, dims), dims; corrected=corrected)
 
 """
     cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
@@ -399,14 +415,14 @@ cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true) =
     covm(x, Base.mean(x), y, Base.mean(y); corrected=corrected)
 
 """
-    cov(X::AbstractVecOrMat, Y::AbstractVecOrMat[, vardim::Int=1]; corrected::Bool=true)
+    cov(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims::Int=1, corrected::Bool=true)
 
 Compute the covariance between the vectors or matrices `X` and `Y` along the dimension
-`vardim`. If `corrected` is `true` (the default) then the sum is scaled with `n-1`, whereas
-the sum is scaled with `n` if `corrected` is `false` where `n = size(X, vardim) = size(Y, vardim)`.
+`dims`. If `corrected` is `true` (the default) then the sum is scaled with `n-1`, whereas
+the sum is scaled with `n` if `corrected` is `false` where `n = size(X, dims) = size(Y, dims)`.
 """
-cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, vardim::Int=1; corrected::Bool=true) =
-    covm(X, _vmean(X, vardim), Y, _vmean(Y, vardim), vardim; corrected=corrected)
+cov(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims::Int=1, corrected::Bool=true) =
+    covm(X, _vmean(X, dims), Y, _vmean(Y, dims), dims; corrected=corrected)
 
 ##### correlation #####
 
@@ -474,11 +490,11 @@ function corzm(x::AbstractMatrix, vardim::Int=1)
     return cov2cor!(c, collect(sqrt(c[i,i]) for i in 1:min(size(c)...)))
 end
 corzm(x::AbstractVector, y::AbstractMatrix, vardim::Int=1) =
-    cov2cor!(unscaled_covzm(x, y, vardim), sqrt(sum(abs2, x)), sqrt!(sum(abs2, y, vardim)))
+    cov2cor!(unscaled_covzm(x, y, vardim), sqrt(sum(abs2, x)), sqrt!(sum(abs2, y, dims=vardim)))
 corzm(x::AbstractMatrix, y::AbstractVector, vardim::Int=1) =
-    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, vardim)), sqrt(sum(abs2, y)))
+    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, dims=vardim)), sqrt(sum(abs2, y)))
 corzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int=1) =
-    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, vardim)), sqrt!(sum(abs2, y, vardim)))
+    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, dims=vardim)), sqrt!(sum(abs2, y, dims=vardim)))
 
 # corm
 
@@ -518,11 +534,11 @@ Return the number one.
 cor(x::AbstractVector) = one(real(eltype(x)))
 
 """
-    cor(X::AbstractMatrix[, vardim::Int=1])
+    cor(X::AbstractMatrix; dims::Int=1)
 
-Compute the Pearson correlation matrix of the matrix `X` along the dimension `vardim`.
+Compute the Pearson correlation matrix of the matrix `X` along the dimension `dims`.
 """
-cor(X::AbstractMatrix, vardim::Int=1) = corm(X, _vmean(X, vardim), vardim)
+cor(X::AbstractMatrix; dims::Int=1) = corm(X, _vmean(X, dims), dims)
 
 """
     cor(x::AbstractVector, y::AbstractVector)
@@ -532,12 +548,12 @@ Compute the Pearson correlation between the vectors `x` and `y`.
 cor(x::AbstractVector, y::AbstractVector) = corm(x, Base.mean(x), y, Base.mean(y))
 
 """
-    cor(X::AbstractVecOrMat, Y::AbstractVecOrMat[, vardim=1])
+    cor(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims=1)
 
-Compute the Pearson correlation between the vectors or matrices `X` and `Y` along the dimension `vardim`.
+Compute the Pearson correlation between the vectors or matrices `X` and `Y` along the dimension `dims`.
 """
-cor(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int=1) =
-    corm(x, _vmean(x, vardim), y, _vmean(y, vardim), vardim)
+cor(x::AbstractVecOrMat, y::AbstractVecOrMat; dims::Int=1) =
+    corm(x, _vmean(x, dims), y, _vmean(y, dims), dims)
 
 ##### median & quantiles #####
 
@@ -615,13 +631,12 @@ function median!(v::AbstractVector)
     end
 end
 median!(v::AbstractArray) = median!(vec(v))
-median(v::AbstractArray{T}) where {T} = median!(copyto!(Array{T,1}(uninitialized, _length(v)), v))
 
 """
-    median(v[, region])
+    median(v; dims)
 
 Compute the median of an entire array `v`, or, optionally,
-along the dimensions in `region`. For an even number of
+along the given dimensions. For an even number of
 elements no exact median element exists, so the result is
 equivalent to calculating mean of two median elements.
 
@@ -629,7 +644,11 @@ equivalent to calculating mean of two median elements.
     Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
     to represent missing values, and the [`skipmissing`](@ref) function to omit them.
 """
-median(v::AbstractArray, region) = mapslices(median!, v, region)
+median(v::AbstractArray; dims=:) = _median(v, dims)
+
+_median(v::AbstractArray, dims) = mapslices(median!, v, dims)
+
+_median(v::AbstractArray{T}, ::Colon) where {T} = median!(copyto!(Array{T,1}(uninitialized, _length(v)), v))
 
 # for now, use the R/S definition of quantile; may want variants later
 # see ?quantile in R -- this is type 7
diff --git a/base/twiceprecision.jl b/base/twiceprecision.jl
index 16bfd22fe6833..ddba86b809c1e 100644
--- a/base/twiceprecision.jl
+++ b/base/twiceprecision.jl
@@ -537,8 +537,8 @@ function sum(r::StepRangeLen)
     np, nn = l - r.offset, r.offset - 1  # positive, negative
     # To prevent overflow in sum(1:n), multiply its factors by the step
     sp, sn = sumpair(np), sumpair(nn)
-    tp = _prod(r.step, sp[1], sp[2])
-    tn = _prod(r.step, sn[1], sn[2])
+    tp = _tp_prod(r.step, sp[1], sp[2])
+    tn = _tp_prod(r.step, sn[1], sn[2])
     s_hi, s_lo = add12(tp.hi, -tn.hi)
     s_lo += tp.lo - tn.lo
     # Add in contributions of ref
@@ -692,11 +692,11 @@ narrow(::Type{Float64}) = Float32
 narrow(::Type{Float32}) = Float16
 narrow(::Type{Float16}) = Float16
 
-function _prod(t::TwicePrecision, x, y...)
+function _tp_prod(t::TwicePrecision, x, y...)
     @_inline_meta
-    _prod(t * x, y...)
+    _tp_prod(t * x, y...)
 end
-_prod(t::TwicePrecision) = t
+_tp_prod(t::TwicePrecision) = t
 <(x::TwicePrecision{T}, y::TwicePrecision{T}) where {T} =
     x.hi < y.hi || ((x.hi == y.hi) & (x.lo < y.lo))
 
diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index d9aa4238df152..f2d56d0e0fdb6 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -1012,7 +1012,8 @@ chol(A::SparseMatrixCSC) = error("Use cholfact() instead of chol() for sparse ma
 lu(A::SparseMatrixCSC) = error("Use lufact() instead of lu() for sparse matrices.")
 eig(A::SparseMatrixCSC) = error("Use IterativeEigensolvers.eigs() instead of eig() for sparse matrices.")
 
-function Base.cov(X::SparseMatrixCSC, vardim::Int=1; corrected::Bool=true)
+function Base.cov(X::SparseMatrixCSC; dims::Int=1, corrected::Bool=true)
+    vardim = dims
     a, b = size(X)
     n, p = vardim == 1 ? (a, b) : (b, a)
 
@@ -1024,7 +1025,7 @@ function Base.cov(X::SparseMatrixCSC, vardim::Int=1; corrected::Bool=true)
     out = Matrix(Base.unscaled_covzm(X, vardim))
 
     # Compute x̄
-    x̄ᵀ = mean(X, vardim)
+    x̄ᵀ = mean(X, dims=vardim)
 
     # Subtract n*x̄*x̄' from X'X
     @inbounds for j in 1:p, i in 1:p
diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl
index a0edc3fb13a14..be7c15c213eff 100644
--- a/stdlib/SparseArrays/test/sparse.jl
+++ b/stdlib/SparseArrays/test/sparse.jl
@@ -487,10 +487,10 @@ end
         for f in (sum, prod, minimum, maximum, var)
             farr = Array(arr)
             @test f(arr) ≈ f(farr)
-            @test f(arr, 1) ≈ f(farr, 1)
-            @test f(arr, 2) ≈ f(farr, 2)
-            @test f(arr, (1, 2)) ≈ [f(farr)]
-            @test isequal(f(arr, 3), f(farr, 3))
+            @test f(arr, dims=1) ≈ f(farr, dims=1)
+            @test f(arr, dims=2) ≈ f(farr, dims=2)
+            @test f(arr, dims=(1, 2)) ≈ [f(farr)]
+            @test isequal(f(arr, dims=3), f(farr, dims=3))
         end
     end
 
@@ -503,9 +503,9 @@ end
         # case where f(0) would throw
         @test f(x->sqrt(x-1), pA .+ 1) ≈ f(sqrt.(pA))
         # these actually throw due to #10533
-        # @test f(x->sqrt(x-1), pA .+ 1, 1) ≈ f(sqrt(pA), 1)
-        # @test f(x->sqrt(x-1), pA .+ 1, 2) ≈ f(sqrt(pA), 2)
-        # @test f(x->sqrt(x-1), pA .+ 1, 3) ≈ f(pA)
+        # @test f(x->sqrt(x-1), pA .+ 1, dims=1) ≈ f(sqrt(pA), dims=1)
+        # @test f(x->sqrt(x-1), pA .+ 1, dims=2) ≈ f(sqrt(pA), dims=2)
+        # @test f(x->sqrt(x-1), pA .+ 1, dims=3) ≈ f(pA)
     end
 
     @testset "empty cases" begin
@@ -516,16 +516,16 @@ end
         @test var(sparse(Int[])) === NaN
 
         for f in (sum, prod, var)
-            @test isequal(f(spzeros(0, 1), 1), f(Matrix{Int}(I, 0, 1), 1))
-            @test isequal(f(spzeros(0, 1), 2), f(Matrix{Int}(I, 0, 1), 2))
-            @test isequal(f(spzeros(0, 1), (1, 2)), f(Matrix{Int}(I, 0, 1), (1, 2)))
-            @test isequal(f(spzeros(0, 1), 3), f(Matrix{Int}(I, 0, 1), 3))
+            @test isequal(f(spzeros(0, 1), dims=1), f(Matrix{Int}(I, 0, 1), dims=1))
+            @test isequal(f(spzeros(0, 1), dims=2), f(Matrix{Int}(I, 0, 1), dims=2))
+            @test isequal(f(spzeros(0, 1), dims=(1, 2)), f(Matrix{Int}(I, 0, 1), dims=(1, 2)))
+            @test isequal(f(spzeros(0, 1), dims=3), f(Matrix{Int}(I, 0, 1), dims=3))
         end
         for f in (minimum, maximum, findmin, findmax)
-            @test_throws ArgumentError f(spzeros(0, 1), 1)
-            @test isequal(f(spzeros(0, 1), 2), f(Matrix{Int}(I, 0, 1), 2))
-            @test_throws ArgumentError f(spzeros(0, 1), (1, 2))
-            @test isequal(f(spzeros(0, 1), 3), f(Matrix{Int}(I, 0, 1), 3))
+            @test_throws ArgumentError f(spzeros(0, 1), dims=1)
+            @test isequal(f(spzeros(0, 1), dims=2), f(Matrix{Int}(I, 0, 1), dims=2))
+            @test_throws ArgumentError f(spzeros(0, 1), dims=(1, 2))
+            @test isequal(f(spzeros(0, 1), dims=3), f(Matrix{Int}(I, 0, 1), dims=3))
         end
     end
 end
@@ -582,9 +582,9 @@ end
     @test minimum(-P) === -3.0
     @test maximum(-P) === 0.0
 
-    @test maximum(P, (1,)) == [1.0 2.0 3.0]
-    @test maximum(P, (2,)) == reshape([1.0,2.0,3.0],3,1)
-    @test maximum(P, (1,2)) == reshape([3.0],1,1)
+    @test maximum(P, dims=(1,)) == [1.0 2.0 3.0]
+    @test maximum(P, dims=(2,)) == reshape([1.0,2.0,3.0],3,1)
+    @test maximum(P, dims=(1,2)) == reshape([3.0],1,1)
 
     @test maximum(sparse(fill(-1,3,3))) == -1
     @test minimum(sparse(fill(1,3,3))) == 1
@@ -1035,7 +1035,7 @@ end
     @test findmin(S) == findmin(A)
     @test findmax(S) == findmax(A)
     for region in [(1,), (2,), (1,2)], m in [findmax, findmin]
-        @test m(S, region) == m(A, region)
+        @test m(S, dims=region) == m(A, dims=region)
     end
 
     S = spzeros(10,8)
@@ -1114,33 +1114,33 @@ end
 
     A = sparse([BigInt(10)])
     for (tup, rval, rind) in [((2,), [BigInt(10)], [1])]
-        @test isequal(findmin(A, tup), (rval, rind))
+        @test isequal(findmin(A, dims=tup), (rval, rind))
     end
 
     for (tup, rval, rind) in [((2,), [BigInt(10)], [1])]
-        @test isequal(findmax(A, tup), (rval, rind))
+        @test isequal(findmax(A, dims=tup), (rval, rind))
     end
 
     A = sparse([BigInt(-10)])
     for (tup, rval, rind) in [((2,), [BigInt(-10)], [1])]
-        @test isequal(findmin(A, tup), (rval, rind))
+        @test isequal(findmin(A, dims=tup), (rval, rind))
     end
 
     for (tup, rval, rind) in [((2,), [BigInt(-10)], [1])]
-        @test isequal(findmax(A, tup), (rval, rind))
+        @test isequal(findmax(A, dims=tup), (rval, rind))
     end
 
     A = sparse([BigInt(10) BigInt(-10)])
     for (tup, rval, rind) in [((2,), reshape([BigInt(-10)], 1, 1), reshape([CartesianIndex(1,2)], 1, 1))]
-        @test isequal(findmin(A, tup), (rval, rind))
+        @test isequal(findmin(A, dims=tup), (rval, rind))
     end
 
     for (tup, rval, rind) in [((2,), reshape([BigInt(10)], 1, 1), reshape([CartesianIndex(1,1)], 1, 1))]
-        @test isequal(findmax(A, tup), (rval, rind))
+        @test isequal(findmax(A, dims=tup), (rval, rind))
     end
 
     A = sparse(["a", "b"])
-    @test_throws MethodError findmin(A, 1)
+    @test_throws MethodError findmin(A, dims=1)
 end
 
 # Support the case when user defined `zero` and `isless` for non-numerical type
@@ -1153,11 +1153,11 @@ Base.isless(x::CustomType, y::CustomType) = isless(x.x, y.x)
     A = sparse([CustomType("a"), CustomType("b")])
 
     for (tup, rval, rind) in [((1,), [CustomType("a")], [1])]
-        @test isequal(findmin(A, tup), (rval, rind))
+        @test isequal(findmin(A, dims=tup), (rval, rind))
     end
 
     for (tup, rval, rind) in [((1,), [CustomType("b")], [2])]
-        @test isequal(findmax(A, tup), (rval, rind))
+        @test isequal(findmax(A, dims=tup), (rval, rind))
     end
 end
 
@@ -2043,8 +2043,8 @@ end
     x_dense  = convert(Matrix{elty}, x_sparse)
     @testset "Test with no Infs and NaNs, vardim=$vardim, corrected=$corrected" for vardim in (1, 2),
                                                                                  corrected in (true, false)
-        @test cov(x_sparse, vardim, corrected=corrected) ≈
-              cov(x_dense , vardim, corrected=corrected)
+        @test cov(x_sparse, dims=vardim, corrected=corrected) ≈
+              cov(x_dense , dims=vardim, corrected=corrected)
     end
 
     @testset "Test with $x11, vardim=$vardim, corrected=$corrected" for x11 in (NaN, Inf),
@@ -2053,8 +2053,8 @@ end
         x_sparse[1,1] = x11
         x_dense[1 ,1] = x11
 
-        cov_sparse = cov(x_sparse, vardim, corrected=corrected)
-        cov_dense  = cov(x_dense , vardim, corrected=corrected)
+        cov_sparse = cov(x_sparse, dims=vardim, corrected=corrected)
+        cov_dense  = cov(x_dense , dims=vardim, corrected=corrected)
         @test cov_sparse[2:end, 2:end] ≈ cov_dense[2:end, 2:end]
         @test isfinite.(cov_sparse) == isfinite.(cov_dense)
         @test isfinite.(cov_sparse) == isfinite.(cov_dense)
@@ -2067,8 +2067,8 @@ end
         x_sparse[2,1] = NaN
         x_dense[2 ,1] = NaN
 
-        cov_sparse = cov(x_sparse, vardim, corrected=corrected)
-        cov_dense  = cov(x_dense , vardim, corrected=corrected)
+        cov_sparse = cov(x_sparse, dims=vardim, corrected=corrected)
+        cov_dense  = cov(x_dense , dims=vardim, corrected=corrected)
         @test cov_sparse[(1 + vardim):end, (1 + vardim):end] ≈
               cov_dense[ (1 + vardim):end, (1 + vardim):end]
         @test isfinite.(cov_sparse) == isfinite.(cov_dense)
@@ -2217,7 +2217,7 @@ end
 # #25943
 @testset "operations on Integer subtypes" begin
     s = sparse(UInt8[1, 2, 3], UInt8[1, 2, 3], UInt8[1, 2, 3])
-    @test sum(s, 2) == reshape([1, 2, 3], 3, 1)
+    @test sum(s, dims=2) == reshape([1, 2, 3], 3, 1)
 end
 
-end # module
\ No newline at end of file
+end # module
diff --git a/test/arrayops.jl b/test/arrayops.jl
index 75660e3ff03e9..fa94995359e3e 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -627,13 +627,13 @@ end
     A1 = reshape(repeat([1,2],1,12),2,3,4)
     A2 = reshape(repeat([1 2 3],2,4),2,3,4)
     A3 = reshape(repeat([1 2 3 4],6,1),2,3,4)
-    @test isequal(cumsum(A,1),A1)
-    @test isequal(cumsum(A,2),A2)
-    @test isequal(cumsum(A,3),A3)
+    @test isequal(cumsum(A,dims=1),A1)
+    @test isequal(cumsum(A,dims=2),A2)
+    @test isequal(cumsum(A,dims=3),A3)
 
     # issue 20112
     A3 = reshape(repeat([1 2 3 4],UInt32(6),UInt32(1)),2,3,4)
-    @test isequal(cumsum(A,3),A3)
+    @test isequal(cumsum(A,dims=3),A3)
     @test repeat([1,2,3,4], UInt32(1)) == [1,2,3,4]
     @test repeat([1 2], UInt32(2)) == repeat([1 2], UInt32(2), UInt32(1))
 
@@ -852,13 +852,13 @@ end
 
     A = rand(4,4)
     for s in Any[A[1:2:4, 1:2:4], view(A, 1:2:4, 1:2:4)]
-        c = cumsum(s, 1)
+        c = cumsum(s, dims=1)
         @test c[1,1] == A[1,1]
         @test c[2,1] == A[1,1]+A[3,1]
         @test c[1,2] == A[1,3]
         @test c[2,2] == A[1,3]+A[3,3]
 
-        c = cumsum(s, 2)
+        c = cumsum(s, dims=2)
         @test c[1,1] == A[1,1]
         @test c[2,1] == A[3,1]
         @test c[1,2] == A[1,1]+A[1,3]
@@ -902,7 +902,7 @@ end
 end
 
 # issue #2342
-@test isequal(cumsum([1 2 3], 1), [1 2 3])
+@test isequal(cumsum([1 2 3], dims=1), [1 2 3])
 
 @testset "set-like operations" begin
     @test isequal(union([1,2,3], [4,3,4]), [1,2,3,4])
@@ -1062,31 +1062,31 @@ end
     @test isless(asc[:,2],asc[:,1])
     @test isless(asc[:,3],asc[:,2])
 
-    as = sort(a, 1)
+    as = sort(a, dims=1)
     @test issorted(as[:,1])
     @test issorted(as[:,2])
     @test issorted(as[:,3])
 
-    as = sort(a, 2)
+    as = sort(a, dims=2)
     @test issorted(as[1,:])
     @test issorted(as[2,:])
     @test issorted(as[3,:])
 
     local b = rand(21,21,2)
 
-    bs = sort(b, 1)
+    bs = sort(b, dims=1)
     for i in 1:21
         @test issorted(bs[:,i,1])
         @test issorted(bs[:,i,2])
     end
 
-    bs = sort(b, 2)
+    bs = sort(b, dims=2)
     for i in 1:21
         @test issorted(bs[i,:,1])
         @test issorted(bs[i,:,2])
     end
 
-    bs = sort(b, 3)
+    bs = sort(b, dims=3)
     @test all(bs[:,:,1] .<= bs[:,:,2])
 end
 
@@ -1851,8 +1851,8 @@ copyto!(S, A)
 @test cat(1, A, B, S) == cat(1, A, A, A)
 @test cat(2, A, B, S) == cat(2, A, A, A)
 
-@test cumsum(A, 1) == cumsum(B, 1) == cumsum(S, 1)
-@test cumsum(A, 2) == cumsum(B, 2) == cumsum(S, 2)
+@test cumsum(A, dims=1) == cumsum(B, dims=1) == cumsum(S, dims=1)
+@test cumsum(A, dims=2) == cumsum(B, dims=2) == cumsum(S, dims=2)
 
 @test mapslices(sort, A, 1) == mapslices(sort, B, 1) == mapslices(sort, S, 1)
 @test mapslices(sort, A, 2) == mapslices(sort, B, 2) == mapslices(sort, S, 2)
@@ -1944,15 +1944,15 @@ let f = OOB_Functor([1,2])
 end
 
 # issue 15654
-@test cumprod([5], 2) == [5]
-@test cumprod([1 2; 3 4], 3) == [1 2; 3 4]
-@test cumprod([1 2; 3 4], 1) == [1 2; 3 8]
-@test cumprod([1 2; 3 4], 2) == [1 2; 3 12]
+@test cumprod([5], dims=2) == [5]
+@test cumprod([1 2; 3 4], dims=3) == [1 2; 3 4]
+@test cumprod([1 2; 3 4], dims=1) == [1 2; 3 8]
+@test cumprod([1 2; 3 4], dims=2) == [1 2; 3 12]
 
-@test cumsum([5], 2) == [5]
-@test cumsum([1 2; 3 4], 1) == [1 2; 4 6]
-@test cumsum([1 2; 3 4], 2) == [1 3; 3 7]
-@test cumsum([1 2; 3 4], 3) == [1 2; 3 4]
+@test cumsum([5], dims=2) == [5]
+@test cumsum([1 2; 3 4], dims=1) == [1 2; 4 6]
+@test cumsum([1 2; 3 4], dims=2) == [1 3; 3 7]
+@test cumsum([1 2; 3 4], dims=3) == [1 2; 3 4]
 
 # issue #18363
 @test_throws DimensionMismatch cumsum!([0,0], 1:4)
@@ -2111,8 +2111,8 @@ end
 
 @testset "accumulate, accumulate!" begin
     @test accumulate(+, [1,2,3]) == [1, 3, 6]
-    @test accumulate(min, [1 2; 3 4], 1) == [1 2; 1 2]
-    @test accumulate(max, [1 2; 3 0], 2) == [1 2; 3 3]
+    @test accumulate(min, [1 2; 3 4], dims=1) == [1 2; 1 2]
+    @test accumulate(max, [1 2; 3 0], dims=2) == [1 2; 3 3]
     @test accumulate(+, Bool[]) == Int[]
     @test accumulate(*, Bool[]) == Bool[]
     @test accumulate(+, Float64[]) == Float64[]
@@ -2120,10 +2120,10 @@ end
     @test accumulate(min, [1, 2, 5, -1, 3, -2]) == [1, 1, 1, -1, -1, -2]
     @test accumulate(max, [1, 2, 5, -1, 3, -2]) == [1, 2, 5, 5, 5, 5]
 
-    @test accumulate(max, [1 0; 0 1], 1) == [1 0; 1 1]
-    @test accumulate(max, [1 0; 0 1], 2) == [1 1; 0 1]
-    @test accumulate(min, [1 0; 0 1], 1) == [1 0; 0 0]
-    @test accumulate(min, [1 0; 0 1], 2) == [1 0; 0 0]
+    @test accumulate(max, [1 0; 0 1], dims=1) == [1 0; 1 1]
+    @test accumulate(max, [1 0; 0 1], dims=2) == [1 1; 0 1]
+    @test accumulate(min, [1 0; 0 1], dims=1) == [1 0; 0 0]
+    @test accumulate(min, [1 0; 0 1], dims=2) == [1 0; 0 0]
 
     @test isa(accumulate(+,     Int[]) , Vector{Int})
     @test isa(accumulate(+, 1., Int[]) , Vector{Float64})
@@ -2139,7 +2139,7 @@ end
             @test accumulate_arr ≈ cumop(arr)
             @test accumulate_arr[end] ≈ reduce(op, arr)
             @test accumulate_arr[1] ≈ arr[1]
-            @test accumulate(op, arr, 10) ≈ arr
+            @test accumulate(op, arr, dims=10) ≈ arr
 
             if eltype(arr) in [Int, Float64] # eltype of out easy
                 out = similar(arr)
@@ -2160,7 +2160,7 @@ end
     # asymmetric operation
     op(x,y) = 2x+y
     @test accumulate(op, [10,20, 30]) == [10, op(10, 20), op(op(10, 20), 30)] == [10, 40, 110]
-    @test accumulate(op, [10 20 30], 2) == [10 op(10, 20) op(op(10, 20), 30)] == [10 40 110]
+    @test accumulate(op, [10 20 30], dims=2) == [10 op(10, 20) op(op(10, 20), 30)] == [10 40 110]
 end
 
 struct F21666{T <: Base.ArithmeticStyle}
diff --git a/test/bitarray.jl b/test/bitarray.jl
index 168c8a600cb91..40ce442c452fa 100644
--- a/test/bitarray.jl
+++ b/test/bitarray.jl
@@ -13,9 +13,9 @@ tc(r1,r2) = false
 bitcheck(b::BitArray) = Test._check_bitarray_consistency(b)
 bitcheck(x) = true
 
-function check_bitop_call(ret_type, func, args...)
-    r1 = func(args...)
-    r2 = func(map(x->(isa(x, BitArray) ? Array(x) : x), args)...)
+function check_bitop_call(ret_type, func, args...; kwargs...)
+    r1 = func(args...; kwargs...)
+    r2 = func(map(x->(isa(x, BitArray) ? Array(x) : x), args)...; kwargs...)
     ret_type ≢ nothing && !isa(r1, ret_type) && @show ret_type, r1
     ret_type ≢ nothing && @test isa(r1, ret_type)
     @test tc(r1, r2)
@@ -1230,9 +1230,9 @@ end
     b1 = bitrand(s1, s2, s3, s4)
     m1 = 1
     m2 = 3
-    @check_bit_operation maximum(b1, (m1, m2)) BitArray{4}
-    @check_bit_operation minimum(b1, (m1, m2)) BitArray{4}
-    @check_bit_operation sum(b1, (m1, m2)) Array{Int,4}
+    @check_bit_operation maximum(b1, dims=(m1, m2)) BitArray{4}
+    @check_bit_operation minimum(b1, dims=(m1, m2)) BitArray{4}
+    @check_bit_operation sum(b1, dims=(m1, m2)) Array{Int,4}
 
     @check_bit_operation maximum(b1) Bool
     @check_bit_operation minimum(b1) Bool
diff --git a/test/broadcast.jl b/test/broadcast.jl
index f094027160e27..c3f5c8be4f81d 100644
--- a/test/broadcast.jl
+++ b/test/broadcast.jl
@@ -213,7 +213,7 @@ end
 # PR #17300: loop fusion
 @test (x->x+1).((x->x+2).((x->x+3).(1:10))) == 7:16
 let A = [sqrt(i)+j for i = 1:3, j=1:4]
-    @test atan2.(log.(A), sum(A,1)) == broadcast(atan2, broadcast(log, A), sum(A, 1))
+    @test atan2.(log.(A), sum(A, dims=1)) == broadcast(atan2, broadcast(log, A), sum(A, dims=1))
 end
 let x = sin.(1:10)
     @test atan2.((x->x+1).(x), (x->x+2).(x)) == broadcast(atan2, x.+1, x.+2)
diff --git a/test/offsetarray.jl b/test/offsetarray.jl
index de21878a1cf72..42e05f1d4c252 100644
--- a/test/offsetarray.jl
+++ b/test/offsetarray.jl
@@ -338,21 +338,21 @@ A = OffsetArray(rand(4,4), (-3,5))
 @test maximum(A) == maximum(parent(A))
 @test minimum(A) == minimum(parent(A))
 @test extrema(A) == extrema(parent(A))
-@test maximum(A, 1) == OffsetArray(maximum(parent(A), 1), (0,A.offsets[2]))
-@test maximum(A, 2) == OffsetArray(maximum(parent(A), 2), (A.offsets[1],0))
-@test maximum(A, 1:2) == maximum(parent(A), 1:2)
+@test maximum(A, dims=1) == OffsetArray(maximum(parent(A), dims=1), (0,A.offsets[2]))
+@test maximum(A, dims=2) == OffsetArray(maximum(parent(A), dims=2), (A.offsets[1],0))
+@test maximum(A, dims=1:2) == maximum(parent(A), dims=1:2)
 C = similar(A)
-cumsum!(C, A, 1)
-@test parent(C) == cumsum(parent(A), 1)
-@test parent(cumsum(A, 1)) == cumsum(parent(A), 1)
-cumsum!(C, A, 2)
-@test parent(C) == cumsum(parent(A), 2)
+cumsum!(C, A, dims=1)
+@test parent(C) == cumsum(parent(A), dims=1)
+@test parent(cumsum(A, dims=1)) == cumsum(parent(A), dims=1)
+cumsum!(C, A, dims=2)
+@test parent(C) == cumsum(parent(A), dims=2)
 R = similar(A, (1:1, 6:9))
 maximum!(R, A)
-@test parent(R) == maximum(parent(A), 1)
+@test parent(R) == maximum(parent(A), dims=1)
 R = similar(A, (-2:1, 1:1))
 maximum!(R, A)
-@test parent(R) == maximum(parent(A), 2)
+@test parent(R) == maximum(parent(A), dims=2)
 amin, iamin = findmin(A)
 pmin, ipmin = findmin(parent(A))
 @test amin == pmin
@@ -372,11 +372,11 @@ I = findall(!iszero, z)
 @test findall(x->x==0, h) == [2]
 @test mean(A_3_3) == median(A_3_3) == 5
 @test mean(x->2x, A_3_3) == 10
-@test mean(A_3_3, 1) == median(A_3_3, 1) == OffsetArray([2 5 8], (0,A_3_3.offsets[2]))
-@test mean(A_3_3, 2) == median(A_3_3, 2) == OffsetArray(reshape([4,5,6],(3,1)), (A_3_3.offsets[1],0))
+@test mean(A_3_3, dims=1) == median(A_3_3, dims=1) == OffsetArray([2 5 8], (0,A_3_3.offsets[2]))
+@test mean(A_3_3, dims=2) == median(A_3_3, dims=2) == OffsetArray(reshape([4,5,6],(3,1)), (A_3_3.offsets[1],0))
 @test var(A_3_3) == 7.5
-@test std(A_3_3, 1) == OffsetArray([1 1 1], (0,A_3_3.offsets[2]))
-@test std(A_3_3, 2) == OffsetArray(reshape([3,3,3], (3,1)), (A_3_3.offsets[1],0))
+@test std(A_3_3, dims=1) == OffsetArray([1 1 1], (0,A_3_3.offsets[2]))
+@test std(A_3_3, dims=2) == OffsetArray(reshape([3,3,3], (3,1)), (A_3_3.offsets[1],0))
 @test sum(OffsetArray(fill(1,3000), -1000)) == 3000
 
 @test vecnorm(v) ≈ vecnorm(parent(v))
@@ -412,8 +412,8 @@ v = OffsetArray(rand(8), (-2,))
 @test sort(v) == OffsetArray(sort(parent(v)), v.offsets)
 @test sortrows(A) == OffsetArray(sortrows(parent(A)), A.offsets)
 @test sortcols(A) == OffsetArray(sortcols(parent(A)), A.offsets)
-@test sort(A, 1) == OffsetArray(sort(parent(A), 1), A.offsets)
-@test sort(A, 2) == OffsetArray(sort(parent(A), 2), A.offsets)
+@test sort(A, dims=1) == OffsetArray(sort(parent(A), dims=1), A.offsets)
+@test sort(A, dims=2) == OffsetArray(sort(parent(A), dims=2), A.offsets)
 
 @test mapslices(sort, A, 1) == OffsetArray(mapslices(sort, parent(A), 1), A.offsets)
 @test mapslices(sort, A, 2) == OffsetArray(mapslices(sort, parent(A), 2), A.offsets)
diff --git a/test/reduce.jl b/test/reduce.jl
index e5071f2745d1d..ed5734dc879b2 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -222,9 +222,9 @@ A = circshift(reshape(1:24,2,3,4), (0,1,1))
 @test extrema(A,(1,3)) == reshape([(5,24),(1,20),(3,22)],1,3,1)
 @test extrema(A,(2,3)) == reshape([(1,23),(2,24)],2,1,1)
 @test extrema(A,(1,2,3)) == reshape([(1,24)],1,1,1)
-@test size(extrema(A,1)) == size(maximum(A,1))
-@test size(extrema(A,(1,2))) == size(maximum(A,(1,2)))
-@test size(extrema(A,(1,2,3))) == size(maximum(A,(1,2,3)))
+@test size(extrema(A,1)) == size(maximum(A,dims=1))
+@test size(extrema(A,(1,2))) == size(maximum(A,dims=(1,2)))
+@test size(extrema(A,(1,2,3))) == size(maximum(A,dims=(1,2,3)))
 
 # any & all
 
diff --git a/test/reducedim.jl b/test/reducedim.jl
index 53435f454d464..e50a66ae22201 100644
--- a/test/reducedim.jl
+++ b/test/reducedim.jl
@@ -51,14 +51,14 @@ for region in Any[
     fill!(r, -1.5)
     @test minimum!(abs, r, Areduc, init=false) ≈ fill!(r2, -1.5)
 
-    @test sum(Areduc, region) ≈ safe_sum(Areduc, region)
-    @test prod(Areduc, region) ≈ safe_prod(Areduc, region)
-    @test maximum(Areduc, region) ≈ safe_maximum(Areduc, region)
-    @test minimum(Areduc, region) ≈ safe_minimum(Areduc, region)
-    @test sum(abs, Areduc, region) ≈ safe_sumabs(Areduc, region)
-    @test sum(abs2, Areduc, region) ≈ safe_sumabs2(Areduc, region)
-    @test maximum(abs, Areduc, region) ≈ safe_maxabs(Areduc, region)
-    @test minimum(abs, Areduc, region) ≈ safe_minabs(Areduc, region)
+    @test sum(Areduc, dims=region) ≈ safe_sum(Areduc, region)
+    @test prod(Areduc, dims=region) ≈ safe_prod(Areduc, region)
+    @test maximum(Areduc, dims=region) ≈ safe_maximum(Areduc, region)
+    @test minimum(Areduc, dims=region) ≈ safe_minimum(Areduc, region)
+    @test sum(abs, Areduc, dims=region) ≈ safe_sumabs(Areduc, region)
+    @test sum(abs2, Areduc, dims=region) ≈ safe_sumabs2(Areduc, region)
+    @test maximum(abs, Areduc, dims=region) ≈ safe_maxabs(Areduc, region)
+    @test minimum(abs, Areduc, dims=region) ≈ safe_minabs(Areduc, region)
 end
 
 # Test reduction along first dimension; this is special-cased for
@@ -68,9 +68,9 @@ r = fill(NaN, map(length, Base.reduced_indices(axes(Breduc), 1)))
 @test sum!(r, Breduc) ≈ safe_sum(Breduc, 1)
 @test sum!(abs, r, Breduc) ≈ safe_sumabs(Breduc, 1)
 @test sum!(abs2, r, Breduc) ≈ safe_sumabs2(Breduc, 1)
-@test sum(Breduc, 1) ≈ safe_sum(Breduc, 1)
-@test sum(abs, Breduc, 1) ≈ safe_sumabs(Breduc, 1)
-@test sum(abs2, Breduc, 1) ≈ safe_sumabs2(Breduc, 1)
+@test sum(Breduc, dims=1) ≈ safe_sum(Breduc, 1)
+@test sum(abs, Breduc, dims=1) ≈ safe_sumabs(Breduc, 1)
+@test sum(abs2, Breduc, dims=1) ≈ safe_sumabs2(Breduc, 1)
 
 fill!(r, 4.2)
 @test sum!(r, Breduc, init=false) ≈ safe_sum(Breduc, 1) .+ 4.2
@@ -91,37 +91,35 @@ let R = [2]
     @test prod!(R, A, init=false) == [1440]
 
     # min/max
-    @test reducedim(max, A, 1) == [3 6]
-    @test reducedim(min, A, 2) == reshape([1,2,3], 3, 1)
+    @test reduce(max, A, dims=1) == [3 6]
+    @test reduce(min, A, dims=2) == reshape([1,2,3], 3, 1)
 end
 
 # Small integers
-@test @inferred(sum(Int8[1], 1)) == [1]
-@test @inferred(sum(UInt8[1], 1)) == [1]
+@test @inferred(sum(Int8[1], dims=1)) == [1]
+@test @inferred(sum(UInt8[1], dims=1)) == [1]
 
 # Complex types
-@test typeof(@inferred(sum([1.0+1.0im], 1))) == Vector{ComplexF64}
-@test typeof(@inferred(Base.sum(abs, [1.0+1.0im], 1))) == Vector{Float64}
-@test typeof(@inferred(Base.sum(abs2, [1.0+1.0im], 1))) == Vector{Float64}
-@test typeof(@inferred(prod([1.0+1.0im], 1))) == Vector{ComplexF64}
-@test typeof(@inferred(Base.prod(abs, [1.0+1.0im], 1))) == Vector{Float64}
-@test typeof(@inferred(Base.prod(abs2, [1.0+1.0im], 1))) == Vector{Float64}
+@test typeof(@inferred(sum([1.0+1.0im], dims=1))) == Vector{ComplexF64}
+@test typeof(@inferred(Base.sum(abs, [1.0+1.0im], dims=1))) == Vector{Float64}
+@test typeof(@inferred(Base.sum(abs2, [1.0+1.0im], dims=1))) == Vector{Float64}
+@test typeof(@inferred(prod([1.0+1.0im], dims=1))) == Vector{ComplexF64}
+@test typeof(@inferred(Base.prod(abs, [1.0+1.0im], dims=1))) == Vector{Float64}
+@test typeof(@inferred(Base.prod(abs2, [1.0+1.0im], dims=1))) == Vector{Float64}
 
 # Heterogeneously typed arrays
-@test sum(Union{Float32, Float64}[1.0], 1) == [1.0]
-@test prod(Union{Float32, Float64}[1.0], 1) == [1.0]
+@test sum(Union{Float32, Float64}[1.0], dims=1) == [1.0]
+@test prod(Union{Float32, Float64}[1.0], dims=1) == [1.0]
 
-@test reducedim((a,b) -> a|b, [true false; false false], 1, false) == [true false]
-let R = reducedim((a,b) -> a+b, [1 2; 3 4], 2, 0.0)
+@test reduce((a,b) -> a|b, false, [true false; false false], dims=1) == [true false]
+let R = reduce((a,b) -> a+b, 0.0, [1 2; 3 4], dims=2)
     @test eltype(R) == Float64
     @test R ≈ [3,7]
 end
-@test reducedim((a,b) -> a+b, [1 2; 3 4], 1, 0) == [4 6]
+@test reduce((a,b) -> a+b, 0, [1 2; 3 4], dims=1) == [4 6]
 
 # inferred return types
-let rt = Base.return_types(reducedim, Tuple{Function, Array{Float64, 3}, Int, Float64})
-    @test length(rt) == 1 && rt[1] == Array{Float64, 3}
-end
+@test typeof(@inferred(reduce(+, 0.0, ones(3,3,3), dims=1))) == Array{Float64, 3}
 
 @testset "empty cases" begin
     A = Matrix{Int}(uninitialized, 0,1)
@@ -131,30 +129,30 @@ end
     @test_throws ArgumentError maximum(A)
     @test var(A) === NaN
 
-    @test isequal(sum(A, 1), zeros(Int, 1, 1))
-    @test isequal(sum(A, 2), zeros(Int, 0, 1))
-    @test isequal(sum(A, (1, 2)), zeros(Int, 1, 1))
-    @test isequal(sum(A, 3), zeros(Int, 0, 1))
-    @test isequal(prod(A, 1), fill(1, 1, 1))
-    @test isequal(prod(A, 2), fill(1, 0, 1))
-    @test isequal(prod(A, (1, 2)), fill(1, 1, 1))
-    @test isequal(prod(A, 3), fill(1, 0, 1))
-    @test isequal(var(A, 1), fill(NaN, 1, 1))
-    @test isequal(var(A, 2), fill(NaN, 0, 1))
-    @test isequal(var(A, (1, 2)), fill(NaN, 1, 1))
-    @test isequal(var(A, 3), fill(NaN, 0, 1))
+    @test isequal(sum(A, dims=1), zeros(Int, 1, 1))
+    @test isequal(sum(A, dims=2), zeros(Int, 0, 1))
+    @test isequal(sum(A, dims=(1, 2)), zeros(Int, 1, 1))
+    @test isequal(sum(A, dims=3), zeros(Int, 0, 1))
+    @test isequal(prod(A, dims=1), fill(1, 1, 1))
+    @test isequal(prod(A, dims=2), fill(1, 0, 1))
+    @test isequal(prod(A, dims=(1, 2)), fill(1, 1, 1))
+    @test isequal(prod(A, dims=3), fill(1, 0, 1))
+    @test isequal(var(A, dims=1), fill(NaN, 1, 1))
+    @test isequal(var(A, dims=2), fill(NaN, 0, 1))
+    @test isequal(var(A, dims=(1, 2)), fill(NaN, 1, 1))
+    @test isequal(var(A, dims=3), fill(NaN, 0, 1))
 
     for f in (minimum, maximum)
-        @test_throws ArgumentError f(A, 1)
-        @test isequal(f(A, 2), zeros(Int, 0, 1))
-        @test_throws ArgumentError f(A, (1, 2))
-        @test isequal(f(A, 3), zeros(Int, 0, 1))
+        @test_throws ArgumentError f(A, dims=1)
+        @test isequal(f(A, dims=2), zeros(Int, 0, 1))
+        @test_throws ArgumentError f(A, dims=(1, 2))
+        @test isequal(f(A, dims=3), zeros(Int, 0, 1))
     end
     for f in (findmin, findmax)
-        @test_throws ArgumentError f(A, 1)
-        @test isequal(f(A, 2), (zeros(Int, 0, 1), zeros(Int, 0, 1)))
-        @test_throws ArgumentError f(A, (1, 2))
-        @test isequal(f(A, 3), (zeros(Int, 0, 1), zeros(Int, 0, 1)))
+        @test_throws ArgumentError f(A, dims=1)
+        @test isequal(f(A, dims=2), (zeros(Int, 0, 1), zeros(Int, 0, 1)))
+        @test_throws ArgumentError f(A, dims=(1, 2))
+        @test isequal(f(A, dims=3), (zeros(Int, 0, 1), zeros(Int, 0, 1)))
     end
 
 end
@@ -165,9 +163,9 @@ A = [1.0 5.0 6.0;
 for (tup, rval, rind) in [((1,), [1.0 2.0 4.0], [CartesianIndex(1,1) CartesianIndex(2,2) CartesianIndex(2,3)]),
                           ((2,), reshape([1.0,2.0], 2, 1), reshape([CartesianIndex(1,1),CartesianIndex(2,2)], 2, 1)),
                           ((1,2), fill(1.0,1,1),fill(CartesianIndex(1,1),1,1))]
-    @test findmin(A, tup) == (rval, rind)
+    @test findmin(A, dims=tup) == (rval, rind)
     @test findmin!(similar(rval), similar(rind), A) == (rval, rind)
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
@@ -175,9 +173,9 @@ end
 for (tup, rval, rind) in [((1,), [5.0 5.0 6.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(1,3)]),
                           ((2,), reshape([6.0,5.0], 2, 1), reshape([CartesianIndex(1,3),CartesianIndex(2,1)], 2, 1)),
                           ((1,2), fill(6.0,1,1),fill(CartesianIndex(1,3),1,1))]
-    @test findmax(A, tup) == (rval, rind)
+    @test findmax(A, dims=tup) == (rval, rind)
     @test findmax!(similar(rval), similar(rind), A) == (rval, rind)
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
@@ -189,9 +187,9 @@ A = [1.0 3.0 6.0;
 for (tup, rval, rind) in [((1,), [NaN 2.0 4.0], [CartesianIndex(2,1) CartesianIndex(2,2) CartesianIndex(2,3)]),
                           ((2,), reshape([1.0, NaN], 2, 1), reshape([CartesianIndex(1,1),CartesianIndex(2,1)], 2, 1)),
                           ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
     @test isequal(Base.reducedim!(min, copy(rval), A), rval)
@@ -200,9 +198,9 @@ end
 for (tup, rval, rind) in [((1,), [NaN 3.0 6.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(1,3)]),
                           ((2,), reshape([6.0, NaN], 2, 1), reshape([CartesianIndex(1,3),CartesianIndex(2,1)], 2, 1)),
                           ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
     @test isequal(Base.reducedim!(max, copy(rval), A), rval)
@@ -213,9 +211,9 @@ A = [1.0 NaN 6.0;
 for (tup, rval, rind) in [((1,), [NaN NaN 4.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(2,3)]),
                           ((2,), reshape([NaN, NaN], 2, 1), reshape([CartesianIndex(1,2),CartesianIndex(2,1)], 2, 1)),
                           ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
@@ -223,9 +221,9 @@ end
 for (tup, rval, rind) in [((1,), [NaN NaN 6.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(1,3)]),
                           ((2,), reshape([NaN, NaN], 2, 1), reshape([CartesianIndex(1,2),CartesianIndex(2,1)], 2, 1)),
                           ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
@@ -235,9 +233,9 @@ A = [Inf -Inf Inf  -Inf;
 for (tup, rval, rind) in [((1,), [Inf -Inf -Inf -Inf], [CartesianIndex(1,1) CartesianIndex(1,2) CartesianIndex(2,3) CartesianIndex(1,4)]),
                           ((2,), reshape([-Inf -Inf], 2, 1), reshape([CartesianIndex(1,2),CartesianIndex(2,3)], 2, 1)),
                           ((1,2), fill(-Inf,1,1),fill(CartesianIndex(1,2),1,1))]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
@@ -245,104 +243,104 @@ end
 for (tup, rval, rind) in [((1,), [Inf Inf Inf -Inf], [CartesianIndex(1,1) CartesianIndex(2,2) CartesianIndex(1,3) CartesianIndex(1,4)]),
                           ((2,), reshape([Inf Inf], 2, 1), reshape([CartesianIndex(1,1),CartesianIndex(2,1)], 2, 1)),
                           ((1,2), fill(Inf,1,1),fill(CartesianIndex(1,1),1,1))]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
 
 A = [BigInt(10)]
 for (tup, rval, rind) in [((2,), [BigInt(10)], [1])]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
 
 for (tup, rval, rind) in [((2,), [BigInt(10)], [1])]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
 
 A = [BigInt(-10)]
 for (tup, rval, rind) in [((2,), [BigInt(-10)], [1])]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
 
 for (tup, rval, rind) in [((2,), [BigInt(-10)], [1])]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
 
 A = [BigInt(10) BigInt(-10)]
 for (tup, rval, rind) in [((2,), reshape([BigInt(-10)], 1, 1), reshape([CartesianIndex(1,2)], 1, 1))]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
 
 for (tup, rval, rind) in [((2,), reshape([BigInt(10)], 1, 1), reshape([CartesianIndex(1,1)], 1, 1))]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
 
 A = ["a", "b"]
 for (tup, rval, rind) in [((1,), ["a"], [1])]
-    @test isequal(findmin(A, tup), (rval, rind))
+    @test isequal(findmin(A, dims=tup), (rval, rind))
     @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(minimum(A, tup), rval)
+    @test isequal(minimum(A, dims=tup), rval)
     @test isequal(minimum!(similar(rval), A), rval)
     @test isequal(minimum!(copy(rval), A, init=false), rval)
 end
 
 for (tup, rval, rind) in [((1,), ["b"], [2])]
-    @test isequal(findmax(A, tup), (rval, rind))
+    @test isequal(findmax(A, dims=tup), (rval, rind))
     @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
-    @test isequal(maximum(A, tup), rval)
+    @test isequal(maximum(A, dims=tup), rval)
     @test isequal(maximum!(similar(rval), A), rval)
     @test isequal(maximum!(copy(rval), A, init=false), rval)
 end
 
 # issue #6672
-@test sum(Real[1 2 3; 4 5.3 7.1], 2) == reshape([6, 16.4], 2, 1)
-@test std(AbstractFloat[1,2,3], 1) == [1.0]
-@test sum(Any[1 2;3 4],1) == [4 6]
-@test sum(Vector{Int}[[1,2],[4,3]], 1)[1] == [5,5]
+@test sum(Real[1 2 3; 4 5.3 7.1], dims=2) == reshape([6, 16.4], 2, 1)
+@test std(AbstractFloat[1,2,3], dims=1) == [1.0]
+@test sum(Any[1 2;3 4], dims=1) == [4 6]
+@test sum(Vector{Int}[[1,2],[4,3]], dims=1)[1] == [5,5]
 
 # issue #10461
 Areduc = rand(3, 4, 5, 6)
 for region in Any[-1, 0, (-1, 2), [0, 1], (1,-2,3), [0 1;
                                                      2 3], "hello"]
-    @test_throws ArgumentError sum(Areduc, region)
-    @test_throws ArgumentError prod(Areduc, region)
-    @test_throws ArgumentError maximum(Areduc, region)
-    @test_throws ArgumentError minimum(Areduc, region)
-    @test_throws ArgumentError sum(abs, Areduc, region)
-    @test_throws ArgumentError sum(abs2, Areduc, region)
-    @test_throws ArgumentError maximum(abs, Areduc, region)
-    @test_throws ArgumentError minimum(abs, Areduc, region)
+    @test_throws ArgumentError sum(Areduc, dims=region)
+    @test_throws ArgumentError prod(Areduc, dims=region)
+    @test_throws ArgumentError maximum(Areduc, dims=region)
+    @test_throws ArgumentError minimum(Areduc, dims=region)
+    @test_throws ArgumentError sum(abs, Areduc, dims=region)
+    @test_throws ArgumentError sum(abs2, Areduc, dims=region)
+    @test_throws ArgumentError maximum(abs, Areduc, dims=region)
+    @test_throws ArgumentError minimum(abs, Areduc, dims=region)
 end
 
 # check type of result
 @testset "type of sum(::Array{$T}" for T in [UInt8, Int8, Int32, Int64, BigInt]
-    result = sum(T[1 2 3; 4 5 6; 7 8 9], 2)
+    result = sum(T[1 2 3; 4 5 6; 7 8 9], dims=2)
     @test result == hcat([6, 15, 24])
     @test eltype(result) === (T <: Base.SmallSigned ? Int :
                               T <: Base.SmallUnsigned ? UInt :
diff --git a/test/statistics.jl b/test/statistics.jl
index eb65090780627..a9868c8b2e1ed 100644
--- a/test/statistics.jl
+++ b/test/statistics.jl
@@ -32,15 +32,15 @@ end
     @test isnan(median([-Inf,Inf]))
 
     X = [2 3 1 -1; 7 4 5 -4]
-    @test all(median(X, 2) .== [1.5, 4.5])
-    @test all(median(X, 1) .== [4.5 3.5 3.0 -2.5])
+    @test all(median(X, dims=2) .== [1.5, 4.5])
+    @test all(median(X, dims=1) .== [4.5 3.5 3.0 -2.5])
     @test X == [2 3 1 -1; 7 4 5 -4] # issue #17153
 
     @test_throws ArgumentError median([])
     @test isnan(median([NaN]))
     @test isnan(median([0.0,NaN]))
     @test isnan(median([NaN,0.0]))
-    @test isequal(median([NaN 0.0; 1.2 4.5], 2), reshape([NaN; 2.85], 2, 1))
+    @test isequal(median([NaN 0.0; 1.2 4.5], dims=2), reshape([NaN; 2.85], 2, 1))
 
     @test median!([1 2 3 4]) == 2.5
     @test median!([1 2; 3 4]) == 2.5
@@ -55,8 +55,8 @@ end
     @test mean([1.]) === 1.
     @test mean([1.,3]) == 2.
     @test mean([1,2,3]) == 2.
-    @test mean([0 1 2; 4 5 6], 1) == [2.  3.  4.]
-    @test mean([1 2 3; 4 5 6], 1) == [2.5 3.5 4.5]
+    @test mean([0 1 2; 4 5 6], dims=1) == [2.  3.  4.]
+    @test mean([1 2 3; 4 5 6], dims=1) == [2.5 3.5 4.5]
     @test mean(i->i+1, 0:2) === 2.
     @test mean(isodd, [3]) === 1.
     @test mean(x->3x, (1,1)) === 3.
@@ -68,7 +68,7 @@ end
     @test isnan(mean([0.,Inf,-Inf]))
     @test isnan(mean([1.,-1.,Inf,-Inf]))
     @test isnan(mean([-Inf,Inf]))
-    @test isequal(mean([NaN 0.0; 1.2 4.5], 2), reshape([NaN; 2.85], 2, 1))
+    @test isequal(mean([NaN 0.0; 1.2 4.5], dims=2), reshape([NaN; 2.85], 2, 1))
 
     # Check that small types are accumulated using wider type
     for T in (Int8, UInt8)
@@ -76,7 +76,7 @@ end
         g = (v for v in x)
         @test mean(x) == mean(g) == typemax(T)
         @test mean(identity, x) == mean(identity, g) == typemax(T)
-        @test mean(x, 2) == [typemax(T)]'
+        @test mean(x, dims=2) == [typemax(T)]'
     end
 end
 
@@ -93,10 +93,10 @@ end
     @test isnan(var(Int[]; mean=2))
     @test isnan(var(Int[]; mean=2, corrected=false))
     # reduction across dimensions
-    @test isequal(var(Int[], 1), [NaN])
-    @test isequal(var(Int[], 1; corrected=false), [NaN])
-    @test isequal(var(Int[], 1; mean=[2]), [NaN])
-    @test isequal(var(Int[], 1; mean=[2], corrected=false), [NaN])
+    @test isequal(var(Int[], dims=1), [NaN])
+    @test isequal(var(Int[], dims=1; corrected=false), [NaN])
+    @test isequal(var(Int[], dims=1; mean=[2]), [NaN])
+    @test isequal(var(Int[], dims=1; mean=[2], corrected=false), [NaN])
 
     # edge case: one-element vector
     # iterable
@@ -110,10 +110,10 @@ end
     @test var([1]; mean=2) === Inf
     @test var([1]; mean=2, corrected=false) === 1.0
     # reduction across dimensions
-    @test isequal(@inferred(var([1], 1)), [NaN])
-    @test var([1], 1; corrected=false) ≈ [0.0]
-    @test var([1], 1; mean=[2]) ≈ [Inf]
-    @test var([1], 1; mean=[2], corrected=false) ≈ [1.0]
+    @test isequal(@inferred(var([1], dims=1)), [NaN])
+    @test var([1], dims=1; corrected=false) ≈ [0.0]
+    @test var([1], dims=1; mean=[2]) ≈ [Inf]
+    @test var([1], dims=1; mean=[2], corrected=false) ≈ [1.0]
 
     @test var(1:8) == 6.
     @test varm(1:8,1) == varm(Vector(1:8),1)
@@ -146,8 +146,8 @@ end
     @test var((1,2,3); mean=0, corrected=false) ≈ 14.0/3
     @test_throws ArgumentError var((1,2,3); mean=())
 
-    @test var([1 2 3 4 5; 6 7 8 9 10], 2) ≈ [2.5 2.5]'
-    @test var([1 2 3 4 5; 6 7 8 9 10], 2; corrected=false) ≈ [2.0 2.0]'
+    @test var([1 2 3 4 5; 6 7 8 9 10], dims=2) ≈ [2.5 2.5]'
+    @test var([1 2 3 4 5; 6 7 8 9 10], dims=2; corrected=false) ≈ [2.0 2.0]'
 
     @test stdm([1,2,3], 2) ≈ 1.
     @test std([1,2,3]) ≈ 1.
@@ -161,8 +161,8 @@ end
     @test std((1,2,3); mean=0) ≈ sqrt(7.0)
     @test std((1,2,3); mean=0, corrected=false) ≈ sqrt(14.0/3)
 
-    @test std([1 2 3 4 5; 6 7 8 9 10], 2) ≈ sqrt.([2.5 2.5]')
-    @test std([1 2 3 4 5; 6 7 8 9 10], 2; corrected=false) ≈ sqrt.([2.0 2.0]')
+    @test std([1 2 3 4 5; 6 7 8 9 10], dims=2) ≈ sqrt.([2.5 2.5]')
+    @test std([1 2 3 4 5; 6 7 8 9 10], dims=2; corrected=false) ≈ sqrt.([2.0 2.0]')
 
     let A = ComplexF64[exp(i*im) for i in 1:10^4]
         @test varm(A, 0.) ≈ sum(map(abs2, A)) / (length(A) - 1)
@@ -170,10 +170,10 @@ end
     end
 
     @test var([1//1, 2//1]) isa Rational{Int}
-    @test var([1//1, 2//1], 1) isa Vector{Rational{Int}}
+    @test var([1//1, 2//1], dims=1) isa Vector{Rational{Int}}
 
     @test std([1//1, 2//1]) isa Float64
-    @test std([1//1, 2//1], 1) isa Vector{Float64}
+    @test std([1//1, 2//1], dims=1) isa Vector{Float64}
 end
 
 function safe_cov(x, y, zm::Bool, cr::Bool)
@@ -226,12 +226,12 @@ Y = [6.0  2.0;
         @test c ≈ Cxx[1,1]
         @inferred cov(x1, corrected=cr)
 
-        @test cov(X) == Base.covm(X, mean(X, 1))
+        @test cov(X) == Base.covm(X, mean(X, dims=1))
         C = zm ? Base.covm(X, 0, vd, corrected=cr) :
-                 cov(X, vd, corrected=cr)
+                 cov(X, dims=vd, corrected=cr)
         @test size(C) == (k, k)
         @test C ≈ Cxx
-        @inferred cov(X, vd, corrected=cr)
+        @inferred cov(X, dims=vd, corrected=cr)
 
         @test cov(x1, y1) == Base.covm(x1, mean(x1), y1, mean(y1))
         c = zm ? Base.covm(x1, 0, y1, 0, corrected=cr) :
@@ -241,29 +241,29 @@ Y = [6.0  2.0;
         @inferred cov(x1, y1, corrected=cr)
 
         if vd == 1
-            @test cov(x1, Y) == Base.covm(x1, mean(x1), Y, mean(Y, 1))
+            @test cov(x1, Y) == Base.covm(x1, mean(x1), Y, mean(Y, dims=1))
         end
         C = zm ? Base.covm(x1, 0, Y, 0, vd, corrected=cr) :
-                 cov(x1, Y, vd, corrected=cr)
+                 cov(x1, Y, dims=vd, corrected=cr)
         @test size(C) == (1, k)
         @test vec(C) ≈ Cxy[1,:]
-        @inferred cov(x1, Y, vd, corrected=cr)
+        @inferred cov(x1, Y, dims=vd, corrected=cr)
 
         if vd == 1
-            @test cov(X, y1) == Base.covm(X, mean(X, 1), y1, mean(y1))
+            @test cov(X, y1) == Base.covm(X, mean(X, dims=1), y1, mean(y1))
         end
         C = zm ? Base.covm(X, 0, y1, 0, vd, corrected=cr) :
-                 cov(X, y1, vd, corrected=cr)
+                 cov(X, y1, dims=vd, corrected=cr)
         @test size(C) == (k, 1)
         @test vec(C) ≈ Cxy[:,1]
-        @inferred cov(X, y1, vd, corrected=cr)
+        @inferred cov(X, y1, dims=vd, corrected=cr)
 
-        @test cov(X, Y) == Base.covm(X, mean(X, 1), Y, mean(Y, 1))
+        @test cov(X, Y) == Base.covm(X, mean(X, dims=1), Y, mean(Y, dims=1))
         C = zm ? Base.covm(X, 0, Y, 0, vd, corrected=cr) :
-                 cov(X, Y, vd, corrected=cr)
+                 cov(X, Y, dims=vd, corrected=cr)
         @test size(C) == (k, k)
         @test C ≈ Cxy
-        @inferred cov(X, Y, vd, corrected=cr)
+        @inferred cov(X, Y, dims=vd, corrected=cr)
     end
 end
 
@@ -306,11 +306,11 @@ end
         @test c ≈ Cxx[1,1]
         @inferred cor(x1)
 
-        @test cor(X) == Base.corm(X, mean(X, 1))
-        C = zm ? Base.corm(X, 0, vd) : cor(X, vd)
+        @test cor(X) == Base.corm(X, mean(X, dims=1))
+        C = zm ? Base.corm(X, 0, vd) : cor(X, dims=vd)
         @test size(C) == (k, k)
         @test C ≈ Cxx
-        @inferred cor(X, vd)
+        @inferred cor(X, dims=vd)
 
         @test cor(x1, y1) == Base.corm(x1, mean(x1), y1, mean(y1))
         c = zm ? Base.corm(x1, 0, y1, 0) : cor(x1, y1)
@@ -319,26 +319,26 @@ end
         @inferred cor(x1, y1)
 
         if vd == 1
-            @test cor(x1, Y) == Base.corm(x1, mean(x1), Y, mean(Y, 1))
+            @test cor(x1, Y) == Base.corm(x1, mean(x1), Y, mean(Y, dims=1))
         end
-        C = zm ? Base.corm(x1, 0, Y, 0, vd) : cor(x1, Y, vd)
+        C = zm ? Base.corm(x1, 0, Y, 0, vd) : cor(x1, Y, dims=vd)
         @test size(C) == (1, k)
         @test vec(C) ≈ Cxy[1,:]
-        @inferred cor(x1, Y, vd)
+        @inferred cor(x1, Y, dims=vd)
 
         if vd == 1
-            @test cor(X, y1) == Base.corm(X, mean(X, 1), y1, mean(y1))
+            @test cor(X, y1) == Base.corm(X, mean(X, dims=1), y1, mean(y1))
         end
-        C = zm ? Base.corm(X, 0, y1, 0, vd) : cor(X, y1, vd)
+        C = zm ? Base.corm(X, 0, y1, 0, vd) : cor(X, y1, dims=vd)
         @test size(C) == (k, 1)
         @test vec(C) ≈ Cxy[:,1]
-        @inferred cor(X, y1, vd)
+        @inferred cor(X, y1, dims=vd)
 
-        @test cor(X, Y) == Base.corm(X, mean(X, 1), Y, mean(Y, 1))
-        C = zm ? Base.corm(X, 0, Y, 0, vd) : cor(X, Y, vd)
+        @test cor(X, Y) == Base.corm(X, mean(X, dims=1), Y, mean(Y, dims=1))
+        C = zm ? Base.corm(X, 0, Y, 0, vd) : cor(X, Y, dims=vd)
         @test size(C) == (k, k)
         @test C ≈ Cxy
-        @inferred cor(X, Y, vd)
+        @inferred cor(X, Y, dims=vd)
     end
 
     @test cor(repeat(1:17, 1, 17))[2] <= 1.0
@@ -374,11 +374,11 @@ end
 
 @testset "variance of complex arrays (#13309)" begin
     z = rand(ComplexF64, 10)
-    @test var(z) ≈ invoke(var, Tuple{Any}, z) ≈ cov(z) ≈ var(z,1)[1] ≈ sum(abs2, z .- mean(z))/9
+    @test var(z) ≈ invoke(var, Tuple{Any}, z) ≈ cov(z) ≈ var(z,dims=1)[1] ≈ sum(abs2, z .- mean(z))/9
     @test isa(var(z), Float64)
     @test isa(invoke(var, Tuple{Any}, z), Float64)
     @test isa(cov(z), Float64)
-    @test isa(var(z,1), Vector{Float64})
+    @test isa(var(z,dims=1), Vector{Float64})
     @test varm(z, 0.0) ≈ invoke(varm, Tuple{Any,Float64}, z, 0.0) ≈ sum(abs2, z)/9
     @test isa(varm(z, 0.0), Float64)
     @test isa(invoke(varm, Tuple{Any,Float64}, z, 0.0), Float64)
@@ -410,22 +410,22 @@ end
 @testset "Issue #17153 and PR #17154" begin
     a = rand(10,10)
     b = deepcopy(a)
-    x = median(a, 1)
+    x = median(a, dims=1)
 
     @test b == a
-    x = median(a, 2)
+    x = median(a, dims=2)
     @test b == a
-    x = mean(a, 1)
+    x = mean(a, dims=1)
     @test b == a
-    x = mean(a, 2)
+    x = mean(a, dims=2)
     @test b == a
-    x = var(a, 1)
+    x = var(a, dims=1)
     @test b == a
-    x = var(a, 2)
+    x = var(a, dims=2)
     @test b == a
-    x = std(a, 1)
+    x = std(a, dims=1)
     @test b == a
-    x = std(a, 2)
+    x = std(a, dims=2)
     @test b == a
 end
 
@@ -443,9 +443,9 @@ using .Main.TestHelpers: Furlong
 
     # Issue #21786
     A = [Furlong{1}(rand(-5:5)) for i in 1:2, j in 1:2]
-    @test mean(mean(A, 1), 2)[1] === mean(A)
-    @test var(A, 1)[1] === var(A[:, 1])
-    @test std(A, 1)[1] === std(A[:, 1])
+    @test mean(mean(A, dims=1), dims=2)[1] === mean(A)
+    @test var(A, dims=1)[1] === var(A[:, 1])
+    @test std(A, dims=1)[1] === std(A[:, 1])
 end
 
 # Issue #22901
@@ -461,9 +461,9 @@ end
 
 @testset "Promotion in covzm. Issue #8080" begin
     A = [1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
-    @test Base.covzm(A) - mean(A, 1)'*mean(A, 1)*size(A, 1)/(size(A, 1) - 1) ≈ cov(A)
+    @test Base.covzm(A) - mean(A, dims=1)'*mean(A, dims=1)*size(A, 1)/(size(A, 1) - 1) ≈ cov(A)
     A = [1//1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
-    @test (A'A - size(A, 1)*Base.mean(A, 1)'*Base.mean(A, 1))/4 == cov(A)
+    @test (A'A - size(A, 1)*Base.mean(A, dims=1)'*Base.mean(A, dims=1))/4 == cov(A)
 end
 
 @testset "Mean along dimension of empty array" begin
@@ -471,15 +471,15 @@ end
     a00 = zeros(0, 0)
     a01 = zeros(0, 1)
     a10 = zeros(1, 0)
-    @test isequal(mean(a0, 1)      , fill(NaN, 1))
-    @test isequal(mean(a00, (1, 2)), fill(NaN, 1, 1))
-    @test isequal(mean(a01, 1)     , fill(NaN, 1, 1))
-    @test isequal(mean(a10, 2)     , fill(NaN, 1, 1))
+    @test isequal(mean(a0, dims=1)      , fill(NaN, 1))
+    @test isequal(mean(a00, dims=(1, 2)), fill(NaN, 1, 1))
+    @test isequal(mean(a01, dims=1)     , fill(NaN, 1, 1))
+    @test isequal(mean(a10, dims=2)     , fill(NaN, 1, 1))
 end
 
 @testset "cov/var/std of Vector{Vector}" begin
     x = [[2,4,6],[4,6,8]]
-    @test var(x) ≈ vec(var([x[1] x[2]], 2))
-    @test std(x) ≈ vec(std([x[1] x[2]], 2))
-    @test cov(x) ≈ cov([x[1] x[2]], 2)
+    @test var(x) ≈ vec(var([x[1] x[2]], dims=2))
+    @test std(x) ≈ vec(std([x[1] x[2]], dims=2))
+    @test cov(x) ≈ cov([x[1] x[2]], dims=2)
 end