diff --git a/NEWS.md b/NEWS.md
index ae4013cb098e6..56239fb1782be 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -83,6 +83,8 @@ New library features
 
 * The `redirect_*` functions can now be called on `IOContext` objects.
 * New constructor `NamedTuple(iterator)` that constructs a named tuple from a key-value pair iterator.
+* A new `reinterpret(reshape, T, a::AbstractArray{S})` reinterprets `a` to have eltype `T` while potentially
+  inserting or consuming the first dimension depending on the ratio of `sizeof(T)` and `sizeof(S)`.
 
 Standard library changes
 ------------------------
diff --git a/base/reinterpretarray.jl b/base/reinterpretarray.jl
index 4a11a4d55ae7d..18df68093441f 100644
--- a/base/reinterpretarray.jl
+++ b/base/reinterpretarray.jl
@@ -3,36 +3,39 @@
 """
 Gives a reinterpreted view (of element type T) of the underlying array (of element type S).
 If the size of `T` differs from the size of `S`, the array will be compressed/expanded in
-the first dimension.
+the first dimension. The variant `reinterpret(reshape, T, a)` instead adds or consumes the first dimension
+depending on the ratio of element sizes.
 """
-struct ReinterpretArray{T,N,S,A<:AbstractArray{S, N}} <: AbstractArray{T, N}
+struct ReinterpretArray{T,N,S,A<:AbstractArray{S},IsReshaped} <: AbstractArray{T, N}
     parent::A
     readable::Bool
     writable::Bool
+
+    function throwbits(S::Type, T::Type, U::Type)
+        @_noinline_meta
+        throw(ArgumentError("cannot reinterpret `$(S)` as `$(T)`, type `$(U)` is not a bits type"))
+    end
+    function throwsize0(S::Type, T::Type, msg)
+        @_noinline_meta
+        throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a $msg size"))
+    end
+
     global reinterpret
     function reinterpret(::Type{T}, a::A) where {T,N,S,A<:AbstractArray{S, N}}
-        function throwbits(::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U}
-            @_noinline_meta
-            throw(ArgumentError("cannot reinterpret `$(S)` `$(T)`, type `$(U)` is not a bits type"))
-        end
-        function throwsize0(::Type{S}, ::Type{T})
-            @_noinline_meta
-            throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a different size"))
-        end
-        function thrownonint(::Type{S}, ::Type{T}, dim)
+        function thrownonint(S::Type, T::Type, dim)
             @_noinline_meta
             throw(ArgumentError("""
                 cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
                 The resulting array would have non-integral first dimension.
                 """))
         end
-        function throwaxes1(::Type{S}, ::Type{T}, ax1)
+        function throwaxes1(S::Type, T::Type, ax1)
             @_noinline_meta
             throw(ArgumentError("cannot reinterpret a `$(S)` array to `$(T)` when the first axis is $ax1. Try reshaping first."))
         end
         isbitstype(T) || throwbits(S, T, T)
         isbitstype(S) || throwbits(S, T, S)
-        (N != 0 || sizeof(T) == sizeof(S)) || throwsize0(S, T)
+        (N != 0 || sizeof(T) == sizeof(S)) || throwsize0(S, T, "different")
         if N != 0 && sizeof(S) != sizeof(T)
             ax1 = axes(a)[1]
             dim = length(ax1)
@@ -41,15 +44,82 @@ struct ReinterpretArray{T,N,S,A<:AbstractArray{S, N}} <: AbstractArray{T, N}
         end
         readable = array_subpadding(T, S)
         writable = array_subpadding(S, T)
-        new{T, N, S, A}(a, readable, writable)
+        new{T, N, S, A, false}(a, readable, writable)
+    end
+
+    # With reshaping
+    function reinterpret(::typeof(reshape), ::Type{T}, a::A) where {T,S,A<:AbstractArray{S}}
+        function throwintmult(S::Type, T::Type)
+            @_noinline_meta
+            throw(ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got $(sizeof(T))) and `sizeof(eltype(a))` (got $(sizeof(S))) be an integer multiple of the other"))
+        end
+        function throwsize1(a::AbstractArray, T::Type)
+            @_noinline_meta
+            throw(ArgumentError("`reinterpret(reshape, $T, a)` where `eltype(a)` is $(eltype(a)) requires that `axes(a, 1)` (got $(axes(a, 1))) be equal to 1:$(sizeof(T) ÷ sizeof(eltype(a))) (from the ratio of element sizes)"))
+        end
+        isbitstype(T) || throwbits(S, T, T)
+        isbitstype(S) || throwbits(S, T, S)
+        if sizeof(S) == sizeof(T)
+            N = ndims(a)
+        elseif sizeof(S) > sizeof(T)
+            rem(sizeof(S), sizeof(T)) == 0 || throwintmult(S, T)
+            N = ndims(a) + 1
+        else
+            rem(sizeof(T), sizeof(S)) == 0 || throwintmult(S, T)
+            N = ndims(a) - 1
+            N > -1 || throwsize0(S, T, "larger")
+            axes(a, 1) == Base.OneTo(sizeof(T) ÷ sizeof(S)) || throwsize1(a, T)
+        end
+        readable = array_subpadding(T, S)
+        writable = array_subpadding(S, T)
+        new{T, N, S, A, true}(a, readable, writable)
     end
 end
 
+ReshapedReinterpretArray{T,N,S,A<:AbstractArray{S}} = ReinterpretArray{T,N,S,A,true}
+NonReshapedReinterpretArray{T,N,S,A<:AbstractArray{S, N}} = ReinterpretArray{T,N,S,A,false}
+
+"""
+    reinterpret(reshape, T, A::AbstractArray{S}) -> B
+
+Change the type-interpretation of `A` while consuming or adding a "channel dimension."
+
+If `sizeof(T) = n*sizeof(S)` for `n>1`, `A`'s first dimension must be
+of size `n` and `B` lacks `A`'s first dimension. Conversely, if `sizeof(S) = n*sizeof(T)` for `n>1`,
+`B` gets a new first dimension of size `n`. The dimensionality is unchanged if `sizeof(T) == sizeof(S)`.
+
+# Examples
+
+```jldoctest
+julia> A = [1 2; 3 4]
+2×2 Matrix{$Int}:
+ 1  2
+ 3  4
+
+julia> reinterpret(reshape, Complex{Int}, A)    # the result is a vector
+2-element reinterpret(reshape, Complex{$Int}, ::Matrix{$Int}):
+ 1 + 3im
+ 2 + 4im
+
+julia> a = [(1,2,3), (4,5,6)]
+2-element Vector{Tuple{$Int, $Int, $Int}}:
+ (1, 2, 3)
+ (4, 5, 6)
+
+julia> reinterpret(reshape, Int, a)             # the result is a matrix
+3×2 reinterpret(reshape, $Int, ::Vector{Tuple{$Int, $Int, $Int}}):
+ 1  4
+ 2  5
+ 3  6
+```
+"""
+reinterpret(::typeof(reshape), T::Type, a::AbstractArray)
+
 reinterpret(::Type{T}, a::ReinterpretArray) where {T} = reinterpret(T, a.parent)
 
 # Definition of StridedArray
 StridedFastContiguousSubArray{T,N,A<:DenseArray} = FastContiguousSubArray{T,N,A}
-StridedReinterpretArray{T,N,A<:Union{DenseArray,StridedFastContiguousSubArray}} = ReinterpretArray{T,N,S,A} where S
+StridedReinterpretArray{T,N,A<:Union{DenseArray,StridedFastContiguousSubArray},IsReshaped} = ReinterpretArray{T,N,S,A,IsReshaped} where S
 StridedReshapedArray{T,N,A<:Union{DenseArray,StridedFastContiguousSubArray,StridedReinterpretArray}} = ReshapedArray{T,N,A}
 StridedSubArray{T,N,A<:Union{DenseArray,StridedReshapedArray,StridedReinterpretArray},
     I<:Tuple{Vararg{Union{RangeIndex, ReshapedUnitRange, AbstractCartesianIndex}}}} = SubArray{T,N,A,I}
@@ -105,31 +175,132 @@ function check_writable(a::ReinterpretArray{T, N, S} where N) where {T,S}
     end
 end
 
-IndexStyle(a::ReinterpretArray) = IndexStyle(a.parent)
+## IndexStyle specializations
+
+# For `reinterpret(reshape, T, a)` where we're adding a channel dimension and with
+# `IndexStyle(a) == IndexLinear()`, it's advantageous to retain pseudo-linear indexing.
+struct IndexSCartesian2{K} <: IndexStyle end   # K = sizeof(S) ÷ sizeof(T), a static-sized 2d cartesian iterator
+
+IndexStyle(::Type{ReinterpretArray{T,N,S,A,false}}) where {T,N,S,A<:AbstractArray{S,N}} = IndexStyle(A)
+function IndexStyle(::Type{ReinterpretArray{T,N,S,A,true}}) where {T,N,S,A<:AbstractArray{S}}
+    if sizeof(T) < sizeof(S)
+        IndexStyle(A) === IndexLinear() && return IndexSCartesian2{sizeof(S) ÷ sizeof(T)}()
+        return IndexCartesian()
+    end
+    return IndexStyle(A)
+end
+IndexStyle(::IndexSCartesian2{K}, ::IndexSCartesian2{K}) where {K} = IndexSCartesian2{K}()
+
+struct SCartesianIndex2{K}   # can't make <:AbstractCartesianIndex without N, and 2 would be a bit misleading
+    i::Int
+    j::Int
+end
+to_index(i::SCartesianIndex2) = i
+
+struct SCartesianIndices2{K,R<:AbstractUnitRange{Int}} <: AbstractMatrix{SCartesianIndex2{K}}
+    indices2::R
+end
+SCartesianIndices2{K}(indices2::AbstractUnitRange{Int}) where {K} = (@assert K::Int > 1; SCartesianIndices2{K,typeof(indices2)}(indices2))
+
+eachindex(::IndexSCartesian2{K}, A::AbstractArray) where {K} = SCartesianIndices2{K}(eachindex(IndexLinear(), parent(A)))
+
+size(iter::SCartesianIndices2{K}) where K = (K, length(iter.indices2))
+axes(iter::SCartesianIndices2{K}) where K = (Base.OneTo(K), iter.indices2)
+
+first(iter::SCartesianIndices2{K}) where {K} = SCartesianIndex2{K}(1, first(iter.indices2))
+last(iter::SCartesianIndices2{K}) where {K}  = SCartesianIndex2{K}(K, last(iter.indices2))
+
+@inline function getindex(iter::SCartesianIndices2{K}, i::Int, j::Int) where {K}
+    @boundscheck checkbounds(iter, i, j)
+    return SCartesianIndex2{K}(i, iter.indices2[j])
+end
+
+function iterate(iter::SCartesianIndices2{K}) where {K}
+    ret = iterate(iter.indices2)
+    ret === nothing && return nothing
+    item2, state2 = ret
+    return SCartesianIndex2{K}(1, item2), (1, item2, state2)
+end
+
+function iterate(iter::SCartesianIndices2{K}, (state1, item2, state2)) where {K}
+    if state1 < K
+        item1 = state1 + 1
+        return SCartesianIndex2{K}(item1, item2), (item1, item2, state2)
+    end
+    ret = iterate(iter.indices2, state2)
+    ret === nothing && return nothing
+    item2, state2 = ret
+    return SCartesianIndex2{K}(1, item2), (1, item2, state2)
+end
+
+SimdLoop.simd_outer_range(iter::SCartesianIndices2) = iter.indices2
+SimdLoop.simd_inner_length(::SCartesianIndices2{K}, ::Any) where K = K
+@inline function SimdLoop.simd_index(::SCartesianIndices2{K}, Ilast::Int, I1::Int) where {K}
+    SCartesianIndex2{K}(I1+1, Ilast)
+end
+
+_maybe_reshape(::IndexSCartesian2, A::ReshapedReinterpretArray, I...) = A
+
+# fallbacks
+function _getindex(::IndexSCartesian2, A::AbstractArray{T,N}, I::Vararg{Int, N}) where {T,N}
+    @_propagate_inbounds_meta
+    getindex(A, I...)
+end
+function _setindex!(::IndexSCartesian2, A::AbstractArray{T,N}, v, I::Vararg{Int, N}) where {T,N}
+    @_propagate_inbounds_meta
+    setindex!(A, v, I...)
+end
+# fallbacks for array types that use "pass-through" indexing (e.g., `IndexStyle(A) = IndexStyle(parent(A))`)
+# but which don't handle SCartesianIndex2
+function _getindex(::IndexSCartesian2, A::AbstractArray{T,N}, ind::SCartesianIndex2) where {T,N}
+    @_propagate_inbounds_meta
+    I = _to_subscript_indices(A, ind.i, ind.j)
+    getindex(A, I...)
+end
+function _setindex!(::IndexSCartesian2, A::AbstractArray{T,N}, v, ind::SCartesianIndex2) where {T,N}
+    @_propagate_inbounds_meta
+    I = _to_subscript_indices(A, ind.i, ind.j)
+    setindex!(A, v, I...)
+end
+
+
+## AbstractArray interface
 
 parent(a::ReinterpretArray) = a.parent
 dataids(a::ReinterpretArray) = dataids(a.parent)
 unaliascopy(a::ReinterpretArray{T}) where {T} = reinterpret(T, unaliascopy(a.parent))
 
-function size(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
+function size(a::NonReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
     psize = size(a.parent)
     size1 = div(psize[1]*sizeof(S), sizeof(T))
     tuple(size1, tail(psize)...)
 end
-size(a::ReinterpretArray{T,0}) where {T} = ()
+function size(a::ReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
+    psize = size(a.parent)
+    sizeof(S) > sizeof(T) && return (div(sizeof(S), sizeof(T)), psize...)
+    sizeof(S) < sizeof(T) && return Base.tail(psize)
+    return psize
+end
+size(a::NonReshapedReinterpretArray{T,0}) where {T} = ()
 
-function axes(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
+function axes(a::NonReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
     paxs = axes(a.parent)
     f, l = first(paxs[1]), length(paxs[1])
     size1 = div(l*sizeof(S), sizeof(T))
     tuple(oftype(paxs[1], f:f+size1-1), tail(paxs)...)
 end
-axes(a::ReinterpretArray{T,0}) where {T} = ()
+function axes(a::ReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
+    paxs = axes(a.parent)
+    sizeof(S) > sizeof(T) && return (Base.OneTo(div(sizeof(S), sizeof(T))), paxs...)
+    sizeof(S) < sizeof(T) && return Base.tail(paxs)
+    return paxs
+end
+axes(a::NonReshapedReinterpretArray{T,0}) where {T} = ()
 
 elsize(::Type{<:ReinterpretArray{T}}) where {T} = sizeof(T)
 unsafe_convert(::Type{Ptr{T}}, a::ReinterpretArray{T,N,S} where N) where {T,S} = Ptr{T}(unsafe_convert(Ptr{S},a.parent))
 
-@inline @propagate_inbounds getindex(a::ReinterpretArray{T,0}) where {T} = reinterpret(T, a.parent[])
+@inline @propagate_inbounds getindex(a::NonReshapedReinterpretArray{T,0}) where {T} = reinterpret(T, a.parent[])
 @inline @propagate_inbounds getindex(a::ReinterpretArray) = a[1]
 
 @inline @propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, inds::Vararg{Int, N}) where {T,N,S}
@@ -148,9 +319,22 @@ end
     _getindex_ra(a, inds[1], tail(inds))
 end
 
+@inline @propagate_inbounds function getindex(a::ReshapedReinterpretArray{T,N,S}, ind::SCartesianIndex2) where {T,N,S}
+    check_readable(a)
+    n = sizeof(S) ÷ sizeof(T)
+    t = Ref{NTuple{n,T}}()
+    s = Ref{S}(a.parent[ind.j])
+    GC.@preserve t s begin
+        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
+        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
+        _memcpy!(tptr, sptr, sizeof(S))
+    end
+    return t[][ind.i]
+end
+
 @inline _memcpy!(dst, src, n) = ccall(:memcpy, Cvoid, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), dst, src, n)
 
-@inline @propagate_inbounds function _getindex_ra(a::ReinterpretArray{T,N,S}, i1::Int, tailinds::TT) where {T,N,S,TT}
+@inline @propagate_inbounds function _getindex_ra(a::NonReshapedReinterpretArray{T,N,S}, i1::Int, tailinds::TT) where {T,N,S,TT}
     # Make sure to match the scalar reinterpret if that is applicable
     if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
         return reinterpret(T, a.parent[i1, tailinds...])
@@ -194,8 +378,55 @@ end
     end
 end
 
+@inline @propagate_inbounds function _getindex_ra(a::ReshapedReinterpretArray{T,N,S}, i1::Int, tailinds::TT) where {T,N,S,TT}
+    # Make sure to match the scalar reinterpret if that is applicable
+    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
+        return reinterpret(T, a.parent[i1, tailinds...])
+    end
+    @boundscheck checkbounds(a, i1, tailinds...)
+    if sizeof(T) >= sizeof(S)
+        t = Ref{T}()
+        s = Ref{S}()
+        GC.@preserve t s begin
+            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
+            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
+            if sizeof(T) > sizeof(S)
+                # Extra dimension in the parent array
+                n = sizeof(T) ÷ sizeof(S)
+                if isempty(tailinds) && IndexStyle(a.parent) === IndexLinear()
+                    offset = n * (i1 - firstindex(a))
+                    for i = 1:n
+                        s[] = a.parent[i + offset]
+                        _memcpy!(tptr + (i-1)*sizeof(S), sptr, sizeof(S))
+                    end
+                else
+                    for i = 1:n
+                        s[] = a.parent[i, i1, tailinds...]
+                        _memcpy!(tptr + (i-1)*sizeof(S), sptr, sizeof(S))
+                    end
+                end
+            else
+                # No extra dimension
+                s[] = a.parent[i1, tailinds...]
+                _memcpy!(tptr, sptr, sizeof(S))
+            end
+        end
+        return t[]
+    end
+    # S is bigger than T and contains an integer number of them
+    n = sizeof(S) ÷ sizeof(T)
+    t = Ref{NTuple{n,T}}()
+    s = Ref{S}()
+    GC.@preserve t s begin
+        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
+        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
+        s[] = a.parent[tailinds...]
+        _memcpy!(tptr, sptr, sizeof(S))
+    end
+    return t[][i1]
+end
 
-@inline @propagate_inbounds setindex!(a::ReinterpretArray{T,0,S} where T, v) where {S} = (a.parent[] = reinterpret(S, v))
+@inline @propagate_inbounds setindex!(a::NonReshapedReinterpretArray{T,0,S} where T, v) where {S} = (a.parent[] = reinterpret(S, v))
 @inline @propagate_inbounds setindex!(a::ReinterpretArray, v) = (a[1] = v)
 
 @inline @propagate_inbounds function setindex!(a::ReinterpretArray{T,N,S}, v, inds::Vararg{Int, N}) where {T,N,S}
@@ -212,7 +443,21 @@ end
     _setindex_ra!(a, v, inds[1], tail(inds))
 end
 
-@inline @propagate_inbounds function _setindex_ra!(a::ReinterpretArray{T,N,S}, v, i1::Int, tailinds::TT) where {T,N,S,TT}
+@inline @propagate_inbounds function setindex!(a::ReshapedReinterpretArray{T,N,S}, v, ind::SCartesianIndex2) where {T,N,S}
+    check_writable(a)
+    v = convert(T, v)::T
+    t = Ref{T}(v)
+    s = Ref{S}(a.parent[ind.j])
+    GC.@preserve t s begin
+        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
+        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
+        _memcpy!(sptr + (ind.i-1)*sizeof(T), tptr, sizeof(T))
+    end
+    a.parent[ind.j] = s[]
+    return a
+end
+
+@inline @propagate_inbounds function _setindex_ra!(a::NonReshapedReinterpretArray{T,N,S}, v, i1::Int, tailinds::TT) where {T,N,S,TT}
     v = convert(T, v)::T
     # Make sure to match the scalar reinterpret if that is applicable
     if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
@@ -273,6 +518,49 @@ end
     return a
 end
 
+@inline @propagate_inbounds function _setindex_ra!(a::ReshapedReinterpretArray{T,N,S}, v, i1::Int, tailinds::TT) where {T,N,S,TT}
+    v = convert(T, v)::T
+    # Make sure to match the scalar reinterpret if that is applicable
+    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
+        return setindex!(a.parent, reinterpret(S, v), i1, tailinds...)
+    end
+    @boundscheck checkbounds(a, i1, tailinds...)
+    t = Ref{T}(v)
+    s = Ref{S}()
+    GC.@preserve t s begin
+        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
+        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
+        if sizeof(T) >= sizeof(S)
+            if sizeof(T) > sizeof(S)
+                # Extra dimension in the parent array
+                n = sizeof(T) ÷ sizeof(S)
+                if isempty(tailinds) && IndexStyle(a.parent) === IndexLinear()
+                    offset = n * (i1 - firstindex(a))
+                    for i = 1:n
+                        _memcpy!(sptr, tptr + (i-1)*sizeof(S), sizeof(S))
+                        a.parent[i + offset] = s[]
+                    end
+                else
+                    for i = 1:n
+                        _memcpy!(sptr, tptr + (i-1)*sizeof(S), sizeof(S))
+                        a.parent[i, i1, tailinds...] = s[]
+                    end
+                end
+            else
+                # No extra dimension
+                _memcpy!(sptr, tptr, sizeof(S))
+                a.parent[i1, tailinds...] = s[]
+            end
+        else
+            # S is bigger than T and contains an integer number of them
+            s[] = a.parent[tailinds...]
+            _memcpy!(sptr + (i1-1)*sizeof(T), tptr, sizeof(T))
+            a.parent[tailinds...] = s[]
+        end
+    end
+    return a
+end
+
 # Padding
 struct Padding
     offset::Int
@@ -366,3 +654,46 @@ using .Iterators: Stateful
     end
     return true
 end
+
+# Reductions with IndexSCartesian2
+
+function _mapreduce(f::F, op::OP, style::IndexSCartesian2{K}, A::AbstractArrayOrBroadcasted) where {F,OP,K}
+    inds = eachindex(style, A)
+    n = size(inds)[2]
+    if n == 0
+        return mapreduce_empty_iter(f, op, A, IteratorEltype(A))
+    else
+        return mapreduce_impl(f, op, A, first(inds), last(inds))
+    end
+end
+
+@noinline function mapreduce_impl(f::F, op::OP, A::AbstractArrayOrBroadcasted,
+                                  ifirst::SCI, ilast::SCI, blksize::Int) where {F,OP,SCI<:SCartesianIndex2{K}} where K
+    if ifirst.j + blksize > ilast.j
+        # sequential portion
+        @inbounds a1 = A[ifirst]
+        @inbounds a2 = A[SCI(2,ifirst.j)]
+        v = op(f(a1), f(a2))
+        @simd for i = ifirst.i + 2 : K
+            @inbounds ai = A[SCI(i,ifirst.j)]
+            v = op(v, f(ai))
+        end
+        # Remaining columns
+        for j = ifirst.j+1 : ilast.j
+            @simd for i = 1:K
+                @inbounds ai = A[SCI(i,j)]
+                v = op(v, f(ai))
+            end
+        end
+        return v
+    else
+        # pairwise portion
+        jmid = (ifirst.j + ilast.j) >> 1
+        v1 = mapreduce_impl(f, op, A, ifirst, SCI(K,jmid), blksize)
+        v2 = mapreduce_impl(f, op, A, SCI(1,jmid+1), ilast, blksize)
+        return op(v1, v2)
+    end
+end
+
+mapreduce_impl(f::F, op::OP, A::AbstractArrayOrBroadcasted, ifirst::SCartesianIndex2, ilast::SCartesianIndex2) where {F,OP} =
+    mapreduce_impl(f, op, A, ifirst, ilast, pairwise_blocksize(f, op))
diff --git a/base/show.jl b/base/show.jl
index 1984a6349edea..2abb094c12acc 100644
--- a/base/show.jl
+++ b/base/show.jl
@@ -2562,12 +2562,18 @@ function showarg(io::IO, r::ReshapedArray, toplevel)
     toplevel && print(io, " with eltype ", eltype(r))
 end
 
-function showarg(io::IO, r::ReinterpretArray{T}, toplevel) where {T}
+function showarg(io::IO, r::NonReshapedReinterpretArray{T}, toplevel) where {T}
     print(io, "reinterpret(", T, ", ")
     showarg(io, parent(r), false)
     print(io, ')')
 end
 
+function showarg(io::IO, r::ReshapedReinterpretArray{T}, toplevel) where {T}
+    print(io, "reinterpret(reshape, ", T, ", ")
+    showarg(io, parent(r), false)
+    print(io, ')')
+end
+
 # printing iterators from Base.Iterators
 
 function show(io::IO, e::Iterators.Enumerate)
diff --git a/test/reinterpretarray.jl b/test/reinterpretarray.jl
index c6990066c1423..3791ad78db8d7 100644
--- a/test/reinterpretarray.jl
+++ b/test/reinterpretarray.jl
@@ -3,29 +3,107 @@
 using Test
 isdefined(Main, :OffsetArrays) || @eval Main include("testhelpers/OffsetArrays.jl")
 using .Main.OffsetArrays
+isdefined(Main, :TSlow) || @eval Main include("testhelpers/arrayindexingtypes.jl")
+using .Main: TSlow, WrapperArray
 
 A = Int64[1, 2, 3, 4]
+As = TSlow(A)
+Ars = Int64[1 3; 2 4]
+Arss = TSlow(Ars)
 B = Complex{Int64}[5+6im, 7+8im, 9+10im]
+Bs = TSlow(B)
+Av = [Int32[1,2], Int32[3,4]]
+
+for Ar in (Ars, Arss)
+    @test @inferred(ndims(reinterpret(reshape, Complex{Int64}, Ar))) == 1
+    @test @inferred(axes(reinterpret(reshape, Complex{Int64}, Ar))) === (Base.OneTo(2),)
+    @test @inferred(size(reinterpret(reshape, Complex{Int64}, Ar))) == (2,)
+end
+for _B in (B, Bs)
+    @test @inferred(ndims(reinterpret(reshape, Int64, _B))) == 2
+    @test @inferred(axes(reinterpret(reshape, Int64, _B))) === (Base.OneTo(2), Base.OneTo(3))
+    @test @inferred(size(reinterpret(reshape, Int64, _B))) == (2, 3)
+    @test @inferred(ndims(reinterpret(reshape, Int128, _B))) == 1
+    @test @inferred(axes(reinterpret(reshape, Int128, _B))) === (Base.OneTo(3),)
+    @test @inferred(size(reinterpret(reshape, Int128, _B))) == (3,)
+end
+
+@test_throws ArgumentError("cannot reinterpret `Int64` as `Vector{Int64}`, type `Vector{Int64}` is not a bits type") reinterpret(Vector{Int64}, A)
+@test_throws ArgumentError("cannot reinterpret `Vector{Int32}` as `Int32`, type `Vector{Int32}` is not a bits type") reinterpret(Int32, Av)
+@test_throws ArgumentError("cannot reinterpret a zero-dimensional `Int64` array to `Int32` which is of a different size") reinterpret(Int32, reshape([Int64(0)]))
+@test_throws ArgumentError("cannot reinterpret a zero-dimensional `Int32` array to `Int64` which is of a different size") reinterpret(Int64, reshape([Int32(0)]))
+@test_throws ArgumentError("""cannot reinterpret an `$Int` array to `Tuple{$Int, $Int}` whose first dimension has size `5`.
+                              The resulting array would have non-integral first dimension.
+                              """) reinterpret(Tuple{Int,Int}, [1,2,3,4,5])
+
+@test_throws ArgumentError("`reinterpret(reshape, Complex{Int64}, a)` where `eltype(a)` is Int64 requires that `axes(a, 1)` (got Base.OneTo(4)) be equal to 1:2 (from the ratio of element sizes)") reinterpret(reshape, Complex{Int64}, A)
+@test_throws ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got 24) and `sizeof(eltype(a))` (got 16) be an integer multiple of the other") reinterpret(reshape, NTuple{3, Int64}, B)
+@test_throws ArgumentError("cannot reinterpret `Int64` as `Vector{Int64}`, type `Vector{Int64}` is not a bits type") reinterpret(reshape, Vector{Int64}, Ars)
+@test_throws ArgumentError("cannot reinterpret a zero-dimensional `UInt8` array to `UInt16` which is of a larger size") reinterpret(reshape, UInt16, reshape([0x01]))
+
 # getindex
-@test reinterpret(Complex{Int64}, A) == [1 + 2im, 3 + 4im]
-@test reinterpret(Float64, A) == reinterpret.(Float64, A)
+for _A in (A, As)
+    @test reinterpret(Complex{Int64}, _A) == [1 + 2im, 3 + 4im]
+    @test reinterpret(Float64, _A) == reinterpret.(Float64, A)
+    @test reinterpret(reshape, Float64, _A) == reinterpret.(Float64, A)
+end
+for Ar in (Ars, Arss)
+    @test reinterpret(reshape, Complex{Int64}, Ar) == [1 + 2im, 3 + 4im]
+    @test reinterpret(reshape, Float64, Ar) == reinterpret.(Float64, Ars)
+end
 
-@test reinterpret(NTuple{3, Int64}, B) == [(5,6,7),(8,9,10)]
+for _B in (B, Bs)
+    @test reinterpret(NTuple{3, Int64}, _B) == [(5,6,7),(8,9,10)]
+    @test reinterpret(reshape, Int64, _B) == [5 7 9; 6 8 10]
+end
 
 # setindex
-let Ac = copy(A), Bc = copy(B)
-    reinterpret(Complex{Int64}, Ac)[2] = -1 - 2im
-    @test Ac == [1, 2, -1, -2]
-    reinterpret(NTuple{3, Int64}, Bc)[2] = (4,5,6)
-    @test Bc == Complex{Int64}[5+6im, 7+4im, 5+6im]
-    reinterpret(NTuple{3, Int64}, Bc)[1] = (1,2,3)
-    @test Bc == Complex{Int64}[1+2im, 3+4im, 5+6im]
+for (_A, Ar, _B) in ((A, Ars, B), (As, Arss, Bs))
+    let Ac = copy(_A), Arsc = copy(Ar), Bc = copy(_B)
+        reinterpret(Complex{Int64}, Ac)[2] = -1 - 2im
+        @test Ac == [1, 2, -1, -2]
+        reinterpret(Complex{Int64}, Arsc)[2] = -1 - 2im
+        @test Arsc == [1 -1; 2 -2]
+        reinterpret(NTuple{3, Int64}, Bc)[2] = (4,5,6)
+        @test Bc == Complex{Int64}[5+6im, 7+4im, 5+6im]
+        reinterpret(NTuple{3, Int64}, Bc)[1] = (1,2,3)
+        @test Bc == Complex{Int64}[1+2im, 3+4im, 5+6im]
+        Bc = copy(_B)
+        Brrs = reinterpret(reshape, Int64, Bc)
+        Brrs[2, 3] = -5
+        @test Bc == Complex{Int64}[5+6im, 7+8im, 9-5im]
+        Brrs[last(eachindex(Brrs))] = 22
+        @test Bc == Complex{Int64}[5+6im, 7+8im, 9+22im]
 
-    A1 = reinterpret(Float64, A)
-    A2 = reinterpret(ComplexF64, A)
-    A1[1] = 1.0
-    @test real(A2[1]) == 1.0
+        A1 = reinterpret(Float64, _A)
+        A2 = reinterpret(ComplexF64, _A)
+        A1[1] = 1.0
+        @test real(A2[1]) == 1.0
+        A1 = reinterpret(reshape, Float64, _A)
+        A1[1] = 2.5
+        @test reinterpret(Float64, _A[1]) == 2.5
+        A1rs = reinterpret(Float64, Ar)
+        A2rs = reinterpret(ComplexF64, Ar)
+        A1rs[1, 1] = 1.0
+        @test real(A2rs[1]) == 1.0
+        A1rs = reinterpret(reshape, Float64, Ar)
+        A2rs = reinterpret(reshape, ComplexF64, Ar)
+        A1rs[1, 1] = 2.5
+        @test real(A2rs[1]) == 2.5
+    end
 end
+A3 = collect(reshape(1:18, 2, 3, 3))
+A3r = reinterpret(reshape, Complex{Int}, A3)
+@test A3r[4] === A3r[1,2] === A3r[CartesianIndex(1, 2)] === 7+8im
+A3r[2,3] = -8-15im
+@test A3[1,2,3] == -8
+@test A3[2,2,3] == -15
+A3r[4] = 100+200im
+@test A3[1,1,2] == 100
+@test A3[2,1,2] == 200
+A3r[CartesianIndex(1,2)] = 300+400im
+@test A3[1,1,2] == 300
+@test A3[2,1,2] == 400
 
 # same-size reinterpret where one of the types is non-primitive
 let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
@@ -33,6 +111,25 @@ let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
     reinterpret(Float32, a)[1] = 2.0
     @test reinterpret(Float32, a)[1] == 2.0
 end
+let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
+    @test reinterpret(reshape, Float32, a)[1] == reinterpret(Float32, 0x04030201)
+    reinterpret(reshape, Float32, a)[1] = 2.0
+    @test reinterpret(reshape, Float32, a)[1] == 2.0
+end
+
+# Pass-through indexing
+B = Complex{Int64}[5+6im, 7+8im, 9+10im]
+Br = reinterpret(reshape, Int64, B)
+W = WrapperArray(Br)
+for (b, w) in zip(5:10, W)
+    @test b == w
+end
+for (i, j) in zip(eachindex(W), 11:16)
+    W[i] = j
+end
+@test B[1] === Complex{Int64}(11+12im)
+@test B[2] === Complex{Int64}(13+14im)
+@test B[3] === Complex{Int64}(15+16im)
 
 # ensure that reinterpret arrays aren't erroneously classified as strided
 let A = reshape(1:20, 5, 4)
@@ -166,7 +263,9 @@ let a = [0.1 0.2; 0.3 0.4], at = reshape([(i,i+1) for i = 1:2:8], 2, 2)
     @test r[0,3] === reinterpret(Int64, v[0,3])
     @test r[1,3] === reinterpret(Int64, v[1,3])
     @test_throws ArgumentError("cannot reinterpret a `Float64` array to `UInt32` when the first axis is Base.IdentityUnitRange(0:1). Try reshaping first.") reinterpret(UInt32, v)
+    @test_throws ArgumentError("`reinterpret(reshape, Tuple{Float64, Float64}, a)` where `eltype(a)` is Float64 requires that `axes(a, 1)` (got Base.IdentityUnitRange(0:1)) be equal to 1:2 (from the ratio of element sizes)") reinterpret(reshape, Tuple{Float64,Float64}, v)
     v = OffsetArray(a, (0, 1))
+    @test axes(reinterpret(reshape, Tuple{Float64,Float64}, v)) === (Base.IdentityUnitRange(2:3),)
     r = reinterpret(UInt32, v)
     axsv = axes(v)
     @test axes(r) === (oftype(axsv[1], 1:4), axsv[2])
@@ -206,8 +305,22 @@ end
 # Test 0-dimensional Arrays
 A = zeros(UInt32)
 B = reinterpret(Int32,A)
-@test size(B) == ()
-@test axes(B) == ()
+Brs = reinterpret(reshape,Int32,A)
+@test size(B) == size(Brs) == ()
+@test axes(B) == axes(Brs) == ()
 B[] = Int32(5)
 @test B[] === Int32(5)
+@test Brs[] === Int32(5)
 @test A[] === UInt32(5)
+
+# reductions
+a = [(1,2,3), (4,5,6)]
+ars = reinterpret(reshape, Int, a)
+@test sum(ars) == 21
+@test sum(ars; dims=1) == [6 15]
+@test sum(ars; dims=2) == reshape([5,7,9], (3, 1))
+@test sum(ars; dims=(1,2)) == reshape([21], (1, 1))
+# also test large sizes for the pairwise algorithm
+a = [(k,k+1,k+2) for k = 1:3:4000]
+ars = reinterpret(reshape, Int, a)
+@test sum(ars) == 8010003
diff --git a/test/testhelpers/arrayindexingtypes.jl b/test/testhelpers/arrayindexingtypes.jl
index b26ffc9ec5f24..408559871ae19 100644
--- a/test/testhelpers/arrayindexingtypes.jl
+++ b/test/testhelpers/arrayindexingtypes.jl
@@ -52,3 +52,15 @@ Base.similar(A::TSlow, ::Type{T}, dims::Dims) where {T} = TSlow(T, dims)
 Base.IndexStyle(::Type{A}) where {A<:TSlow} = IndexCartesian()
 Base.getindex(A::TSlow{T,N}, i::Vararg{Int,N}) where {T,N} = get(A.data, i, zero(T))
 Base.setindex!(A::TSlow{T,N}, v, i::Vararg{Int,N}) where {T,N} = (A.data[i] = v)
+
+# An array type that just passes through to the parent
+struct WrapperArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
+    parent::A
+end
+Base.IndexStyle(::Type{WrapperArray{T,N,A}}) where {T,N,A<:AbstractArray{T,N}} = IndexStyle(A)
+Base.parent(A::WrapperArray) = A.parent
+Base.size(A::WrapperArray) = size(A.parent)
+Base.axes(A::WrapperArray) = axes(A.parent)
+Base.getindex(A::WrapperArray, i::Int...) = A.parent[i...]
+Base.setindex!(A::WrapperArray, v, i::Int...) = A.parent[i...] = v
+Base.similar(A::WrapperArray, ::Type{T}, dims::Dims) where T = similar(A.parent, T, dims)