diff --git a/NEWS.md b/NEWS.md
index 46fa26c0721c7..c7dc89a38227c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -211,6 +211,13 @@ This section lists changes that do not have deprecation warnings.
     longer present. Use `first(R)` and `last(R)` to obtain
     start/stop. ([#20974])
 
+  * `CartesianRange` inherits from AbstractArray and construction with an
+    `AbstractArray` argument constructs the indices for that array. Consequently,
+    linear indexing can be used to provide linear-to-cartesian conversion ([#24715])
+
+  * The type `CartesianToLinear` has been added, providing conversion from
+    cartesian incices to linear indices using the normal indexing operation. ([#24715])
+
   * The `Diagonal`, `Bidiagonal`, `Tridiagonal` and `SymTridiagonal` type definitions have
     changed from `Diagonal{T}`, `Bidiagonal{T}`, `Tridiagonal{T}` and `SymTridiagonal{T}`
     to `Diagonal{T,V<:AbstractVector{T}}`, `Bidiagonal{T,V<:AbstractVector{T}}`,
@@ -755,6 +762,8 @@ Deprecated or removed
   * `cumsum`, `cumprod`, `accumulate`, and their mutating versions now require a `dim`
     argument instead of defaulting to using the first dimension ([#24684]).
 
+  * `sub2ind` and `ind2sub` are deprecated in favor of using `CartesianRange` and `CartesianToLinear` ([#24715]).
+
   * The `sum_kbn` and `cumsum_kbn` functions have been moved to the
     [KahanSummation](https://github.com/JuliaMath/KahanSummation.jl) package ([#24869]).
 
@@ -1748,5 +1757,6 @@ Command-line option changes
 [#24396]: https://github.com/JuliaLang/julia/issues/24396
 [#24413]: https://github.com/JuliaLang/julia/issues/24413
 [#24653]: https://github.com/JuliaLang/julia/issues/24653
+[#24715]: https://github.com/JuliaLang/julia/issues/24715
 [#24869]: https://github.com/JuliaLang/julia/issues/24869
 [#25021]: https://github.com/JuliaLang/julia/issues/25021
diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 4c1ad65d92dd3..da7aabe1d39da 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -308,9 +308,8 @@ type:
 The default is `IndexCartesian()`.
 
 Julia's internal indexing machinery will automatically (and invisibly)
-convert all indexing operations into the preferred style using
-[`sub2ind`](@ref) or [`ind2sub`](@ref). This allows users to access
-elements of your array using any indexing style, even when explicit
+convert all indexing operations into the preferred style. This allows users
+to access elements of your array using any indexing style, even when explicit
 methods have not been provided.
 
 If you define both styles of indexing for your `AbstractArray`, this
@@ -962,7 +961,7 @@ end
 _to_linear_index(A::AbstractArray, i::Int) = i
 _to_linear_index(A::AbstractVector, i::Int, I::Int...) = i
 _to_linear_index(A::AbstractArray) = 1
-_to_linear_index(A::AbstractArray, I::Int...) = (@_inline_meta; sub2ind(A, I...))
+_to_linear_index(A::AbstractArray, I::Int...) = (@_inline_meta; _sub2ind(A, I...))
 
 ## IndexCartesian Scalar indexing: Canonical method is full dimensionality of Ints
 function _getindex(::IndexCartesian, A::AbstractArray, I::Vararg{Int,M}) where M
@@ -996,8 +995,8 @@ _to_subscript_indices(A, J::Tuple, Jrem::Tuple) = J # already bounds-checked, sa
 _to_subscript_indices(A::AbstractArray{T,N}, I::Vararg{Int,N}) where {T,N} = I
 _remaining_size(::Tuple{Any}, t::Tuple) = t
 _remaining_size(h::Tuple, t::Tuple) = (@_inline_meta; _remaining_size(tail(h), tail(t)))
-_unsafe_ind2sub(::Tuple{}, i) = () # ind2sub may throw(BoundsError()) in this case
-_unsafe_ind2sub(sz, i) = (@_inline_meta; ind2sub(sz, i))
+_unsafe_ind2sub(::Tuple{}, i) = () # _ind2sub may throw(BoundsError()) in this case
+_unsafe_ind2sub(sz, i) = (@_inline_meta; _ind2sub(sz, i))
 
 ## Setindex! is defined similarly. We first dispatch to an internal _setindex!
 # function that allows dispatch on array storage
@@ -1584,73 +1583,43 @@ function (==)(A::AbstractArray, B::AbstractArray)
     return anymissing ? missing : true
 end
 
-# sub2ind and ind2sub
+# _sub2ind and _ind2sub
 # fallbacks
-function sub2ind(A::AbstractArray, I...)
+function _sub2ind(A::AbstractArray, I...)
     @_inline_meta
-    sub2ind(axes(A), I...)
+    _sub2ind(axes(A), I...)
 end
 
-"""
-    ind2sub(a, index) -> subscripts
-
-Return a tuple of subscripts into array `a` corresponding to the linear index `index`.
-
-# Examples
-```jldoctest
-julia> A = ones(5,6,7);
-
-julia> ind2sub(A,35)
-(5, 1, 2)
-
-julia> ind2sub(A,70)
-(5, 2, 3)
-```
-"""
-function ind2sub(A::AbstractArray, ind)
+function _ind2sub(A::AbstractArray, ind)
     @_inline_meta
-    ind2sub(axes(A), ind)
+    _ind2sub(axes(A), ind)
 end
 
 # 0-dimensional arrays and indexing with []
-sub2ind(::Tuple{}) = 1
-sub2ind(::DimsInteger) = 1
-sub2ind(::Indices) = 1
-sub2ind(::Tuple{}, I::Integer...) = (@_inline_meta; _sub2ind((), 1, 1, I...))
-# Generic cases
-
-"""
-    sub2ind(dims, i, j, k...) -> index
-
-The inverse of [`ind2sub`](@ref), return the linear index corresponding to the provided subscripts.
-
-# Examples
-```jldoctest
-julia> sub2ind((5,6,7),1,2,3)
-66
+_sub2ind(::Tuple{}) = 1
+_sub2ind(::DimsInteger) = 1
+_sub2ind(::Indices) = 1
+_sub2ind(::Tuple{}, I::Integer...) = (@_inline_meta; _sub2ind_recurse((), 1, 1, I...))
 
-julia> sub2ind((5,6,7),1,6,3)
-86
-```
-"""
-sub2ind(dims::DimsInteger, I::Integer...) = (@_inline_meta; _sub2ind(dims, 1, 1, I...))
-sub2ind(inds::Indices, I::Integer...) = (@_inline_meta; _sub2ind(inds, 1, 1, I...))
+# Generic cases
+_sub2ind(dims::DimsInteger, I::Integer...) = (@_inline_meta; _sub2ind_recurse(dims, 1, 1, I...))
+_sub2ind(inds::Indices, I::Integer...) = (@_inline_meta; _sub2ind_recurse(inds, 1, 1, I...))
 # In 1d, there's a question of whether we're doing cartesian indexing
 # or linear indexing. Support only the former.
-sub2ind(inds::Indices{1}, I::Integer...) =
+_sub2ind(inds::Indices{1}, I::Integer...) =
     throw(ArgumentError("Linear indexing is not defined for one-dimensional arrays"))
-sub2ind(inds::Tuple{OneTo}, I::Integer...) = (@_inline_meta; _sub2ind(inds, 1, 1, I...)) # only OneTo is safe
-sub2ind(inds::Tuple{OneTo}, i::Integer)    = i
+_sub2ind(inds::Tuple{OneTo}, I::Integer...) = (@_inline_meta; _sub2ind_recurse(inds, 1, 1, I...)) # only OneTo is safe
+_sub2ind(inds::Tuple{OneTo}, i::Integer)    = i
 
-_sub2ind(::Any, L, ind) = ind
-function _sub2ind(::Tuple{}, L, ind, i::Integer, I::Integer...)
+_sub2ind_recurse(::Any, L, ind) = ind
+function _sub2ind_recurse(::Tuple{}, L, ind, i::Integer, I::Integer...)
     @_inline_meta
-    _sub2ind((), L, ind+(i-1)*L, I...)
+    _sub2ind_recurse((), L, ind+(i-1)*L, I...)
 end
-function _sub2ind(inds, L, ind, i::Integer, I::Integer...)
+function _sub2ind_recurse(inds, L, ind, i::Integer, I::Integer...)
     @_inline_meta
     r1 = inds[1]
-    _sub2ind(tail(inds), nextL(L, r1), ind+offsetin(i, r1)*L, I...)
+    _sub2ind_recurse(tail(inds), nextL(L, r1), ind+offsetin(i, r1)*L, I...)
 end
 
 nextL(L, l::Integer) = L*l
@@ -1658,42 +1627,23 @@ nextL(L, r::AbstractUnitRange) = L*unsafe_length(r)
 offsetin(i, l::Integer) = i-1
 offsetin(i, r::AbstractUnitRange) = i-first(r)
 
-ind2sub(::Tuple{}, ind::Integer) = (@_inline_meta; ind == 1 ? () : throw(BoundsError()))
-
-"""
-    ind2sub(dims, index) -> subscripts
-
-Return a tuple of subscripts into an array with dimensions `dims`,
-corresponding to the linear index `index`.
-
-# Examples
-```jldoctest
-julia> ind2sub((3,4),2)
-(2, 1)
-
-julia> ind2sub((3,4),3)
-(3, 1)
-
-julia> ind2sub((3,4),4)
-(1, 2)
-```
-"""
-ind2sub(dims::DimsInteger, ind::Integer) = (@_inline_meta; _ind2sub(dims, ind-1))
-ind2sub(inds::Indices, ind::Integer)     = (@_inline_meta; _ind2sub(inds, ind-1))
-ind2sub(inds::Indices{1}, ind::Integer) =
+_ind2sub(::Tuple{}, ind::Integer) = (@_inline_meta; ind == 1 ? () : throw(BoundsError()))
+_ind2sub(dims::DimsInteger, ind::Integer) = (@_inline_meta; _ind2sub_recurse(dims, ind-1))
+_ind2sub(inds::Indices, ind::Integer)     = (@_inline_meta; _ind2sub_recurse(inds, ind-1))
+_ind2sub(inds::Indices{1}, ind::Integer) =
     throw(ArgumentError("Linear indexing is not defined for one-dimensional arrays"))
-ind2sub(inds::Tuple{OneTo}, ind::Integer) = (ind,)
+_ind2sub(inds::Tuple{OneTo}, ind::Integer) = (ind,)
 
-_ind2sub(::Tuple{}, ind) = (ind+1,)
-function _ind2sub(indslast::NTuple{1}, ind)
+_ind2sub_recurse(::Tuple{}, ind) = (ind+1,)
+function _ind2sub_recurse(indslast::NTuple{1}, ind)
     @_inline_meta
     (_lookup(ind, indslast[1]),)
 end
-function _ind2sub(inds, ind)
+function _ind2sub_recurse(inds, ind)
     @_inline_meta
     r1 = inds[1]
     indnext, f, l = _div(ind, r1)
-    (ind-l*indnext+f, _ind2sub(tail(inds), indnext)...)
+    (ind-l*indnext+f, _ind2sub_recurse(tail(inds), indnext)...)
 end
 
 _lookup(ind, d::Integer) = ind+1
@@ -1702,12 +1652,12 @@ _div(ind, d::Integer) = div(ind, d), 1, d
 _div(ind, r::AbstractUnitRange) = (d = unsafe_length(r); (div(ind, d), first(r), d))
 
 # Vectorized forms
-function sub2ind(inds::Indices{1}, I1::AbstractVector{T}, I::AbstractVector{T}...) where T<:Integer
+function _sub2ind(inds::Indices{1}, I1::AbstractVector{T}, I::AbstractVector{T}...) where T<:Integer
     throw(ArgumentError("Linear indexing is not defined for one-dimensional arrays"))
 end
-sub2ind(inds::Tuple{OneTo}, I1::AbstractVector{T}, I::AbstractVector{T}...) where {T<:Integer} =
+_sub2ind(inds::Tuple{OneTo}, I1::AbstractVector{T}, I::AbstractVector{T}...) where {T<:Integer} =
     _sub2ind_vecs(inds, I1, I...)
-sub2ind(inds::Union{DimsInteger,Indices}, I1::AbstractVector{T}, I::AbstractVector{T}...) where {T<:Integer} =
+_sub2ind(inds::Union{DimsInteger,Indices}, I1::AbstractVector{T}, I::AbstractVector{T}...) where {T<:Integer} =
     _sub2ind_vecs(inds, I1, I...)
 function _sub2ind_vecs(inds, I::AbstractVector...)
     I1 = I[1]
@@ -1723,21 +1673,21 @@ end
 function _sub2ind!(Iout, inds, Iinds, I)
     @_noinline_meta
     for i in Iinds
-        # Iout[i] = sub2ind(inds, map(Ij -> Ij[i], I)...)
+        # Iout[i] = _sub2ind(inds, map(Ij -> Ij[i], I)...)
         Iout[i] = sub2ind_vec(inds, i, I)
     end
     Iout
 end
 
-sub2ind_vec(inds, i, I) = (@_inline_meta; sub2ind(inds, _sub2ind_vec(i, I...)...))
+sub2ind_vec(inds, i, I) = (@_inline_meta; _sub2ind(inds, _sub2ind_vec(i, I...)...))
 _sub2ind_vec(i, I1, I...) = (@_inline_meta; (I1[i], _sub2ind_vec(i, I...)...))
 _sub2ind_vec(i) = ()
 
-function ind2sub(inds::Union{DimsInteger{N},Indices{N}}, ind::AbstractVector{<:Integer}) where N
+function _ind2sub(inds::Union{DimsInteger{N},Indices{N}}, ind::AbstractVector{<:Integer}) where N
     M = length(ind)
     t = ntuple(n->similar(ind),Val(N))
     for (i,idx) in pairs(IndexLinear(), ind)
-        sub = ind2sub(inds, idx)
+        sub = _ind2sub(inds, idx)
         for j = 1:N
             t[j][i] = sub[j]
         end
@@ -1745,17 +1695,6 @@ function ind2sub(inds::Union{DimsInteger{N},Indices{N}}, ind::AbstractVector{<:I
     t
 end
 
-function ind2sub!(sub::Array{T}, dims::Tuple{Vararg{T}}, ind::T) where T<:Integer
-    ndims = length(dims)
-    for i=1:ndims-1
-        ind2 = div(ind-1,dims[i])+1
-        sub[i] = ind - dims[i]*(ind2-1)
-        ind = ind2
-    end
-    sub[ndims] = ind
-    return sub
-end
-
 ## iteration utilities ##
 
 """
diff --git a/base/array.jl b/base/array.jl
index 04b3e90b37be6..3d80f6ad72cc1 100644
--- a/base/array.jl
+++ b/base/array.jl
@@ -135,7 +135,7 @@ sizeof(a::Array) = Core.sizeof(a)
 
 function isassigned(a::Array, i::Int...)
     @_inline_meta
-    ii = (sub2ind(size(a), i...) % UInt) - 1
+    ii = (_sub2ind(size(a), i...) % UInt) - 1
     @boundscheck ii < length(a) % UInt || return false
     ccall(:jl_array_isassigned, Cint, (Any, UInt), a, ii) == 1
 end
diff --git a/base/deprecated.jl b/base/deprecated.jl
index b617ec96ff0dd..202d97eb04602 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -2997,6 +2997,27 @@ end
 
 @deprecate merge!(repo::LibGit2.GitRepo, args...; kwargs...) LibGit2.merge!(repo, args...; kwargs...)
 
+# sub2ind and ind2sub deprecation (PR #24715)
+@deprecate ind2sub(A::AbstractArray, ind) CartesianRange(A)[ind]
+@deprecate ind2sub(::Tuple{}, ind::Integer) CartesianRange()[ind]
+@deprecate ind2sub(dims::Tuple{Vararg{Integer,N}} where N, ind::Integer) CartesianRange(dims)[ind]
+@deprecate ind2sub(inds::Tuple{Base.OneTo}, ind::Integer) CartesianRange(inds)[ind]
+@deprecate ind2sub(inds::Tuple{AbstractUnitRange}, ind::Integer) CartesianRange(inds)[ind]
+@deprecate ind2sub(inds::Tuple{Vararg{AbstractUnitRange,N}} where N, ind::Integer) CartesianRange(inds)[ind]
+@deprecate ind2sub(inds::Union{DimsInteger{N},Indices{N}}  where N, ind::AbstractVector{<:Integer}) CartesianRange(inds)[ind]
+
+@deprecate sub2ind(A::AbstractArray, I...) CartesianToLinear(A)[I...]
+@deprecate sub2ind(dims::Tuple{}) CartesianToLinear(dims)[]
+@deprecate sub2ind(dims::DimsInteger) CartesianToLinear(dims)[]
+@deprecate sub2ind(dims::Indices) CartesianToLinear(dims)[]
+@deprecate sub2ind(dims::Tuple{}, I::Integer...) CartesianToLinear(dims)[I...]
+@deprecate sub2ind(dims::DimsInteger, I::Integer...) CartesianToLinear(dims)[I...]
+@deprecate sub2ind(inds::Indices, I::Integer...) CartesianToLinear(inds)[I...]
+@deprecate sub2ind(inds::Tuple{OneTo}, I::Integer...) CartesianToLinear(inds)[I...]
+@deprecate sub2ind(inds::Tuple{OneTo}, i::Integer) CartesianToLinear(inds)[i]
+@deprecate sub2ind(inds::Tuple{OneTo}, I1::AbstractVector{T}, I::AbstractVector{T}...) where {T<:Integer} CartesianToLinear(inds)[CartesianIndex.(I1, I...)]
+@deprecate sub2ind(inds::Union{DimsInteger,Indices}, I1::AbstractVector{T}, I::AbstractVector{T}...) where {T<:Integer} CartesianToLinear(inds)[CartesianIndex.(I1, I...)]
+
 # 24490 - warnings and messages
 const log_info_to = Dict{Tuple{Union{Module,Void},Union{Symbol,Void}},IO}()
 const log_warn_to = Dict{Tuple{Union{Module,Void},Union{Symbol,Void}},IO}()
diff --git a/base/exports.jl b/base/exports.jl
index 3d2c355994b1d..316211868a9fc 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -39,6 +39,7 @@ export
     BufferStream,
     CartesianIndex,
     CartesianRange,
+    CartesianToLinear,
     Channel,
     Cmd,
     Colon,
@@ -453,7 +454,6 @@ export
     flipdim,
     hcat,
     hvcat,
-    ind2sub,
     indexin,
     indmax,
     indmin,
@@ -522,7 +522,6 @@ export
     step,
     stride,
     strides,
-    sub2ind,
     sum!,
     sum,
     to_indices,
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index f5f79853eacd1..4fd556539402b 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -11,7 +11,7 @@ module IteratorsMD
     using Base: IndexLinear, IndexCartesian, AbstractCartesianIndex, fill_to_length, tail
     using Base.Iterators: Reverse
 
-    export CartesianIndex, CartesianRange
+    export CartesianIndex, CartesianRange, CartesianToLinear
 
     """
         CartesianIndex(i, j, k...)   -> I
@@ -172,6 +172,11 @@ module IteratorsMD
     Consequently these can be useful for writing algorithms that
     work in arbitrary dimensions.
 
+        CartesianRange(A::AbstractArray) -> R
+
+    As a convenience, constructing a CartesianRange from an array makes a
+    range of its indices.
+
     # Examples
     ```jldoctest
     julia> foreach(println, CartesianRange((2, 2, 2)))
@@ -183,9 +188,32 @@ module IteratorsMD
     CartesianIndex(2, 1, 2)
     CartesianIndex(1, 2, 2)
     CartesianIndex(2, 2, 2)
+
+    julia> CartesianRange(ones(2,3))
+    2×3 CartesianRange{2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}:
+      CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
+      CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
+    ```
+
+    ## Conversion between linear and cartesian indices
+
+    Linear index to cartesian index conversion exploits the fact that a
+    `CartesianRange` is an `AbstractArray` and can be indexed linearly:
+
+    ```jldoctest subarray
+    julia> cartesian = CartesianRange(1:3,1:2)
+    3×2 CartesianRange{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:
+     CartesianIndex(1, 1)  CartesianIndex(1, 2)
+     CartesianIndex(2, 1)  CartesianIndex(2, 2)
+     CartesianIndex(3, 1)  CartesianIndex(3, 2)
+
+    julia> cartesian[4]
+    CartesianIndex(1, 2)
     ```
+
+    For cartesian to linear index conversion, see [`CartesianToLinear`](@ref).
     """
-    struct CartesianRange{N,R<:NTuple{N,AbstractUnitRange{Int}}}
+    struct CartesianRange{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{CartesianIndex{N},N}
         indices::R
     end
 
@@ -204,6 +232,8 @@ module IteratorsMD
     CartesianRange(inds::NTuple{N,Union{<:Integer,AbstractUnitRange{<:Integer}}}) where {N} =
         CartesianRange(map(i->first(i):last(i), inds))
 
+    CartesianRange(A::AbstractArray) = CartesianRange(indices(A))
+
     convert(::Type{Tuple{}}, R::CartesianRange{0}) = ()
     convert(::Type{NTuple{N,AbstractUnitRange{Int}}}, R::CartesianRange{N}) where {N} =
         R.indices
@@ -222,6 +252,10 @@ module IteratorsMD
     convert(::Type{Tuple{Vararg{UnitRange}}}, R::CartesianRange) =
         convert(Tuple{Vararg{UnitRange{Int}}}, R)
 
+    # AbstractArray implementation
+    Base.IndexStyle(::Type{CartesianRange{N,R}}) where {N,R} = IndexCartesian()
+    @inline Base.getindex(iter::CartesianRange{N,R}, I::Vararg{Int, N}) where {N,R} = CartesianIndex(first.(iter.indices) .- 1 .+ I)
+
     ndims(R::CartesianRange) = ndims(typeof(R))
     ndims(::Type{CartesianRange{N}}) where {N} = N
     ndims(::Type{CartesianRange{N,TT}}) where {N,TT} = N
@@ -347,6 +381,53 @@ module IteratorsMD
     start(iter::Reverse{<:CartesianRange{0}}) = false
     next(iter::Reverse{<:CartesianRange{0}}, state) = CartesianIndex(), true
     done(iter::Reverse{<:CartesianRange{0}}, state) = state
+
+    """
+        CartesianToLinear(inds::CartesianRange) -> R
+        CartesianToLinear(sz::Dims) -> R
+        CartesianToLinear(istart:istop, jstart:jstop, ...) -> R
+
+    Define a mapping between cartesian indices and the corresponding linear index into a CartesianRange
+
+    # Example
+
+    The main purpose of this type is intuitive conversion from cartesian to linear indexing:
+
+    ```jldoctest subarray
+    julia> linear = CartesianToLinear(1:3,1:2)
+    CartesianToLinear{2,Tuple{UnitRange{Int64},UnitRange{Int64}}} with indices 1:3×1:2:
+      1  4
+      2  5
+      3  6
+
+    julia> linear[1,2]
+    4
+    ```
+    """
+    struct CartesianToLinear{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{Int,N}
+        indices::R
+    end
+
+    CartesianToLinear(inds::CartesianRange{N,R}) where {N,R} = CartesianToLinear{N,R}(inds.indices)
+    CartesianToLinear(::Tuple{}) = CartesianToLinear(CartesianRange(()))
+    CartesianToLinear(inds::NTuple{N,AbstractUnitRange{Int}}) where {N} = CartesianToLinear(CartesianRange(inds))
+    CartesianToLinear(inds::Vararg{AbstractUnitRange{Int},N}) where {N} = CartesianToLinear(CartesianRange(inds))
+    CartesianToLinear(inds::NTuple{N,AbstractUnitRange{<:Integer}}) where {N} = CartesianToLinear(CartesianRange(inds))
+    CartesianToLinear(inds::Vararg{AbstractUnitRange{<:Integer},N}) where {N} = CartesianToLinear(CartesianRange(inds))
+    CartesianToLinear(index::CartesianIndex) = CartesianToLinear(CartesianRange(index))
+    CartesianToLinear(sz::NTuple{N,<:Integer}) where {N} = CartesianToLinear(CartesianRange(sz))
+    CartesianToLinear(inds::NTuple{N,Union{<:Integer,AbstractUnitRange{<:Integer}}}) where {N} = CartesianToLinear(CartesianRange(inds))
+    CartesianToLinear(A::AbstractArray) = CartesianToLinear(CartesianRange(A))
+
+    # AbstractArray implementation
+    Base.IndexStyle(::Type{CartesianToLinear{N,R}}) where {N,R} = IndexCartesian()
+    Base.indices(iter::CartesianToLinear{N,R}) where {N,R} = iter.indices
+    @inline function Base.getindex(iter::CartesianToLinear{N,R}, I::Vararg{Int, N}) where {N,R}
+        dims = length.(iter.indices)
+        #without the inbounds, this is slower than Base._sub2ind(iter.indices, I...)
+        @inbounds result = reshape(1:prod(dims), dims)[(I .- first.(iter.indices) .+ 1)...]
+        return result
+    end
 end  # IteratorsMD
 
 
diff --git a/base/precompile.jl b/base/precompile.jl
index 5d5bbea92f4a7..ebc5313d40aa7 100644
--- a/base/precompile.jl
+++ b/base/precompile.jl
@@ -1367,9 +1367,9 @@ precompile(Tuple{typeof(Base.print_matrix_vdots), Base.IOContext{Base.Terminals.
 precompile(Tuple{typeof(Base.print_matrix), Base.IOContext{Base.Terminals.TTYTerminal}, Array{Int64, 1}, String, String, String, String, String, String, Int64, Int64})
 precompile(Tuple{typeof(Base.alignment), Base.IOContext{Base.Terminals.TTYTerminal}, Array{Int64, 1}, Array{Int64, 1}, Array{Int64, 1}, Int64, Int64, Int64})
 precompile(Tuple{getfield(Base, Symbol("#kw##sprint")), Array{Any, 1}, typeof(Base.sprint), Int64, typeof(Base.show), Int64})
-precompile(Tuple{typeof(Base.sub2ind), Tuple{Int64}, Int64, Int64})
-precompile(Tuple{typeof(Base._sub2ind), Tuple{Int64}, Int64, Int64, Int64, Int64})
-precompile(Tuple{typeof(Base._sub2ind), Tuple{}, Int64, Int64, Int64})
+precompile(Tuple{typeof(Base._sub2ind), Tuple{Int64}, Int64, Int64})
+precompile(Tuple{typeof(Base._sub2ind_recurse), Tuple{Int64}, Int64, Int64, Int64, Int64})
+precompile(Tuple{typeof(Base._sub2ind_recurse), Tuple{}, Int64, Int64, Int64})
 precompile(Tuple{typeof(Base.first), Array{Int64, 1}})
 precompile(Tuple{typeof(Base.print_matrix_row), Base.IOContext{Base.Terminals.TTYTerminal}, Array{Int64, 1}, Array{Tuple{Int64, Int64}, 1}, Int64, Array{Int64, 1}, String})
 precompile(Tuple{typeof(Base.print), Base.GenericIOBuffer{Array{UInt8, 1}}, Base.OneTo{Int64}})
diff --git a/base/reshapedarray.jl b/base/reshapedarray.jl
index e4f7daa2c5645..cef74e1eb6582 100644
--- a/base/reshapedarray.jl
+++ b/base/reshapedarray.jl
@@ -204,7 +204,7 @@ end
 end
 
 @inline function _unsafe_getindex(A::ReshapedArray{T,N}, indexes::Vararg{Int,N}) where {T,N}
-    i = sub2ind(size(A), indexes...)
+    i = Base._sub2ind(size(A), indexes...)
     I = ind2sub_rs(A.mi, i)
     _unsafe_getindex_rs(parent(A), I)
 end
@@ -227,7 +227,7 @@ end
 end
 
 @inline function _unsafe_setindex!(A::ReshapedArray{T,N}, val, indexes::Vararg{Int,N}) where {T,N}
-    @inbounds parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indexes...))...] = val
+    @inbounds parent(A)[ind2sub_rs(A.mi, Base._sub2ind(size(A), indexes...))...] = val
     val
 end
 
diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl
index f0a80e8d9da62..e05dbc2482d7a 100644
--- a/base/sparse/sparsematrix.jl
+++ b/base/sparse/sparsematrix.jl
@@ -1279,7 +1279,7 @@ function find(p::Function, S::SparseMatrixCSC)
     end
     sz = size(S)
     I, J = _findn(p, S)
-    return sub2ind(sz, I, J)
+    return Base._sub2ind(sz, I, J)
 end
 
 findn(S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = _findn(x->true, S)
@@ -2273,10 +2273,10 @@ function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractArray) where {Tv,Ti}
 
     for i in 1:n
         ((I[i] < 1) | (I[i] > nA)) && throw(BoundsError())
-        row,col = ind2sub(szA, I[i])
+        row,col = Base._ind2sub(szA, I[i])
         for r in colptrA[col]:(colptrA[col+1]-1)
             @inbounds if rowvalA[r] == row
-                rowB,colB = ind2sub(szB, i)
+                rowB,colB = Base._ind2sub(szB, i)
                 colptrB[colB+1] += 1
                 rowvalB[idxB] = rowB
                 nzvalB[idxB] = nzvalA[r]
@@ -2778,7 +2778,7 @@ function setindex!(A::SparseMatrixCSC, x, I::AbstractVector{<:Real})
         sxidx = S[xidx]
         (sxidx < n) && (I[sxidx] == I[sxidx+1]) && continue
 
-        row,col = ind2sub(szA, I[sxidx])
+        row,col = Base._ind2sub(szA, I[sxidx])
         v = isa(x, AbstractArray) ? x[sxidx] : x
 
         if col > lastcol
@@ -3471,7 +3471,7 @@ function hash(A::SparseMatrixCSC{T}, h::UInt) where T
         for j = colptr[col]:colptr[col+1]-1
             nz = nzval[j]
             isequal(nz, zero(T)) && continue
-            idx = sub2ind(sz, rowval[j], col)
+            idx = Base._sub2ind(sz, rowval[j], col)
             if idx != lastidx+1 || !isequal(nz, lastnz)  # Run is over
                 h = hashrun(lastnz, runlength, h)        # Hash previous run
                 h = hashrun(0, idx-lastidx-1, h)         # Hash intervening zeros
diff --git a/base/sparse/sparsevector.jl b/base/sparse/sparsevector.jl
index 293a9d5925954..cdb758f09a3fa 100644
--- a/base/sparse/sparsevector.jl
+++ b/base/sparse/sparsevector.jl
@@ -618,8 +618,8 @@ function getindex(A::SparseMatrixCSC{Tv}, I::AbstractUnitRange) where Tv
     rowvalB = Vector{Int}(uninitialized, nnzB)
     nzvalB = Vector{Tv}(uninitialized, nnzB)
 
-    rowstart,colstart = ind2sub(szA, first(I))
-    rowend,colend = ind2sub(szA, last(I))
+    rowstart,colstart = Base._ind2sub(szA, first(I))
+    rowend,colend = Base._ind2sub(szA, last(I))
 
     idxB = 1
     @inbounds for col in colstart:colend
@@ -628,7 +628,7 @@ function getindex(A::SparseMatrixCSC{Tv}, I::AbstractUnitRange) where Tv
         for r in colptrA[col]:(colptrA[col+1]-1)
             rowA = rowvalA[r]
             if minrow <= rowA <= maxrow
-                rowvalB[idxB] = sub2ind(szA, rowA, col) - first(I) + 1
+                rowvalB[idxB] = Base._sub2ind(szA, rowA, col) - first(I) + 1
                 nzvalB[idxB] = nzvalA[r]
                 idxB += 1
             end
@@ -656,7 +656,7 @@ function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector) where {Tv,Ti}
     idxB = 1
     for i in 1:n
         ((I[i] < 1) | (I[i] > nA)) && throw(BoundsError(A, I))
-        row,col = ind2sub(szA, I[i])
+        row,col = Base._ind2sub(szA, I[i])
         for r in colptrA[col]:(colptrA[col+1]-1)
             @inbounds if rowvalA[r] == row
                 if idxB <= nnzB
diff --git a/base/subarray.jl b/base/subarray.jl
index 9625161781a38..4d60e48262fe3 100644
--- a/base/subarray.jl
+++ b/base/subarray.jl
@@ -329,7 +329,7 @@ pointer(V::FastSubArray, i::Int) = pointer(V.parent, V.offset1 + V.stride1*i)
 pointer(V::FastContiguousSubArray, i::Int) = pointer(V.parent, V.offset1 + i)
 pointer(V::SubArray, i::Int) = _pointer(V, i)
 _pointer(V::SubArray{<:Any,1}, i::Int) = pointer(V, (i,))
-_pointer(V::SubArray, i::Int) = pointer(V, ind2sub(axes(V), i))
+_pointer(V::SubArray, i::Int) = pointer(V, Base._ind2sub(axes(V), i))
 
 function pointer(V::SubArray{T,N,<:Array,<:Tuple{Vararg{RangeIndex}}}, is::Tuple{Vararg{Int}}) where {T,N}
     index = first_index(V)
diff --git a/doc/src/devdocs/offset-arrays.md b/doc/src/devdocs/offset-arrays.md
index 4e6f77224faa1..c26ac5870372e 100644
--- a/doc/src/devdocs/offset-arrays.md
+++ b/doc/src/devdocs/offset-arrays.md
@@ -81,7 +81,7 @@ Some algorithms are most conveniently (or efficiently) written in terms of a sin
 
 For this reason, your best option may be to iterate over the array with `eachindex(A)`, or, if you require the indices to be sequential integers, to get the index range by calling `linearindices(A)`. This will return `axes(A, 1)` if A is an AbstractVector, and the equivalent of `1:length(A)` otherwise.
 
-By this definition, 1-dimensional arrays always use Cartesian indexing with the array's native indices. To help enforce this, it's worth noting that sub2ind(shape, i...) and ind2sub(shape, ind) will throw an error if shape indicates a 1-dimensional array with unconventional indexing (i.e., is a `Tuple{UnitRange}` rather than a tuple of `OneTo`). For arrays with conventional indexing, these functions continue to work the same as always.
+By this definition, 1-dimensional arrays always use Cartesian indexing with the array's native indices. To help enforce this, it's worth noting that the index conversion functions will throw an error if shape indicates a 1-dimensional array with unconventional indexing (i.e., is a `Tuple{UnitRange}` rather than a tuple of `OneTo`). For arrays with conventional indexing, these functions continue to work the same as always.
 
 Using `indices` and `linearindices`, here is one way you could rewrite `mycopy!`:
 
diff --git a/doc/src/devdocs/subarrays.md b/doc/src/devdocs/subarrays.md
index d564b0b1369a1..62421f289584f 100644
--- a/doc/src/devdocs/subarrays.md
+++ b/doc/src/devdocs/subarrays.md
@@ -22,9 +22,8 @@ computation (such as interpolation), and the type under discussion here, `SubArr
 For these types, the underlying information is more naturally described in terms of
 cartesian indexes.
 
-You can manually convert from a cartesian index to a linear index with `sub2ind`, and vice versa
-using `ind2sub`.  `getindex` and `setindex!` functions for `AbstractArray` types may include similar
-operations.
+The `getindex` and `setindex!` functions for `AbstractArray` types may include automatic conversion
+between indexing types. For explicit conversion, [`CartesianRange`](@ref) can be used.
 
 While converting from a cartesian index to a linear index is fast (it's just multiplication and
 addition), converting from a linear index to a cartesian index is very slow: it relies on the
diff --git a/doc/src/manual/metaprogramming.md b/doc/src/manual/metaprogramming.md
index b14ace6264800..e11f111d95e54 100644
--- a/doc/src/manual/metaprogramming.md
+++ b/doc/src/manual/metaprogramming.md
@@ -1264,7 +1264,7 @@ to build some more advanced (and valid) functionality...
 
 ### An advanced example
 
-Julia's base library has a [`sub2ind`](@ref) function to calculate a linear index into an n-dimensional
+Julia's base library has a an internal `sub2ind` function to calculate a linear index into an n-dimensional
 array, based on a set of n multilinear indices - in other words, to calculate the index `i` that
 can be used to index into an array `A` using `A[i]`, instead of `A[x,y,z,...]`. One possible implementation
 is the following:
diff --git a/doc/src/stdlib/arrays.md b/doc/src/stdlib/arrays.md
index 81d4d7d62e9f9..e8c1e841e7da1 100644
--- a/doc/src/stdlib/arrays.md
+++ b/doc/src/stdlib/arrays.md
@@ -48,8 +48,6 @@ Base.IndexStyle
 Base.conj!
 Base.stride
 Base.strides
-Base.ind2sub
-Base.sub2ind
 Base.LinAlg.checksquare
 ```
 
@@ -89,6 +87,7 @@ Base.isassigned
 Base.Colon
 Base.CartesianIndex
 Base.CartesianRange
+Base.CartesianToLinear
 Base.to_indices
 Base.checkbounds
 Base.checkindex
diff --git a/test/abstractarray.jl b/test/abstractarray.jl
index b74b8bc6b15fa..e1a8b40eecd68 100644
--- a/test/abstractarray.jl
+++ b/test/abstractarray.jl
@@ -112,62 +112,64 @@ end
     @test checkbounds(Bool, A, [CartesianIndex((5, 4))], 4) == false
 end
 
-@testset "sub2ind & ind2sub" begin
+@testset "index conversion" begin
     @testset "0-dimensional" begin
-        for i = 1:4
-            @test sub2ind((), i) == i
-        end
-        @test sub2ind((), 2, 2) == 3
-        @test ind2sub((), 1) == ()
-        @test_throws BoundsError ind2sub((), 2)
+        @test CartesianToLinear()[1] == 1
+        @test_throws BoundsError CartesianToLinear()[2]
+        @test CartesianToLinear()[1,1] == 1
+        @test CartesianRange()[1] == CartesianIndex()
+        @test_throws BoundsError CartesianRange()[2]
     end
 
     @testset "1-dimensional" begin
-        for i = 1:4
-            @test sub2ind((3,), i) == i
-            @test ind2sub((3,), i) == (i,)
+        for i = 1:3
+            @test CartesianToLinear((3,))[i] == i
+            @test CartesianRange((3,))[i] == CartesianIndex(i,)
         end
-        @test sub2ind((3,), 2, 2) == 5
-        @test_throws MethodError ind2sub((3,), 2, 2)
+        @test CartesianToLinear((3,))[2,1] == 2
+        @test_throws BoundsError CartesianRange((3,))[2,2]
         #   ambiguity btw cartesian indexing and linear indexing in 1d when
         #   indices may be nontraditional
-        @test_throws ArgumentError sub2ind((1:3,), 2)
-        @test_throws ArgumentError ind2sub((1:3,), 2)
+        @test_throws ArgumentError Base._sub2ind((1:3,), 2)
+        @test_throws ArgumentError Base._ind2sub((1:3,), 2)
     end
 
     @testset "2-dimensional" begin
         k = 0
+        cartesian = CartesianRange((4,3))
+        linear = CartesianToLinear(cartesian)
         for j = 1:3, i = 1:4
-            @test sub2ind((4,3), i, j) == (k+=1)
-            @test ind2sub((4,3), k) == (i,j)
-            @test sub2ind((1:4,1:3), i, j) == k
-            @test ind2sub((1:4,1:3), k) == (i,j)
-            @test sub2ind((0:3,3:5), i-1, j+2) == k
-            @test ind2sub((0:3,3:5), k) == (i-1, j+2)
+            @test linear[i,j] == (k+=1)
+            @test cartesian[k] == CartesianIndex(i,j)
+            @test CartesianToLinear(0:3,3:5)[i-1,j+2] == k
+            @test CartesianRange(0:3,3:5)[k] == CartesianIndex(i-1,j+2)
         end
     end
 
     @testset "3-dimensional" begin
         l = 0
         for k = 1:2, j = 1:3, i = 1:4
-            @test sub2ind((4,3,2), i, j, k) == (l+=1)
-            @test ind2sub((4,3,2), l) == (i,j,k)
-            @test sub2ind((1:4,1:3,1:2), i, j, k) == l
-            @test ind2sub((1:4,1:3,1:2), l) == (i,j,k)
-            @test sub2ind((0:3,3:5,-101:-100), i-1, j+2, k-102) == l
-            @test ind2sub((0:3,3:5,-101:-100), l) == (i-1, j+2, k-102)
+            @test CartesianToLinear((4,3,2))[i,j,k] == (l+=1)
+            @test CartesianRange((4,3,2))[l] == CartesianIndex(i,j,k)
+            @test CartesianToLinear(1:4,1:3,1:2)[i,j,k] == l
+            @test CartesianRange(1:4,1:3,1:2)[l] == CartesianIndex(i,j,k)
+            @test CartesianToLinear(0:3,3:5,-101:-100)[i-1,j+2,k-102] == l
+            @test CartesianRange(0:3,3:5,-101:-100)[l] == CartesianIndex(i-1, j+2, k-102)
         end
 
         local A = reshape(collect(1:9), (3,3))
-        @test ind2sub(size(A), 6) == (3,2)
-        @test sub2ind(size(A), 3, 2) == 6
-        @test ind2sub(A, 6) == (3,2)
-        @test sub2ind(A, 3, 2) == 6
+        @test CartesianRange(size(A))[6] == CartesianIndex(3,2)
+        @test CartesianToLinear(size(A))[3, 2] == 6
+        @test CartesianRange(A)[6] == CartesianIndex(3,2)
+        @test CartesianToLinear(A)[3, 2] == 6
+        for i in 1:length(A)
+            @test CartesianToLinear(A)[CartesianRange(A)[i]] == i
+        end
 
         @testset "PR #9256" begin
             function pr9256()
                 m = [1 2 3; 4 5 6; 7 8 9]
-                ind2sub(m, 6)
+                Base._ind2sub(m, 6)
             end
             @test pr9256() == (3,2)
         end
@@ -616,10 +618,9 @@ function test_ind2sub(::Type{TestAbstractArray})
     dims = tuple(rand(1:5, n)...)
     len = prod(dims)
     A = reshape(collect(1:len), dims...)
-    I = ind2sub(dims, [1:len...])
+    I = CartesianRange(dims)
     for i in 1:len
-        idx = [ I[j][i] for j in 1:n ]
-        @test A[idx...] == A[i]
+        @test A[I[i]] == A[i]
     end
 end
 
@@ -762,8 +763,8 @@ end
     @test Base.copymutable((1,2,3)) == [1,2,3]
 end
 
-@testset "sub2ind for empty tuple" begin
-    @test sub2ind(()) == 1
+@testset "_sub2ind for empty tuple" begin
+    @test Base._sub2ind(()) == 1
 end
 
 @testset "to_shape" begin
@@ -815,6 +816,27 @@ end
     @test Z == collect(Z) == copy(Z)
 end
 
+@testset "CartesianRange" begin
+    xrng = 2:4
+    yrng = 1:5
+    CR = CartesianRange((xrng,yrng))
+
+    for (i,i_idx) in enumerate(xrng)
+        for (j,j_idx) in enumerate(yrng)
+            @test CR[i,j] == CartesianIndex(i_idx,j_idx)
+        end
+    end
+
+    for i_lin in linearindices(CR)
+        i = (i_lin-1) % length(xrng) + 1
+        j = (i_lin-i) ÷ length(xrng) + 1
+        @test CR[i_lin] == CartesianIndex(xrng[i],yrng[j])
+    end
+
+    @test CartesianRange(ones(2,3)) == CartesianRange((2,3))
+    @test CartesianToLinear((2,3)) == [1 3 5; 2 4 6]
+end
+
 @testset "empty" begin
     @test isempty([])
     v = [1, 2, 3]
diff --git a/test/arrayops.jl b/test/arrayops.jl
index be890b7dec886..e03392c7b8351 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1358,7 +1358,7 @@ end
 # issue #7197
 function i7197()
     S = [1 2 3; 4 5 6; 7 8 9]
-    ind2sub(size(S), 5)
+    Base._ind2sub(size(S), 5)
 end
 @test i7197() == (2,2)
 
diff --git a/test/euler.jl b/test/euler.jl
index 52b9bd0a3702c..7e81d05b635a2 100644
--- a/test/euler.jl
+++ b/test/euler.jl
@@ -65,13 +65,14 @@ end
 #11: 70600674
 function euler11(grid,n)
     m = typemin(eltype(grid))
+    tolinear = CartesianToLinear(size(grid))
     for i = n:size(grid,1)-n+1,
         j = n:size(grid,2)-n+1,
         di = -1:1, dj = -1:1
         di == dj == 0 && continue
-        idx = sub2ind(size(grid),
-                      di==0 ? fill(i,n) : range(i,di,n),
-                      dj==0 ? fill(j,n) : range(j,dj,n))
+        i_idxs = di==0 ? fill(i,n) : range(i,di,n)
+        j_idxs = dj==0 ? fill(j,n) : range(j,dj,n)
+        idx = tolinear[CartesianIndex.(i_idxs, j_idxs)]
         m = max(m,prod(grid[idx]))
     end
     return m