diff --git a/NEWS.md b/NEWS.md
index f05535e65c7605..3df074f68e603c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -194,6 +194,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}}`,
@@ -388,11 +395,6 @@ Library improvements
     This supersedes the old behavior of reinterpret on Arrays. As a result, reinterpreting
     arrays with different alignment requirements (removed in 0.6) is once again allowed ([#23750]).
 
-  * `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
 -----------------------------
 
@@ -668,6 +670,8 @@ Deprecated or removed
   * `a:b` is deprecated for constructing a `StepRange` when `a` and `b` have physical units
     (Dates and Times). Use `a:s:b`, where `s = Dates.Day(1)` or `s = Dates.Second(1)`.
 
+  * `sub2ind` and `ind2sub` are deprecated in favor of using `CartesianRange` and `CartesianToLinear` ([#24715])
+
 Command-line option changes
 ---------------------------
 
@@ -1629,3 +1633,4 @@ Command-line option changes
 [#24320]: https://github.com/JuliaLang/julia/issues/24320
 [#24396]: https://github.com/JuliaLang/julia/issues/24396
 [#24413]: https://github.com/JuliaLang/julia/issues/24413
+[#24715]: https://github.com/JuliaLang/julia/issues/24715
diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index a030cf2dcf84bf..02e5bba45c8d62 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
@@ -811,17 +810,6 @@ if all inputs have fast linear indexing, a [`CartesianRange`](@ref)
 otherwise).
 If the arrays have different sizes and/or dimensionalities, `eachindex` returns an
 iterable that spans the largest range along each dimension.
-
-For a CartesianRange, this returns a reshaped range of the linear indices into
-the range, e.g.:
-
-```jldoctest
-julia> eachindex(CartesianRange((1:2,1:3)))
-2×3 reshape(::Base.OneTo{Int64}, 2, 3) with eltype Int64:
- 1  3  5
- 2  4  6
-```
-
 """
 eachindex(A::AbstractArray) = (@_inline_meta(); eachindex(IndexStyle(A), A))
 
@@ -955,7 +943,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
@@ -989,8 +977,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
@@ -1573,73 +1561,43 @@ function (==)(A::AbstractArray, B::AbstractArray)
     return true
 end
 
-# sub2ind and ind2sub
+# _sub2ind and _ind2sub
 # fallbacks
-function sub2ind(A::AbstractArray, I...)
+function _sub2ind(A::AbstractArray, I...)
     @_inline_meta
-    sub2ind(indices(A), I...)
+    _sub2ind(indices(A), I...)
 end
 
-"""
-    ind2sub(a, index) -> subscripts
-
-Returns 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(indices(A), ind)
+    _ind2sub(indices(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), returns the linear index corresponding to the provided subscripts.
+_sub2ind(::Tuple{}) = 1
+_sub2ind(::DimsInteger) = 1
+_sub2ind(::Indices) = 1
+_sub2ind(::Tuple{}, I::Integer...) = (@_inline_meta; _sub2ind_recurse((), 1, 1, I...))
 
-# Examples
-```jldoctest
-julia> sub2ind((5,6,7),1,2,3)
-66
-
-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
@@ -1647,42 +1605,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
-
-Returns 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
@@ -1691,12 +1630,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]
@@ -1712,21 +1651,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
@@ -1734,17 +1673,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 bd385c1e154648..943243a4792001 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 88c604a38d9234..bdd32f4c94c9df 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -2116,6 +2116,27 @@ end
 finalizer(f::Ptr{Void}, o::Ptr{Void}) = invoke(finalizer, Tuple{Ptr{Void}, Any}, f, o)
 finalizer(f::Ptr{Void}, o::Function) = invoke(finalizer, Tuple{Ptr{Void}, Any}, f, o)
 
+# 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 275fc68be1c175..8273516a8968c3 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,
     indices,
     indmax,
@@ -522,7 +522,6 @@ export
     step,
     stride,
     strides,
-    sub2ind,
     sum!,
     sum,
     sum_kbn,
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 3f4ca3c8815508..8aa716838289fd 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
@@ -211,22 +211,7 @@ module IteratorsMD
     CartesianIndex(1, 2)
     ```
 
-    For cartesian to linear index conversion, [`eachindex`](@ref) returns a
-    reshaped version of the linear indices when called on a `CartesianRange`:
-
-    ```jldoctest subarray
-    julia> linear = eachindex(cartesian)
-    3×2 reshape(::Base.OneTo{Int64}, 3, 2) with eltype Int64:
-     1  4
-     2  5
-     3  6
-
-    julia> linear[1,2]
-    4
-
-    julia> linear[cartesian[4]]
-    4
-    ```
+    For cartesian to linear index conversion, see [`CartesianToLinear`](@ref).
     """
     struct CartesianRange{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{CartesianIndex{N},N}
         indices::R
@@ -270,7 +255,6 @@ module IteratorsMD
     # 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)
-    Base.eachindex(iter::CartesianRange) = reshape(linearindices(iter), size(iter))
 
     ndims(R::CartesianRange) = ndims(typeof(R))
     ndims(::Type{CartesianRange{N}}) where {N} = N
@@ -397,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 b047275fbd82f0..6b885c54b79b45 100644
--- a/base/precompile.jl
+++ b/base/precompile.jl
@@ -1628,9 +1628,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 fef634061165bc..304b84786f604d 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 d9af41a9ab18b6..720f0b826552c7 100644
--- a/base/sparse/sparsematrix.jl
+++ b/base/sparse/sparsematrix.jl
@@ -1264,7 +1264,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->x!=0, S)
@@ -2272,10 +2272,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]
@@ -2777,7 +2777,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
@@ -3470,7 +3470,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 dea981211f70b4..356bfe834df3d5 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}(nnzB)
     nzvalB = Vector{Tv}(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 1e62e6abb3f75a..5171cf527bb8d4 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(indices(V), i))
+_pointer(V::SubArray, i::Int) = pointer(V, Base._ind2sub(indices(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 e369986126d6a2..ea0b3e7530042e 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 `indices(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/manual/metaprogramming.md b/doc/src/manual/metaprogramming.md
index b14ace62648009..e11f111d95e54f 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 383b886ccb8952..493dc073bc210f 100644
--- a/doc/src/stdlib/arrays.md
+++ b/doc/src/stdlib/arrays.md
@@ -46,8 +46,6 @@ Base.IndexStyle
 Base.conj!
 Base.stride
 Base.strides
-Base.ind2sub
-Base.sub2ind
 Base.LinAlg.checksquare
 ```
 
@@ -76,6 +74,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 db36f74bf7ad99..b21c94b29f3ef8 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
 
@@ -825,8 +826,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
@@ -896,5 +897,5 @@ end
     end
 
     @test CartesianRange(ones(2,3)) == CartesianRange((2,3))
-    @test eachindex(CartesianRange((2,3))) == [1 3 5; 2 4 6]
+    @test CartesianToLinear((2,3)) == [1 3 5; 2 4 6]
 end
diff --git a/test/arrayops.jl b/test/arrayops.jl
index f81d51103335e4..dde20a13102992 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1380,7 +1380,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 52b9bd0a3702c8..7e81d05b635a26 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