diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index fc873f441385f..5fe12ec075051 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -174,12 +174,6 @@ similar{T}(a::AbstractArray{T}, dims::Dims)   = similar(a, T, dims)
 similar{T}(a::AbstractArray{T}, dims::Int...) = similar(a, T, dims)
 similar   (a::AbstractArray, T, dims::Int...) = similar(a, T, dims)
 
-function reshape(a::AbstractArray, dims::Dims)
-    if prod(dims) != length(a)
-        throw(ArgumentError("dimensions must be consistent with array size"))
-    end
-    copy!(similar(a, dims), a)
-end
 reshape(a::AbstractArray, dims::Int...) = reshape(a, dims)
 
 vec(a::AbstractArray) = reshape(a,length(a))
diff --git a/base/int.jl b/base/int.jl
index a79ee16a563c1..399574b6cfd1e 100644
--- a/base/int.jl
+++ b/base/int.jl
@@ -514,3 +514,140 @@ checked_mul(x::Int128, y::Int128) = x * y
 checked_add(x::UInt128, y::UInt128) = x + y
 checked_sub(x::UInt128, y::UInt128) = x - y
 checked_mul(x::UInt128, y::UInt128) = x * y
+
+
+## A faster alternative to div if the dividend is reused many times
+
+unsigned_type(::Int8) = UInt8
+unsigned_type(::Int16) = UInt16
+unsigned_type(::Int32) = UInt32
+unsigned_type(::Int64) = UInt64
+unsigned_type(::Int128) = UInt128
+
+abstract FastDivInteger{T}
+
+immutable SignedFastDivInteger{T<:Signed} <: FastDivInteger{T}
+    divisor::T
+    multiplier::T
+    addmul::Int8
+    shift::UInt8
+
+    function SignedFastDivInteger(d::T)
+        ut = unsigned_type(d)
+        signedmin = reinterpret(ut, typemin(d))
+
+        ad::ut = abs(d)
+        ad <= 1 && error("cannot compute magic for d == $d")
+        t::ut = signedmin + signbit(d)
+        anc::ut = t - 1 - rem(t, ad)   # absolute value of nc
+        p = sizeof(d)*8 - 1            # initialize p
+        q1::ut, r1::ut = divrem(signedmin, anc)
+        q2::ut, r2::ut = divrem(signedmin, ad)
+        while true
+            p += 1
+            q1 *= 2                    # update q1 = 2p/abs(nc)
+            r1 *= 2                    # update r1 = rem(2p/abs(nc))
+            if r1 >= anc               # must be unsigned comparison
+                q1 += 1
+                r1 -= anc
+            end
+            q2 *= 2                    # update q2 = 2p/abs(d)
+            r2 *= 2                    # update r2 = rem(2p/abs(d))
+            if r2 >= ad                # must be unsigned comparison
+                q2 += 1
+                r2 -= ad
+            end
+            delta::ut = ad - r2
+            (q1 < delta || (q1 == delta && r1 == 0)) || break
+        end
+
+        m = flipsign((q2 + 1) % T, d)    # resulting magic number
+        s = p - sizeof(d)*8                   # resulting shift
+        new(d, m, d > 0 && m < 0 ? Int8(1) : d < 0 && m > 0 ? Int8(-1) : Int8(0), UInt8(s))
+    end
+end
+SignedFastDivInteger(x::Signed) = SignedFastDivInteger{typeof(x)}(x)
+
+immutable UnsignedFastDivInteger{T<:Unsigned} <: FastDivInteger{T}
+    divisor::T
+    multiplier::T
+    add::Bool
+    shift::UInt8
+
+    function UnsignedFastDivInteger(d::T)
+        (d == 0 || d == 1) && error("cannot compute magic for d == $d")
+        u2 = convert(T, 2)
+        add = false
+        signedmin::typeof(d) = one(d) << (sizeof(d)*8-1)
+        signedmax::typeof(d) = signedmin - 1
+        allones = (zero(d) - 1) % T
+
+        nc::typeof(d) = allones - rem(convert(T, allones - d), d)
+        p = 8*sizeof(d) - 1                    # initialize p
+        q1::typeof(d), r1::typeof(d) = divrem(signedmin, nc)
+        q2::typeof(d), r2::typeof(d) = divrem(signedmax, d)
+        while true
+            p += 1
+            if r1 >= convert(T, nc - r1)
+                q1 = q1 + q1 + T(1)     # update q1
+                r1 = r1 + r1 - nc       # update r1
+            else
+                q1 = q1 + q1            # update q1
+                r1 = r1 + r1            # update r1
+            end
+            if convert(T, r2 + T(1)) >= convert(T, d - r2)
+                add |= q2 >= signedmax
+                q2 = q2 + q2 + 1        # update q2
+                r2 = r2 + r2 + T(1) - d # update r2
+            else
+                add |= q2 >= signedmin
+                q2 = q2 + q2            # update q2
+                r2 = r2 + r2 + T(1)     # update r2
+            end
+            delta::typeof(d) = d - 1 - r2
+            (p < sizeof(d)*16 && (q1 < delta || (q1 == delta && r1 == 0))) || break
+        end
+        m = q2 + 1                   # resulting magic number
+        s = p - sizeof(d)*8 - add    # resulting shift
+        new(d, m % T, add, s % UInt8)
+    end
+end
+UnsignedFastDivInteger(x::Unsigned) = UnsignedFastDivInteger{typeof(x)}(x)
+
+# Special type to handle div by 1
+immutable FastDivInteger1{T} <: FastDivInteger{T} end
+
+(*)(a::FastDivInteger, b::FastDivInteger) = a.divisor*b.divisor
+(*)(a::Number, b::FastDivInteger) = a*b.divisor
+(*)(a::FastDivInteger, b::Number) = a.divisor*b
+
+# div{T}(a::Integer, b::SignedFastDivInteger{T})   = div(convert(T, a), b)
+# div{T}(a::Integer, b::UnsignedFastDivInteger{T}) = div(convert(T, a), b)
+# rem{T}(a::Integer, b::FastDivInteger{T})         = rem(convert(T, a), b)
+# divrem{T}(a::Integer, b::FastDivInteger{T})      = divrem(convert(T, a), b)
+
+function div{T}(a::T, b::SignedFastDivInteger{T})
+    x = ((widen(a)*b.multiplier) >>> sizeof(a)*8) % T
+    x += (a*b.addmul) % T
+    (signbit(x) + (x >> b.shift)) % T
+end
+function div{T}(a::T, b::UnsignedFastDivInteger{T})
+    x = ((widen(a)*b.multiplier) >>> sizeof(a)*8) % T
+    x = ifelse(b.add, convert(T, convert(T, (convert(T, a - x) >>> 1)) + x), x)
+    x >>> b.shift
+end
+div{T}(a, ::FastDivInteger1{T}) = convert(T, a)
+
+rem{T}(a::T, b::FastDivInteger{T}) =
+    a - div(a, b)*b.divisor
+rem{T}(a, ::FastDivInteger1{T}) = zero(T)
+
+function divrem{T}(a::T, b::FastDivInteger{T})
+    d = div(a, b)
+    (d, a - d*b.divisor)
+end
+divrem{T}(a, ::FastDivInteger1{T}) = (convert(T, a), zero(T))
+
+# Type unstable!
+call(::Type{FastDivInteger}, x::Signed) = x == 1 ? FastDivInteger1{typeof(x)}() : SignedFastDivInteger(x)
+call(::Type{FastDivInteger}, x::Unsigned) = x == 1 ? FastDivInteger1{typeof(x)}() : UnsignedFastDivInteger(x)
diff --git a/base/reshapedarray.jl b/base/reshapedarray.jl
new file mode 100644
index 0000000000000..2a0288d8aa24d
--- /dev/null
+++ b/base/reshapedarray.jl
@@ -0,0 +1,159 @@
+module Reshaped
+
+import Base: getindex, ind2sub, linearindexing, reshape, similar, size
+# just using, not overloading:
+import Base: LinearFast, LinearSlow, FastDivInteger, FastDivInteger1, tail
+
+# Remapping a block of In dimensions of the parent array to Out dimensions of the view
+immutable IndexMD{In,Out,MI<:(FastDivInteger...)}
+    invstride_parent::MI
+    dims_view::NTuple{Out,Int}
+end
+IndexMD(mi, dims) = IndexMD{length(mi),length(dims),typeof(mi)}(mi, dims)
+
+typealias ReshapeIndex Union(Colon, IndexMD)
+
+immutable ReshapedArray{T,N,P<:AbstractArray,I<:(ReshapeIndex...)} <: AbstractArray{T,N}
+    parent::P
+    indexes::I
+    dims::NTuple{N,Int}
+end
+
+function ReshapedArray(parent::AbstractArray, indexes::(ReshapeIndex...), dims)
+    ReshapedArray{eltype(parent),length(dims),typeof(parent),typeof(indexes)}(parent, indexes, dims)
+end
+
+# Since ranges are immutable and are frequently used to build arrays, treat them as a special case.
+function reshape(a::Range, dims::Dims)
+    if prod(dims) != length(a)
+        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))"))
+    end
+    A = Array(eltype(a), dims)
+    k = 0
+    for item in a
+        A[k+=1] = item
+    end
+    A
+end
+
+reshape(parent::AbstractArray, dims::Dims) = reshape((parent, linearindexing(parent)), dims)
+
+function reshape(p::(AbstractArray,LinearSlow), dims::Dims)
+    parent = p[1]
+    # Split on dimensions where the strides line up
+    stridep = dimstrides(size(parent))
+    stridev = dimstrides(dims)
+    stridep[end] == stridev[end] || throw(DimensionMismatch("Must have the same number of elements"))
+    indexes = Any[]
+    iplast = ivlast = 1
+    ip = iv = 2
+    while ip <= length(stridep) || iv <= length(stridev)
+        if stridep[ip] == stridev[iv]
+            if ip-iplast <= 1 && iv-ivlast == 1
+                push!(indexes, Colon())
+            else
+                mi = map(FastDivInteger, size(parent)[iplast:ip-1])
+                imd = IndexMD(mi, dims[ivlast:iv-1])
+                push!(indexes, imd)
+            end
+            iplast = ip
+            ivlast = iv
+            ip += (ip < length(stridep) || iv == length(stridev))  # for consuming trailing 1s in dims
+            iv += 1
+        elseif stridep[ip] < stridev[iv]
+            ip += 1
+        else
+            iv += 1
+        end
+    end
+    ReshapedArray(parent, tuple(indexes...), dims)
+end
+
+function reshape(p::(AbstractArray,LinearFast), dims::Dims)
+    parent = p[1]
+    prod(dims) == length(parent) || throw(DimensionMismatch("Must have the same number of elements"))
+    ReshapedArray(parent, (IndexMD((FastDivInteger1{Int}(),), dims),), dims)
+end
+
+reshape(parent::ReshapedArray, dims::Dims) = reshape(parent.parent, dims)
+
+size(A::ReshapedArray) = A.dims
+similar{T}(A::ReshapedArray, ::Type{T}, dims::Dims) = Array(T, dims)
+linearindexing(A::ReshapedArray) = linearindexing(A.parent)
+
+size(index::IndexMD) = index.dims_view
+
+@inline function getindex(indx::IndexMD{1}, indexes::Int...)
+    sub2ind(indx.dims_view, indexes...)
+end
+@inline function getindex(indx::IndexMD, indexes::Int...)
+    ind2sub(indx.invstride_parent, sub2ind(indx.dims_view, indexes...))
+end
+
+consumes{In,Out,MI}(::Type{IndexMD{In,Out,MI}}) = Out
+consumes(::Type{Colon}) = 1
+produces{In,Out,MI}(::Type{IndexMD{In,Out,MI}}) = In
+produces(::Type{Colon}) = 1
+
+getindex(A::ReshapedArray) = A.parent[1]
+getindex(A::ReshapedArray, indx::Real) = A.parent[indx]
+
+stagedfunction getindex{T,N,P,I}(A::ReshapedArray{T,N,P,I}, indexes::Real...)
+    length(indexes) == N || throw(DimensionMismatch("Must index with all $N indexes"))
+    c = map(consumes, I)
+    breaks = [0;cumsum([c...])]
+    argbreaks = Any[]
+    for i = 1:length(c)
+        ab = Expr[]
+        for j = breaks[i]+1:breaks[i+1]
+            push!(ab, :(indexes[$j]))
+        end
+        push!(argbreaks, ab)
+    end
+    argsin = Expr[:(getindex(A.indexes[$i], $(argbreaks[i]...))) for i = 1:length(c)]
+    np = map(produces, I)
+    npc = [0;cumsum([np...])]
+    argsout = Expr[]
+    if length(argsin) > 1
+        for i = 1:npc[end]
+            j = findlast(npc .< i)
+            di = i - npc[j]
+            push!(argsout, :(tindex[$j][$di]))
+        end
+    else
+        for i = 1:npc[end]
+            push!(argsout, :(tindex[$i]))
+        end
+    end
+    meta = Expr(:meta, :inline)
+    ex = length(argsin) == 1 ?
+        quote
+            $meta
+            tindex = $(argsin...)
+            getindex(A.parent, $(argsout...))
+        end :
+        quote
+            $meta
+            tindex = tuple($(argsin...))
+            getindex(A.parent, $(argsout...))
+        end
+    ex
+end
+
+# Like strides, but operates on a Dims tuple, and returns one extra element (the total size)
+dimstrides(::()) = ()
+dimstrides(s::Dims) = dimstrides((1,), s)
+dimstrides(t::Tuple, ::()) = t
+@inline dimstrides(t::Tuple, sz::Dims) = dimstrides(tuple(t..., t[end]*sz[1]), tail(sz))
+
+ind2sub(dims::(FastDivInteger,), ind::Integer) = ind
+@inline ind2sub(dims::(FastDivInteger,FastDivInteger), ind::Integer) = begin
+    dv, rm = divrem(ind-1,dims[1])
+    rm+1, dv+1
+end
+@inline ind2sub(dims::(FastDivInteger,FastDivInteger,FastDivInteger...), ind::Integer) = begin
+    dv, rm = divrem(ind-1,dims[1])
+    tuple(rm+1, ind2sub(tail(dims),dv+1)...)
+end
+
+end
diff --git a/base/sysimg.jl b/base/sysimg.jl
index c3a90ff207fb6..229ee45ac2e3a 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -63,6 +63,7 @@ include("array.jl")
 include("subarray2.jl")
 include("functors.jl")
 include("bitarray.jl")
+include("reshapedarray.jl")
 include("intset.jl")
 include("dict.jl")
 include("set.jl")
diff --git a/test/choosetests.jl b/test/choosetests.jl
index c7d757b38387f..e33e1703add7b 100644
--- a/test/choosetests.jl
+++ b/test/choosetests.jl
@@ -15,7 +15,7 @@ function choosetests(choices = [])
     testnames = [
         "linalg", "core", "keywordargs", "numbers", "strings",
         "dates", "dict", "hashing", "remote", "iobuffer", "staged",
-        "arrayops", "subarray", "reduce", "reducedim", "random",
+        "arrayops", "subarray", "reshapedarray", "reduce", "reducedim", "random",
         "intfuncs", "simdloop", "blas", "fft", "dsp", "sparse",
         "bitarray", "copy", "math", "fastmath", "functional",
         "operators", "path", "ccall",
diff --git a/test/numbers.jl b/test/numbers.jl
index a55230465834f..becbee8e8162c 100644
--- a/test/numbers.jl
+++ b/test/numbers.jl
@@ -2389,3 +2389,21 @@ for T in Any[Int16, Int32, UInt32, Int64, UInt64, BigInt]
 end
 
 @test_throws InexactError UInt128(-1)
+
+# multiplicative inverses
+for divisor in 1:254
+    divisor == 128 && continue
+    for T in (UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64)
+        (divisor == 128 && T == Int8) && continue  # abs will fail
+        mi = Base.FastDivInteger(divisor % T)
+        divisorT = divisor % T
+        for dividend in -101:101
+            dividendT = dividend % T
+            @test div(dividendT, mi) == div(dividendT, divisorT)
+            @test rem(dividendT, mi) == rem(dividendT, divisorT)
+            @test divrem(dividendT, mi) == divrem(dividendT, divisorT)
+        end
+    end
+end
+
+@test_throws ErrorException Base.FastDivInteger(0)
diff --git a/test/perf/Makefile b/test/perf/Makefile
index 848dfbcbe3f28..f305a77bd3790 100644
--- a/test/perf/Makefile
+++ b/test/perf/Makefile
@@ -1,9 +1,9 @@
 JULIAHOME = $(abspath ../..)
 include ../../Make.inc
 
-all: micro kernel cat shootout blas lapack simd sort spell sparse
+all: micro kernel cat shootout blas lapack simd sort spell sparse array
 
-micro kernel cat shootout blas lapack simd sort spell sparse:
+micro kernel cat shootout blas lapack simd sort spell sparse array:
 	@$(MAKE) $(QUIET_MAKE) -C shootout
 ifneq ($(OS),WINNT)
 	@$(call spawn,$(JULIA_EXECUTABLE)) $@/perf.jl | perl -nle '@_=split/,/; printf "%-18s %8.3f %8.3f %8.3f %8.3f\n", $$_[1], $$_[2], $$_[3], $$_[4], $$_[5]'
@@ -31,4 +31,4 @@ clean:
 	$(MAKE) -C micro $@
 	$(MAKE) -C shootout $@
 
-.PHONY: micro kernel cat shootout blas lapack simd sort spell sparse clean
+.PHONY: micro kernel cat shootout blas lapack simd sort spell sparse array clean
diff --git a/test/perf/array/indexing.jl b/test/perf/array/indexing.jl
new file mode 100644
index 0000000000000..ed6087328ce2e
--- /dev/null
+++ b/test/perf/array/indexing.jl
@@ -0,0 +1,84 @@
+# Performance testing
+function sumelt(A, n)
+    s = zero(eltype(A))
+    @inbounds for k = 1:n
+        for a in A
+            s += a
+        end
+    end
+    s
+end
+
+function sumeach(A, n)
+    s = zero(eltype(A))
+    for k = 1:n
+        @simd for I in eachindex(A)
+            @inbounds val = A[I]
+            s += val
+        end
+    end
+    s
+end
+
+function sumfast(A, n)
+    s = zero(eltype(A))
+    for k = 1:n
+        @simd for I in fastindex(A)
+            @inbounds val = A[I]
+            s += val
+        end
+    end
+    s
+end
+
+@inline fastindex(A::AbstractArray) = fastindex((A, Base.linearindexing(A)))
+@inline fastindex(AL::(AbstractArray, Base.LinearFast)) = 1:length(AL[1])
+@inline fastindex(AL::(AbstractArray, Base.LinearSlow)) = eachindex(AL[1])
+
+
+abstract MyArray{T,N} <: AbstractArray{T,N}
+
+immutable ArrayLS{T,N} <: MyArray{T,N}  # LinearSlow
+    data::Array{T,N}
+end
+immutable ArrayLF{T,N} <: MyArray{T,N}  # LinearFast
+    data::Array{T,N}
+end
+immutable ArrayStrides{T,N} <: MyArray{T,N}
+    data::Array{T,N}
+    strides::NTuple{N,Int}
+end
+ArrayStrides(A::Array) = ArrayStrides(A, strides(A))
+
+immutable ArrayStrides1{T} <: MyArray{T,2}
+    data::Matrix{T}
+    stride1::Int
+end
+ArrayStrides1(A::Array) = ArrayStrides1(A, size(A,1))
+
+Base.size(A::MyArray) = size(A.data)
+Base.getindex(A::MyArray, indx::Real) = getindex(A.data, indx)
+Base.getindex(A::MyArray, i::Real, j::Real) = getindex(A.data, i, j)
+
+Base.getindex{T}(A::ArrayStrides{T,2}, i::Real, j::Real) = getindex(A.data, 1+A.strides[1]*(i-1)+A.strides[2]*(j-1))
+Base.getindex(A::ArrayStrides1, i::Real, j::Real) = getindex(A.data, i + A.stride1*(j-1))
+
+Base.linearindexing(A::ArrayLF) = Base.LinearFast()
+
+function makearrays{T}(::Type{T}, sz)
+    L = prod(sz)
+    A = reshape(convert(Vector{T}, [1:L;]), sz)
+    AS = ArrayLS(A)
+    AF = ArrayLF(A)
+    Astrd = ArrayStrides(A)
+    Astrd1 = ArrayStrides1(A)
+    outersz = (sz[1]+1,sz[2]+2)
+    B = reshape(convert(Vector{T}, [1:prod(outersz);]), outersz)
+    Asub = sub(B, 1:sz[1], 2:sz[2]+1)
+    RS = reshape(ArrayLS([1:L;]), sz)
+    RF = reshape(ArrayLF([1:L;]), sz)
+    RLS = reshape(ArrayLS(A), (L,))
+    RLF = reshape(ArrayLF(A), (L,))
+    (A, AS, AF, Astrd, Astrd1, Asub, RS, RF, RLS, RLF)
+end
+
diff --git a/test/perf/array/perf.jl b/test/perf/array/perf.jl
new file mode 100644
index 0000000000000..2fea23690d3b5
--- /dev/null
+++ b/test/perf/array/perf.jl
@@ -0,0 +1,38 @@
+include("../perfutil.jl")
+
+include("indexing.jl")
+
+briefname(A) = typeof(A).name.name
+briefname(A::Base.Reshaped.ReshapedArray) = string(typeof(A).name.name, '{', typeof(A.parent).name.name, '}')
+
+# Small array tests
+sz = (3,5)
+Alist = makearrays(Int, sz)
+for Ar in Alist
+    @timeit sumelt(Ar, 10^5) string("sumeltIs ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
+    @timeit sumeach(Ar, 10^5) string("sumeachIs ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
+    @timeit sumfast(Ar, 10^5) string("sumfastIs ", briefname(Ar)) string("for I in fastindex(A), ", briefname(Ar)) sz
+end
+
+Alist = makearrays(Float32, sz)   # SIMD-able
+for Ar in Alist
+    @timeit sumelt(Ar, 10^5) string("sumeltFs ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
+    @timeit sumeach(Ar, 10^5) string("sumeachFs ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
+    @timeit sumfast(Ar, 10^5) string("sumfastFs ", briefname(Ar)) string("for I in fastindex(A), ", briefname(Ar)) sz
+end
+
+# Big array tests
+sz = (300,500)
+Alist = makearrays(Int, sz)
+for Ar in Alist
+    @timeit sumelt(Ar, 100) string("sumeltIb ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
+    @timeit sumeach(Ar, 100) string("sumeachIb ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
+    @timeit sumfast(Ar, 100) string("sumfastIb ", briefname(Ar)) string("for I in fastindex(A), ", briefname(Ar)) sz
+end
+
+Alist = makearrays(Float32, sz)   # SIMD-able
+for Ar in Alist
+    @timeit sumelt(Ar, 100) string("sumeltFb ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
+    @timeit sumeach(Ar, 100) string("sumeachFb ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
+    @timeit sumfast(Ar, 100) string("sumfastFb ", briefname(Ar)) string("for I in fastindex(A), ", briefname(Ar)) sz
+end
diff --git a/test/reshapedarray.jl b/test/reshapedarray.jl
new file mode 100644
index 0000000000000..97ab9ea20a379
--- /dev/null
+++ b/test/reshapedarray.jl
@@ -0,0 +1,111 @@
+module TestReshaped
+
+using Base.Test
+import Base.Reshaped
+
+abstract MyArray{T,N} <: AbstractArray{T,N}
+
+immutable ArrayLS{T,N} <: MyArray{T,N}  # LinearSlow
+    data::Array{T,N}
+end
+immutable ArrayLF{T,N} <: MyArray{T,N}  # LinearFast
+    data::Array{T,N}
+end
+Base.linearindexing(A::ArrayLF) = Base.LinearFast()
+
+Base.size(A::MyArray) = size(A.data)
+Base.getindex(A::MyArray, indx::Real) = getindex(A.data, indx)
+Base.getindex(A::MyArray, i::Real, j::Real) = getindex(A.data, i, j)
+Base.getindex(A::MyArray, i::Real, j::Real, k::Real) = getindex(A.data, i, j, k)
+
+function myeq(A, B)
+    size(A) == size(B) || return false
+    for I in eachindex(B)
+        A[I] == B[I] || return false
+    end
+    true
+end
+indextype{T,N,P,I}(::Reshaped.ReshapedArray{T,N,P,I}) = I
+
+AA = reshape([1:15;], (3, 5))
+CC = reshape([1:16;], (2,4,2))
+@test isa(AA, Array)
+@test isa(CC, Array)
+
+for AT in (ArrayLF, ArrayLS)
+    A = reshape(AT([1:15;]), (3, 5))
+    @test indextype(A) <: (Reshaped.IndexMD{1,2},)
+    @test A == AA
+    @test myeq(A, AA)
+
+    B = reshape(AT(AA), (15,))
+    din = AT == ArrayLF ? 1 : 2
+    @test indextype(B) <: (Reshaped.IndexMD{din,1},)
+    @test B == [1:15;]
+    @test myeq(B, [1:15;])
+
+    CR = reshape(AT(CC), (4,4))
+    din = AT == ArrayLF ? 1 : ndims(CC)
+    @test indextype(CR) <: (Reshaped.IndexMD{din,2},)
+    @test CR == reshape([1:16;], 4, 4)
+    @test myeq(CR, reshape([1:16;], 4, 4))
+    CR2 = reshape(AT(CC), (8,2))
+    @test indextype(CR) <: (Reshaped.IndexMD{din,2},)
+    @test CR2 == reshape([1:16;], 8, 2)
+    @test myeq(CR2, reshape([1:16;], 8, 2))
+
+    V = sub(CR, 2:3,2:4)
+    @test V == [6 10 14; 7 11 15]
+    VR = reshape(V, (3,2))
+    @test indextype(VR) <: (Reshaped.IndexMD{2,2},)
+    @test VR == [6 11; 7 14; 10 15]
+    @test myeq(VR, [6 11; 7 14; 10 15])
+
+    # Cases with dimension sizes of 1, some of which require FastDivInteger1
+    A1 = reshape([1:15;], (1,1,15))
+    R = reshape(AT(A1), (15,))
+    din = AT == ArrayLF ? 1 : ndims(A1)
+    @test indextype(R) <: (Reshaped.IndexMD{din,1},)
+    @test R == [1:15;]
+    @test myeq(R, [1:15;])
+    R = reshape(AT([1:15;]), (1,1,15,))
+    @test indextype(R) <: (Reshaped.IndexMD{1,3},)
+    @test A1 == R
+    @test myeq(A1, R)
+    R = reshape(AT([1:15;]), (15,1))
+    if AT == ArrayLF
+        @test indextype(R) <: (Reshaped.IndexMD{1,2},)
+    else
+        @test indextype(R) == (Colon,Colon)
+    end
+    A2 = reshape([1:15;], 15, 1)
+    @test R == A2
+    @test myeq(R, A2)
+    R = reshape(AT([1:15;]), (3,1,5))
+    @test indextype(R) <: (Reshaped.IndexMD{1,3},)
+    A315 = reshape([1:15;], (3,1,5))
+    @test R == A315
+    @test myeq(R, A315)
+end
+
+A3 = reshape([1:5*7*8;], 5, 7, 8)
+S = sub(A3, 2:5, 1:3:7, 3:5)
+R = reshape(S, (length(S),))
+@test indextype(R) <: (Reshaped.IndexMD{3,1},)
+Sc = copy(S)
+@test R == Sc[:]
+@test myeq(R, Sc[:])
+R = reshape(S, (12,3))
+@test indextype(R) <: (Reshaped.IndexMD{2,1},Colon)
+@test R == reshape(Sc, (12,3))
+@test myeq(R, reshape(Sc, (12,3)))
+R = reshape(S, (4,9))
+@test indextype(R) <: (Colon,Reshaped.IndexMD{2,1})
+@test R == reshape(Sc, (4,9))
+@test myeq(R, reshape(Sc, (4,9)))
+R = reshape(S, (6,6))
+@test indextype(R) <: (Reshaped.IndexMD{3,2},)
+@test R == reshape(Sc, (6,6))
+@test myeq(R, reshape(Sc, (6,6)))
+
+end