diff --git a/NEWS.md b/NEWS.md
index 1baa1294bcd2c..0a6feb53ced2d 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}}`,
@@ -445,6 +452,11 @@ Library improvements
     defined, linear-algebra function `transpose`. Similarly,
     `permutedims(v::AbstractVector)` will create a row matrix ([#24839]).
 
+  * `CartesianRange` changes ([#24715]):
+    - Inherits from `AbstractArray`
+    - Constructor taking an array
+    - `eachindex` returns the linear indices into a reshaped array, as `sub2ind` alternative
+
 Compiler/Runtime improvements
 -----------------------------
 
@@ -792,6 +804,8 @@ Deprecated or removed
     and `unsafe_get`/`get` can be dropped or replaced with `coalesce`.
     `NullException` has been removed.
 
+  * `sub2ind` and `ind2sub` are deprecated in favor of using `CartesianRange` and `CartesianToLinear` ([#24715]).
+
 Command-line option changes
 ---------------------------
 
@@ -1768,6 +1782,7 @@ Command-line option changes
 [#24413]: https://github.com/JuliaLang/julia/issues/24413
 [#24653]: https://github.com/JuliaLang/julia/issues/24653
 [#24714]: https://github.com/JuliaLang/julia/issues/24714
+[#24715]: https://github.com/JuliaLang/julia/issues/24715
 [#24869]: https://github.com/JuliaLang/julia/issues/24869
 [#25021]: https://github.com/JuliaLang/julia/issues/25021
 [#25088]: https://github.com/JuliaLang/julia/issues/25088
diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 9e1e9bd370d2f..c393f3fed70d7 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
@@ -959,7 +958,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
@@ -993,8 +992,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
@@ -1581,73 +1580,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
@@ -1655,42 +1624,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
@@ -1699,12 +1649,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]
@@ -1720,21 +1670,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
@@ -1742,17 +1692,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 c3ee01b4a4ad3..289382bc7ae04 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 87e29088041a6..558d2b7acb08f 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -3262,6 +3262,27 @@ end
 @deprecate_moved isnull "Nullables"
 @deprecate_moved unsafe_get "Nullables"
 
+# 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...)]
+
 # END 0.7 deprecations
 
 # BEGIN 1.0 deprecations
diff --git a/base/exports.jl b/base/exports.jl
index 29b110cebfdc7..0eb27e0866f27 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -39,6 +39,7 @@ export
     BufferStream,
     CartesianIndex,
     CartesianRange,
+    CartesianToLinear,
     Channel,
     Cmd,
     Colon,
@@ -452,7 +453,6 @@ export
     flipdim,
     hcat,
     hvcat,
-    ind2sub,
     indexin,
     indmax,
     indmin,
@@ -521,7 +521,6 @@ export
     step,
     stride,
     strides,
-    sub2ind,
     sum!,
     sum,
     to_indices,
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 5e6039ab08a57..3b008a6699b75 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(axes(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
@@ -343,6 +377,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.axes(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 1be887564e958..fec501e2bf697 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 070d25ffc20b7..005be471253b4 100644
--- a/base/reshapedarray.jl
+++ b/base/reshapedarray.jl
@@ -204,7 +204,7 @@ end
 end
 
 @inline function _unsafe_getindex(A::ReshapedArray{T,N}, indices::Vararg{Int,N}) where {T,N}
-    i = sub2ind(size(A), indices...)
+    i = Base._sub2ind(size(A), indices...)
     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, indices::Vararg{Int,N}) where {T,N}
-    @inbounds parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indices...))...] = val
+    @inbounds parent(A)[ind2sub_rs(A.mi, Base._sub2ind(size(A), indices...))...] = val
     val
 end
 
diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl
index 4d726ab4bf581..11ac69ab72b59 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)
@@ -2249,10 +2249,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]
@@ -2754,7 +2754,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
@@ -3447,7 +3447,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 c08b1cd2eeab1..83af34955370e 100644
--- a/base/sparse/sparsevector.jl
+++ b/base/sparse/sparsevector.jl
@@ -617,8 +617,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
@@ -627,7 +627,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
@@ -655,7 +655,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 97afeddeaaf17..3dd49f66a7d44 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 ac69ef0f151c9..62fc64413a39a 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 indices.
 
-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 ec4bdd2b4f1c4..b546a62a79b37 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 e894850f4f23d..2a5bf5f6d69d1 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
@@ -826,3 +827,24 @@ end
     @test isempty(v2::Vector{Int})
     @test isempty(v3::Vector{Float64})
 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
diff --git a/test/arrayops.jl b/test/arrayops.jl
index a49eef100f142..2165bcdef2eb4 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