diff --git a/base/abstractarray.jl b/base/abstractarray.jl
index 3e82aeb2245fa..0cfea070c3cba 100644
--- a/base/abstractarray.jl
+++ b/base/abstractarray.jl
@@ -6,8 +6,6 @@ typealias AbstractVector{T} AbstractArray{T,1}
 typealias AbstractMatrix{T} AbstractArray{T,2}
 typealias AbstractVecOrMat{T} Union{AbstractVector{T}, AbstractMatrix{T}}
 typealias RangeIndex Union{Int, Range{Int}, AbstractUnitRange{Int}, Colon}
-typealias Indices{N} NTuple{N,AbstractUnitRange}
-typealias IndicesOne{N} NTuple{N,OneTo}
 typealias DimOrInd Union{Integer, AbstractUnitRange}
 typealias DimsOrInds{N} NTuple{N,DimOrInd}
 
@@ -29,6 +27,7 @@ end
 
 size{T,N}(t::AbstractArray{T,N}, d) = d <= N ? size(t)[d] : 1
 size{N}(x, d1::Integer, d2::Integer, dx::Vararg{Integer, N}) = (size(x, d1), size(x, d2), ntuple(k->size(x, dx[k]), Val{N})...)
+
 """
     indices(A, d)
 
@@ -66,7 +65,7 @@ is `indices(A, 1)`.
 Calling this function is the "safe" way to write algorithms that
 exploit linear indexing.
 """
-linearindices(A)                 = (@_inline_meta; 1:length(A))
+linearindices(A)                 = (@_inline_meta; OneTo(_length(A)))
 linearindices(A::AbstractVector) = (@_inline_meta; indices1(A))
 eltype{T}(::Type{AbstractArray{T}}) = T
 eltype{T,N}(::Type{AbstractArray{T,N}}) = T
@@ -74,7 +73,8 @@ elsize{T}(::AbstractArray{T}) = sizeof(T)
 ndims{T,N}(::AbstractArray{T,N}) = N
 ndims{T,N}(::Type{AbstractArray{T,N}}) = N
 ndims{T<:AbstractArray}(::Type{T}) = ndims(supertype(T))
-length(t::AbstractArray) = prod(size(t))::Int
+length(t::AbstractArray) = prod(size(t))
+_length(A::AbstractArray) = prod(map(unsafe_length, indices(A))) # circumvent missing size
 endof(a::AbstractArray) = length(a)
 first(a::AbstractArray) = a[first(eachindex(a))]
 
@@ -132,6 +132,10 @@ function trailingsize(A, n)
     end
     return s
 end
+function trailingsize(inds::Indices)
+    @_inline_meta
+    prod(map(unsafe_length, inds))
+end
 
 ## Traits for array types ##
 
@@ -230,7 +234,7 @@ function checkbounds_indices(::Type{Bool}, IA::Tuple{Any}, I::Tuple{Any})
 end
 function checkbounds_indices(::Type{Bool}, IA::Tuple, I::Tuple{Any})
     @_inline_meta
-    checkindex(Bool, 1:prod(map(dimlength, IA)), I[1])  # linear indexing
+    checkindex(Bool, OneTo(trailingsize(IA)), I[1])  # linear indexing
 end
 
 """
@@ -264,16 +268,6 @@ end
 
 throw_boundserror(A, I) = (@_noinline_meta; throw(BoundsError(A, I)))
 
-@generated function trailingsize{T,N,n}(A::AbstractArray{T,N}, ::Type{Val{n}})
-    (isa(n, Int) && isa(N, Int)) || error("Must have concrete type")
-    n > N && return 1
-    ex = :(size(A, $n))
-    for m = n+1:N
-        ex = :($ex * size(A, $m))
-    end
-    Expr(:block, Expr(:meta, :inline), ex)
-end
-
 # check along a single dimension
 """
     checkindex(Bool, inds::AbstractUnitRange, index)
@@ -357,7 +351,7 @@ to_shape(dims::DimsOrInds) = map(to_shape, dims)
 to_shape(i::Int) = i
 to_shape(i::Integer) = Int(i)
 to_shape(r::OneTo) = Int(last(r))
-to_shape(r::UnitRange) = convert(UnitRange{Int}, r)
+to_shape(r::AbstractUnitRange) = r
 
 """
     similar(storagetype, indices)
@@ -492,7 +486,7 @@ function copy!(::LinearIndexing, dest::AbstractArray, ::LinearSlow, src::Abstrac
 end
 
 function copy!(dest::AbstractArray, dstart::Integer, src::AbstractArray)
-    copy!(dest, dstart, src, first(linearindices(src)), length(src))
+    copy!(dest, dstart, src, first(linearindices(src)), _length(src))
 end
 
 function copy!(dest::AbstractArray, dstart::Integer, src::AbstractArray, sstart::Integer)
@@ -618,7 +612,7 @@ function _maxlength(A, B, C...)
     max(length(A), _maxlength(B, C...))
 end
 
-isempty(a::AbstractArray) = (length(a) == 0)
+isempty(a::AbstractArray) = (_length(a) == 0)
 
 ## Conversions ##
 
diff --git a/base/abstractarraymath.jl b/base/abstractarraymath.jl
index d4b2e48e885d8..43efa5063b48c 100644
--- a/base/abstractarraymath.jl
+++ b/base/abstractarraymath.jl
@@ -11,7 +11,7 @@ transpose(a::AbstractArray) = error("transpose not implemented for $(typeof(a)).
 
 ## Constructors ##
 
-vec(a::AbstractArray) = reshape(a,length(a))
+vec(a::AbstractArray) = reshape(a,_length(a))
 vec(a::AbstractVector) = a
 
 _sub(::Tuple{}, ::Tuple{}) = ()
@@ -74,22 +74,23 @@ function flipdim(A::AbstractArray, d::Integer)
     if d > nd || isempty(A)
         return copy(A)
     end
+    inds = indices(A)
     B = similar(A)
     nnd = 0
     for i = 1:nd
-        nnd += Int(size(A,i)==1 || i==d)
+        nnd += Int(length(inds[i])==1 || i==d)
     end
-    inds = indices(A, d)
-    sd = first(inds)+last(inds)
+    indsd = inds[d]
+    sd = first(indsd)+last(indsd)
     if nnd==nd
         # flip along the only non-singleton dimension
-        for i in inds
+        for i in indsd
             B[i] = A[sd-i]
         end
         return B
     end
     alli = [ indices(B,n) for n in 1:nd ]
-    for i in inds
+    for i in indsd
         B[[ n==d ? sd-i : alli[n] for n in 1:nd ]...] = slicedim(A, d, i)
     end
     return B
diff --git a/base/arraymath.jl b/base/arraymath.jl
index 785735d37bcb0..7ec3eacf6e5b2 100644
--- a/base/arraymath.jl
+++ b/base/arraymath.jl
@@ -252,19 +252,20 @@ end
 
 const transposebaselength=64
 function transpose_f!(f,B::AbstractMatrix,A::AbstractMatrix)
-    indices(B,1) == indices(A,2) && indices(B,2) == indices(A,1) || throw(DimensionMismatch(string(f)))
+    inds = indices(A)
+    indices(B,1) == inds[2] && indices(B,2) == inds[1] || throw(DimensionMismatch(string(f)))
 
-    m, n = size(A)
+    m, n = length(inds[1]), length(inds[2])
     if m*n<=4*transposebaselength
         @inbounds begin
-            for j = indices(A,2)
-                for i = indices(A,1)
+            for j = inds[2]
+                for i = inds[1]
                     B[j,i] = f(A[i,j])
                 end
             end
         end
     else
-        transposeblock!(f,B,A,m,n,first(indices(A,1))-1,first(indices(A,2))-1)
+        transposeblock!(f,B,A,m,n,first(inds[1])-1,first(inds[2])-1)
     end
     return B
 end
diff --git a/base/bitarray.jl b/base/bitarray.jl
index e0754de709666..ecb6efd9ae16a 100644
--- a/base/bitarray.jl
+++ b/base/bitarray.jl
@@ -8,8 +8,7 @@ type BitArray{N} <: DenseArray{Bool, N}
     chunks::Vector{UInt64}
     len::Int
     dims::NTuple{N,Int}
-    function BitArray(dims::Int...)
-        length(dims) == N || throw(ArgumentError("number of dimensions must be $N, got $(length(dims))"))
+    function BitArray(dims::Vararg{Int,N})
         n = 1
         i = 1
         for d in dims
diff --git a/base/broadcast.jl b/base/broadcast.jl
index f3eebb2c0cac3..cca42e514c214 100644
--- a/base/broadcast.jl
+++ b/base/broadcast.jl
@@ -3,7 +3,7 @@
 module Broadcast
 
 using Base.Cartesian
-using Base: promote_op, promote_eltype, promote_eltype_op, @get!, _msk_end, unsafe_bitgetindex, linearindices, to_shape, tail, dimlength, OneTo
+using Base: promote_op, promote_eltype, promote_eltype_op, @get!, _msk_end, unsafe_bitgetindex, linearindices, tail, OneTo, to_shape
 import Base: .+, .-, .*, ./, .\, .//, .==, .<, .!=, .<=, .÷, .%, .<<, .>>, .^
 export broadcast, broadcast!, bitbroadcast
 export broadcast_getindex, broadcast_setindex!
@@ -20,7 +20,7 @@ broadcast(f, x::Number...) = f(x...)
 broadcast_shape() = ()
 broadcast_shape(A) = indices(A)
 @inline broadcast_shape(A, B...) = broadcast_shape((), indices(A), map(indices, B)...)
-# shape inputs
+# shape (i.e., tuple-of-indices) inputs
 broadcast_shape(shape::Tuple) = shape
 @inline broadcast_shape(shape::Tuple, shape1::Tuple, shapes::Tuple...) = broadcast_shape(_bcs((), shape, shape1), shapes...)
 # _bcs consolidates two shapes into a single output shape
@@ -37,7 +37,7 @@ _bcs1(a::Integer, b) = a == 1 ? b : (first(b) == 1 && last(b) == a ? b : throw(D
 _bcs1(a, b::Integer) = _bcs1(b, a)
 _bcs1(a, b) = _bcsm(b, a) ? b : (_bcsm(a, b) ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size")))
 # _bcsm tests whether the second index is consistent with the first
-_bcsm(a, b) = a == b || (dimlength(b) == 1 && first(b) == first(a))
+_bcsm(a, b) = a == b || length(b) == 1
 _bcsm(a, b::Number) = b == 1
 _bcsm(a::Number, b::Number) = a == b || b == 1
 
@@ -62,22 +62,37 @@ end
 end
 
 ## Indexing manipulations
-# newindex(I, rule) replaces a CartesianIndex with something that is
-# appropriate for a particular array/scalar. `rule` is a tuple that
-# describes the manipulations that should be made.
-@inline newindex(I::CartesianIndex, ::Tuple{}) = 1    # for scalars
-@inline newindex(I::CartesianIndex, indexmap) = CartesianIndex(_newindex((), I.I, indexmap...))
-@inline _newindex(out, I) = out  # can truncate if indexmap is shorter than I
-@inline _newindex(out, I, keep::Bool, indexmap...) = _newindex((out..., ifelse(keep, I[1], 1)), tail(I), indexmap...)
-
-newindexer(sz, x::Number) = ()
-@inline newindexer(sz, A) = _newindexer(sz, size(A))
-@inline _newindexer(sz, szA::Tuple{}) = ()
-@inline _newindexer(sz, szA) = (sz[1] == szA[1], _newindexer(tail(sz), tail(szA))...)
-
-# map(x->newindexer(sz, x), As), but see #15276
-map_newindexer(sz, ::Tuple{}) = ()
-@inline map_newindexer(sz, As) = (newindexer(sz, As[1]), map_newindexer(sz, tail(As))...)
+
+# newindex(I, keep, Idefault) replaces a CartesianIndex `I` with something that
+# is appropriate for a particular broadcast array/scalar. `keep` is a
+# NTuple{N,Bool}, where keep[d] == true means that one should preserve
+# I[d]; if false, replace it with Idefault[d].
+@inline newindex(I::CartesianIndex, keep, Idefault) = CartesianIndex(_newindex(I.I, keep, Idefault))
+@inline _newindex(I, keep, Idefault) =
+    (ifelse(keep[1], I[1], Idefault[1]), _newindex(tail(I), tail(keep), tail(Idefault))...)
+@inline _newindex(I, keep::Tuple{}, Idefault) = ()  # truncate if keep is shorter than I
+
+# newindexer(shape, A) generates `keep` and `Idefault` (for use by
+# `newindex` above) for a particular array `A`, given the
+# broadcast_shape `shape`
+# `keep` is equivalent to map(==, indices(A), shape) (but see #17126)
+newindexer(shape, x::Number) = (), ()
+@inline newindexer(shape, A) = newindexer(shape, indices(A))
+@inline newindexer(shape, indsA::Tuple{}) = (), ()
+@inline function newindexer(shape, indsA::Tuple)
+    ind1 = indsA[1]
+    keep, Idefault = newindexer(tail(shape), tail(indsA))
+    (shape[1] == ind1, keep...), (first(ind1), Idefault...)
+end
+
+# Equivalent to map(x->newindexer(shape, x), As) (but see #17126)
+map_newindexer(shape, ::Tuple{}) = (), ()
+@inline function map_newindexer(shape, As)
+    A1 = As[1]
+    keeps, Idefaults = map_newindexer(shape, tail(As))
+    keep, Idefault = newindexer(shape, A1)
+    (keep, keeps...), (Idefault, Idefaults...)
+end
 
 # For output BitArrays
 const bitcache_chunks = 64 # this can be changed
@@ -88,16 +103,17 @@ dumpbitcache(Bc::Vector{UInt64}, bind::Int, C::Vector{Bool}) =
 
 ## Broadcasting core
 # nargs encodes the number of As arguments (which matches the number
-# of indexmaps). The first two type parameters are to ensure specialization.
-@generated function _broadcast!{M,AT,nargs}(f, B::AbstractArray, indexmaps::M, As::AT, ::Type{Val{nargs}})
+# of keeps). The first two type parameters are to ensure specialization.
+@generated function _broadcast!{K,ID,AT,nargs}(f, B::AbstractArray, keeps::K, Idefaults::ID, As::AT, ::Type{Val{nargs}})
     quote
         $(Expr(:meta, :noinline))
-        # destructure the indexmaps and As tuples
+        # destructure the keeps and As tuples
         @nexprs $nargs i->(A_i = As[i])
-        @nexprs $nargs i->(imap_i = indexmaps[i])
+        @nexprs $nargs i->(keep_i = keeps[i])
+        @nexprs $nargs i->(Idefault_i = Idefaults[i])
         @simd for I in CartesianRange(indices(B))
             # reverse-broadcast the indices
-            @nexprs $nargs i->(I_i = newindex(I, imap_i))
+            @nexprs $nargs i->(I_i = newindex(I, keep_i, Idefault_i))
             # extract array values
             @nexprs $nargs i->(@inbounds val_i = A_i[I_i])
             # call the function and store the result
@@ -108,19 +124,20 @@ end
 
 # For BitArray outputs, we cache the result in a "small" Vector{Bool},
 # and then copy in chunks into the output
-@generated function _broadcast!{M,AT,nargs}(f, B::BitArray, indexmaps::M, As::AT, ::Type{Val{nargs}})
+@generated function _broadcast!{K,ID,AT,nargs}(f, B::BitArray, keeps::K, Idefaults::ID, As::AT, ::Type{Val{nargs}})
     quote
         $(Expr(:meta, :noinline))
-        # destructure the indexmaps and As tuples
+        # destructure the keeps and As tuples
         @nexprs $nargs i->(A_i = As[i])
-        @nexprs $nargs i->(imap_i = indexmaps[i])
+        @nexprs $nargs i->(keep_i = keeps[i])
+        @nexprs $nargs i->(Idefault_i = Idefaults[i])
         C = Vector{Bool}(bitcache_size)
         Bc = B.chunks
         ind = 1
         cind = 1
         @simd for I in CartesianRange(indices(B))
             # reverse-broadcast the indices
-            @nexprs $nargs i->(I_i = newindex(I, imap_i))
+            @nexprs $nargs i->(I_i = newindex(I, keep_i, Idefault_i))
             # extract array values
             @nexprs $nargs i->(@inbounds val_i = A_i[I_i])
             # call the function and store the result
@@ -140,25 +157,26 @@ end
 end
 
 @inline function broadcast!{nargs}(f, B::AbstractArray, As::Vararg{Any,nargs})
-    check_broadcast_shape(indices(B), As...)
-    sz = size(B)
-    mapindex = map(x->newindexer(sz, x), As)
-    _broadcast!(f, B, mapindex, As, Val{nargs})
+    shape = indices(B)
+    check_broadcast_shape(shape, As...)
+    keeps, Idefaults = map_newindexer(shape, As)
+    _broadcast!(f, B, keeps, Idefaults, As, Val{nargs})
     B
 end
 
 # broadcast with computed element type
 
-@generated function _broadcast!{M,AT,nargs}(f, B::AbstractArray, indexmaps::M, As::AT, ::Type{Val{nargs}}, iter, st, count)
+@generated function _broadcast!{K,ID,AT,nargs}(f, B::AbstractArray, keeps::K, Idefaults::ID, As::AT, ::Type{Val{nargs}}, iter, st, count)
     quote
         $(Expr(:meta, :noinline))
-        # destructure the indexmaps and As tuples
+        # destructure the keeps and As tuples
         @nexprs $nargs i->(A_i = As[i])
-        @nexprs $nargs i->(imap_i = indexmaps[i])
+        @nexprs $nargs i->(keep_i = keeps[i])
+        @nexprs $nargs i->(Idefault_i = Idefaults[i])
         while !done(iter, st)
             I, st = next(iter, st)
             # reverse-broadcast the indices
-            @nexprs $nargs i->(I_i = newindex(I, imap_i))
+            @nexprs $nargs i->(I_i = newindex(I, keep_i, Idefault_i))
             # extract array values
             @nexprs $nargs i->(@inbounds val_i = A_i[I_i])
             # call the function
@@ -174,7 +192,7 @@ end
                     new[II] = B[II]
                 end
                 new[I] = V
-                return _broadcast!(f, new, indexmaps, As, Val{nargs}, iter, st, count+1)
+                return _broadcast!(f, new, keeps, Idefaults, As, Val{nargs}, iter, st, count+1)
             end
             count += 1
         end
@@ -183,20 +201,19 @@ end
 end
 
 function broadcast_t(f, ::Type{Any}, As...)
-    shp = broadcast_shape(As...)
-    iter = CartesianRange(shp)
+    shape = broadcast_shape(As...)
+    iter = CartesianRange(shape)
     if isempty(iter)
-        return similar(Array{Union{}}, shp)
+        return similar(Array{Union{}}, shape)
     end
     nargs = length(As)
-    sz = size(iter)
-    indexmaps = map(x->newindexer(sz, x), As)
+    keeps, Idefaults = map_newindexer(shape, As)
     st = start(iter)
     I, st = next(iter, st)
-    val = f([ As[i][newindex(I, indexmaps[i])] for i=1:nargs ]...)
-    B = similar(Array{typeof(val)}, shp)
+    val = f([ As[i][newindex(I, keeps[i], Idefaults[i])] for i=1:nargs ]...)
+    B = similar(Array{typeof(val)}, shape)
     B[I] = val
-    return _broadcast!(f, B, indexmaps, As, Val{nargs}, iter, st, 1)
+    return _broadcast!(f, B, keeps, Idefaults, As, Val{nargs}, iter, st, 1)
 end
 
 @inline broadcast_t(f, T, As...) = broadcast!(f, similar(Array{T}, broadcast_shape(As...)), As...)
@@ -206,25 +223,24 @@ end
 # alternate, more compact implementation; unfortunately slower.
 # also the `collect` machinery doesn't yet support arbitrary index bases.
 #=
-@generated function _broadcast{nargs}(f, indexmaps, As, ::Type{Val{nargs}}, iter)
+@generated function _broadcast{nargs}(f, keeps, As, ::Type{Val{nargs}}, iter)
     quote
-        collect((@ncall $nargs f i->As[i][newindex(I, indexmaps[i])]) for I in iter)
+        collect((@ncall $nargs f i->As[i][newindex(I, keeps[i])]) for I in iter)
     end
 end
 
 function broadcast(f, As...)
-    shp = broadcast_shape(As...)
-    iter = CartesianRange(shp)
-    sz = size(iter)
-    indexmaps = map(x->newindexer(sz, x), As)
+    shape = broadcast_shape(As...)
+    iter = CartesianRange(shape)
+    keeps, Idefaults = map_newindexer(shape, As)
     naT = Val{nfields(As)}
-    _broadcast(f, indexmaps, As, naT, iter)
+    _broadcast(f, keeps, Idefaults, As, naT, iter)
 end
 =#
 
 @inline bitbroadcast(f, As...) = broadcast!(f, similar(BitArray, broadcast_shape(As...)), As...)
 
-broadcast_getindex(src::AbstractArray, I::AbstractArray...) = broadcast_getindex!(Array{eltype(src)}(to_shape(broadcast_shape(I...))), src, I...)
+broadcast_getindex(src::AbstractArray, I::AbstractArray...) = broadcast_getindex!(similar(Array{eltype(src)}, broadcast_shape(I...)), src, I...)
 @generated function broadcast_getindex!(dest::AbstractArray, src::AbstractArray, I::AbstractArray...)
     N = length(I)
     Isplat = Expr[:(I[$d]) for d = 1:N]
@@ -232,7 +248,8 @@ broadcast_getindex(src::AbstractArray, I::AbstractArray...) = broadcast_getindex
         @nexprs $N d->(I_d = I[d])
         check_broadcast_shape(indices(dest), $(Isplat...))  # unnecessary if this function is never called directly
         checkbounds(src, $(Isplat...))
-        @nloops $N i dest d->(@nexprs $N k->(j_d_k = size(I_k, d) == 1 ? 1 : i_d)) begin
+        @nexprs $N d->(@nexprs $N k->(Ibcast_d_k = indices(I_k, d) == OneTo(1)))
+        @nloops $N i dest d->(@nexprs $N k->(j_d_k = Ibcast_d_k ? 1 : i_d)) begin
             @nexprs $N k->(@inbounds J_k = @nref $N I_k d->j_d_k)
             @inbounds (@nref $N dest i) = (@nref $N src J)
         end
@@ -248,18 +265,19 @@ end
         checkbounds(A, $(Isplat...))
         shape = broadcast_shape($(Isplat...))
         @nextract $N shape d->(length(shape) < d ? OneTo(1) : shape[d])
+        @nexprs $N d->(@nexprs $N k->(Ibcast_d_k = indices(I_k, d) == 1:1))
         if !isa(x, AbstractArray)
             xA = convert(eltype(A), x)
-            @nloops $N i d->shape_d d->(@nexprs $N k->(j_d_k = size(I_k, d) == 1 ? 1 : i_d)) begin
+            @nloops $N i d->shape_d d->(@nexprs $N k->(j_d_k = Ibcast_d_k ? 1 : i_d)) begin
                 @nexprs $N k->(@inbounds J_k = @nref $N I_k d->j_d_k)
                 @inbounds (@nref $N A J) = xA
             end
         else
             X = x
-            @nexprs $N d->(shapelen_d = dimlength(shape_d))
+            @nexprs $N d->(shapelen_d = length(shape_d))
             @ncall $N Base.setindex_shape_check X shapelen
             Xstate = start(X)
-            @inbounds @nloops $N i d->shape_d d->(@nexprs $N k->(j_d_k = size(I_k, d) == 1 ? 1 : i_d)) begin
+            @inbounds @nloops $N i d->shape_d d->(@nexprs $N k->(j_d_k = Ibcast_d_k ? 1 : i_d)) begin
                 @nexprs $N k->(J_k = @nref $N I_k d->j_d_k)
                 x_el, Xstate = next(X, Xstate)
                 (@nref $N A J) = x_el
@@ -290,11 +308,11 @@ end
 function broadcast_bitarrays(scalarf, bitf, A::AbstractArray{Bool}, B::AbstractArray{Bool})
     local shape
     try
-        shape = promote_shape(size(A), size(B))
+        shape = promote_shape(indices(A), indices(B))
     catch
         return bitbroadcast(scalarf, A, B)
     end
-    F = BitArray(shape)
+    F = BitArray(to_shape(shape))
     Fc = F.chunks
     Ac = BitArray(A).chunks
     Bc = BitArray(B).chunks
@@ -381,11 +399,11 @@ end
 function (.^){T<:Integer}(A::BitArray, B::Array{T})
     local shape
     try
-        shape = promote_shape(size(A), size(B))
+        shape = promote_shape(indices(A), indices(B))
     catch
         return bitbroadcast(^, A, B)
     end
-    F = BitArray(shape)
+    F = BitArray(to_shape(shape))
     l = length(F)
     l == 0 && return F
     Ac = A.chunks
diff --git a/base/deprecated.jl b/base/deprecated.jl
index 9f0247a809e21..816bcea7f0b66 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -752,16 +752,26 @@ function checkbounds{N,T}(::Type{Bool}, sz::NTuple{N,Integer}, I1::T, I...)
 end
 
 function first(::Colon)
-    depwarn("first(:) is no longer unambiguous, call Base._first(:, A, dim)", :first)
+    depwarn("first(:) is deprecated, see http://docs.julialang.org/en/latest/devdocs/offset-arrays/", :first)
     1
 end
+function _first(i, A, d)
+    depwarn("_first is deprecated, see http://docs.julialang.org/en/latest/devdocs/offset-arrays/", :_first)
+    __first(i, A, d)
+end
+__first(::Colon, P, ::Colon) = first(linearindices(P))
+__first(i, P, ::Colon) = first(i)
+__first(::Colon, P, d) = first(indices(P, d))
+__first(i, P, d) = first(i)
 
-# Not exported, but may be useful just in case
+# Not exported, but deprecation may be useful just in case
 function Broadcast.check_broadcast_shape(sz::Dims, As::Union{AbstractArray,Number}...)
     depwarn("check_broadcast_shape(size(A), B...) should be replaced with check_broadcast_shape(indices(A), B...)", :check_broadcast_shape)
     Broadcast.check_broadcast_shape(map(OneTo, sz), As...)
 end
 
+@deprecate trailingsize{n}(A::AbstractArray, ::Type{Val{n}}) trailingsize(A, n)
+
 @deprecate slice view
 @deprecate sub view
 
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 199ef18cf5afc..959034ccdc374 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -183,13 +183,25 @@ function checkindex{N}(::Type{Bool}, inds::Tuple, I::AbstractArray{CartesianInde
     b
 end
 
+# combined dimensionality of all indices, including CartesianIndex and
+# AbstractArray{CartesianIndex}
+# rather than returning N, it returns an NTuple{N,Bool} so the result is inferrable
+@inline index_ndims(i1, I...) = (true, index_ndims(I...)...)
+@inline function index_ndims{N}(i1::CartesianIndex{N}, I...)
+    (map(x->true, i1.I)..., index_ndims(I...)...)
+end
+@inline function index_ndims{N}(i1::AbstractArray{CartesianIndex{N}}, I...)
+    (ntuple(x->true, Val{N})..., index_ndims(I...)...)
+end
+index_ndims() = ()
+
 # Recursively compute the lengths of a list of indices, without dropping scalars
 # These need to be inlined for more than 3 indexes
 index_lengths(A::AbstractArray, I::Colon) = (length(A),)
 @inline index_lengths(A::AbstractArray, I...) = index_lengths_dim(A, 1, I...)
 index_lengths_dim(A, dim) = ()
 index_lengths_dim(A, dim, ::Colon) = (trailingsize(A, dim),)
-@inline index_lengths_dim(A, dim, ::Colon, i, I...) = (size(A, dim), index_lengths_dim(A, dim+1, i, I...)...)
+@inline index_lengths_dim(A, dim, ::Colon, i, I...) = (_length(indices(A, dim)), index_lengths_dim(A, dim+1, i, I...)...)
 @inline index_lengths_dim(A, dim, ::Real, I...) = (1, index_lengths_dim(A, dim+1, I...)...)
 @inline index_lengths_dim{N}(A, dim, ::CartesianIndex{N}, I...) = (1, index_lengths_dim(A, dim+N, I...)...)
 @inline index_lengths_dim(A, dim, i::AbstractArray, I...) = (length(i), index_lengths_dim(A, dim+1, I...)...)
@@ -197,34 +209,49 @@ index_lengths_dim(A, dim, ::Colon) = (trailingsize(A, dim),)
 @inline index_lengths_dim{N}(A, dim, i::AbstractArray{CartesianIndex{N}}, I...) = (length(i), index_lengths_dim(A, dim+N, I...)...)
 
 # shape of array to create for getindex() with indexes I, dropping scalars
-# Rather than use an Integer dim, we grow a tuple (true, true, ...)
-# whose length is equal to the dimension we're to process next. This
-# allows us to dispatch, which is important for the type-stability of
-# the lines involving Colon as the final index.
-index_shape(A::AbstractVector, I::Colon) = indices(A)
-index_shape(A::AbstractArray,  I::Colon) = (length(A),)
-@inline index_shape(A::AbstractArray, I...) = index_shape_dim(A, (true,), I...)
-@inline index_shape_dim(A, dim, ::Colon) = (trailingsize(A, length(dim)),)
-@inline index_shape_dim{T,N}(A::AbstractArray{T,N}, dim::NTuple{N}, ::Colon) = (indices(A, N),)
-@inline index_shape_dim(A, dim, I::Real...) = ()
-@inline index_shape_dim(A, dim, ::Colon, i, I...) = (indices(A, length(dim)), index_shape_dim(A, (dim...,true), i, I...)...)
-@inline index_shape_dim(A, dim, ::Real, I...) = (index_shape_dim(A, (dim...,true), I...)...)
-@inline index_shape_dim{N}(A, dim, ::CartesianIndex{N}, I...) = (index_shape_dim(A, (dim...,ntuple(d->true,Val{N})...), I...)...)
-@inline index_shape_dim(A, dim, i::AbstractArray, I...) = (indices(i)..., index_shape_dim(A, (dim...,true), I...)...)
-@inline index_shape_dim(A, dim, i::AbstractArray{Bool}, I...) = (sum(i), index_shape_dim(A, (dim...,true), I...)...)
-@inline index_shape_dim{N}(A, dim, i::AbstractArray{CartesianIndex{N}}, I...) = (indices(i)..., index_shape_dim(A, (dim...,ntuple(d->true,Val{N})...), I...)...)
-
-@inline decolon(A::AbstractVector, ::Colon) = (indices(A,1),)
-@inline decolon(A::AbstractArray,  ::Colon) = (1:length(A),)
-@inline decolon(A::AbstractArray, I...)     = decolon_dim(A, (true,), I...)
-@inline decolon_dim(A::AbstractArray, dim)  = ()
-@inline decolon_dim{T,N}(A::AbstractArray{T,N}, dim::NTuple{N}, ::Colon) = (indices(A, N),)
-@inline decolon_dim(A, dim, ::Colon)        = (1:trailingsize(A, length(dim)),)
-@inline decolon_dim(A::AbstractArray, dim, i1, I...)      = (i1, decolon_dim(A, (dim...,true), I...)...)
-@inline decolon_dim{N}(A::AbstractArray, dim, i1::AbstractArray{CartesianIndex{N}}, I...) = (i1, decolon_dim(A, (dim...,ntuple(d->true,Val{N})...), I...)...)
-@inline decolon_dim(A::AbstractArray, dim, ::Colon, I...) = (indices(A, length(dim)), decolon_dim(A, (dim...,true), I...)...)
+# returns a Tuple{Vararg{AbstractUnitRange}} of indices
+index_shape(A::AbstractArray,  I::Colon)    = (linearindices(A),)
+@inline index_shape(A::AbstractArray, I...) = index_shape_dim(indices(A), I...)
+@inline index_shape_dim(inds::Tuple{Any}, ::Colon)          = inds
+@inline index_shape_dim(inds,             ::Colon)          = (OneTo(trailingsize(inds)),)
+@inline function index_shape_dim(inds,    ::Colon, i, I...)
+    inds1, indstail = IteratorsMD.split(inds, Val{1})
+    (inds1..., index_shape_dim(indstail, i, I...)...)
+end
+@inline index_shape_dim(inds,  ::Real...)             = ()
+@inline index_shape_dim(inds,  ::Real, I...)          = index_shape_dim(safe_tail(inds), I...)
+@inline index_shape_dim(inds, i::AbstractArray, I...) =
+    (indices(i)..., index_shape_dim(safe_tail(inds), I...)...)
+@inline index_shape_dim(inds, i::AbstractArray{Bool}, I...) =
+    (OneTo(sum(i)), index_shape_dim(safe_tail(inds), I...)...)
+# single CartesianIndex version not needed because of call to flatten in _getindex...
+# ...but array of CartesianIndex is not covered
+@inline function index_shape_dim{N}(inds, i::AbstractArray{CartesianIndex{N}}, I...)
+    indsN, indstail = IteratorsMD.split(inds, Val{N})
+    (indices(i)..., index_shape_dim(indstail, I...)...)
+end
+
+# Convert Colon indices into explicit indices
+@inline decolon(A::AbstractArray, ::Colon) = (linearindices(A),)
+@inline decolon(A::AbstractArray, I...) = decolon_dim(indices(A), I...)
+@inline decolon_dim(inds)  = ()
+@inline decolon_dim(inds::Tuple{Any}, ::Colon)       = inds
+@inline decolon_dim(inds,             ::Colon)       = (OneTo(trailingsize(inds)),)
+@inline function decolon_dim(inds,    ::Colon, I...)
+    inds1, indstail = IteratorsMD.split(inds, Val{1})
+    (maybe_oneto(inds1...), decolon_dim(indstail, I...)...)
+end
+@inline decolon_dim(inds, i1, I...) = (i1, decolon_dim(safe_tail(inds), I...)...)
+@inline function decolon_dim{N}(inds, i1::AbstractArray{CartesianIndex{N}}, I...)
+    indsN, indstail = IteratorsMD.split(inds, Val{N})
+    (i1, decolon_dim(indstail, I...)...)
+end
+maybe_oneto(i) = i
+maybe_oneto() = OneTo(1)
 
 ### From abstractarray.jl: Internal multidimensional indexing definitions ###
+getindex(x::Number, i::CartesianIndex{0}) = x
+
 # These are not defined on directly on getindex to avoid
 # ambiguities for AbstractArray subtypes. See the note in abstractarray.jl
 
@@ -236,16 +263,18 @@ end
 # Explicitly allow linear indexing with one non-scalar index
 @inline function _getindex(l::LinearIndexing, A::AbstractArray, i::Union{Real, AbstractArray, Colon})
     @boundscheck checkbounds(A, i)
-    _unsafe_getindex(l, _maybe_linearize(l, A), i)
+    _unsafe_getindex(l, _maybe_reshape(l, A, (i,)), i)
 end
 # But we can speed up LinearSlow arrays by reshaping them to vectors:
-_maybe_linearize(::LinearFast, A::AbstractArray) = A
-_maybe_linearize(::LinearSlow, A::AbstractVector) = A
-_maybe_linearize(::LinearSlow, A::AbstractArray) = reshape(A, length(A))
+_maybe_reshape(::LinearFast, A::AbstractArray, i) = A
+_maybe_reshape(::LinearSlow, A::AbstractVector, i) = A
+@inline _maybe_reshape(::LinearSlow, A::AbstractArray, i) = _maybe_reshape(LinearSlow(), index_ndims(i...), A)
+@inline _maybe_reshape{T,N}(::LinearIndexing, ::NTuple{N}, A::AbstractArray{T,N}) = A
+@inline _maybe_reshape{N}(::LinearIndexing, ::NTuple{N}, A) = reshape(A, Val{N})
 
 @inline function _getindex{N}(l::LinearIndexing, A::AbstractArray, I::Vararg{Union{Real, AbstractArray, Colon},N}) # TODO: DEPRECATE FOR #14770
     @boundscheck checkbounds(A, I...)
-    _unsafe_getindex(l, reshape(A, Val{N}), I...)
+    _unsafe_getindex(l, _maybe_reshape(l, A, I), I...)
 end
 
 @generated function _unsafe_getindex(::LinearIndexing, A::AbstractArray, I::Union{Real, AbstractArray, Colon}...)
@@ -255,7 +284,7 @@ end
         @nexprs $N d->(I_d = to_index(I[d]))
         shape = @ncall $N index_shape A I
         dest = similar(A, shape)
-        size(dest) == map(dimlength, shape) || throw_checksize_error(dest, shape)
+        map(unsafe_length, indices(dest)) == map(unsafe_length, shape) || throw_checksize_error(dest, shape)
         @ncall $N _unsafe_getindex! dest A I
     end
 end
@@ -264,7 +293,7 @@ end
 function _unsafe_getindex(::LinearIndexing, src::AbstractArray, I::AbstractArray{Bool})
     shape = index_shape(src, I)
     dest = similar(src, shape)
-    size(dest) == map(dimlength, shape) || throw_checksize_error(dest, shape)
+    map(unsafe_length, indices(dest)) == map(unsafe_length, shape) || throw_checksize_error(dest, shape)
 
     D = eachindex(dest)
     Ds = start(D)
@@ -281,7 +310,7 @@ end
 function _unsafe_getindex(::LinearFast, src::AbstractArray, I::AbstractArray{Bool})
     shape = index_shape(src, I)
     dest = similar(src, shape)
-    size(dest) == shape || throw_checksize_error(dest, shape)
+    map(unsafe_length, indices(dest)) == map(unsafe_length, shape) || throw_checksize_error(dest, shape)
 
     D = eachindex(dest)
     Ds = start(D)
@@ -328,12 +357,12 @@ _iterable(v) = repeated(v)
 end
 @inline function _setindex!(l::LinearIndexing, A::AbstractArray, x, j::Union{Real,AbstractArray,Colon})
     @boundscheck checkbounds(A, j)
-    _unsafe_setindex!(l, _maybe_linearize(l, A), x, j)
+    _unsafe_setindex!(l, _maybe_reshape(l, A, (j,)), x, j)
     A
 end
 @inline function _setindex!{N}(l::LinearIndexing, A::AbstractArray, x, J::Vararg{Union{Real, AbstractArray, Colon},N}) # TODO: DEPRECATE FOR #14770
     @boundscheck checkbounds(A, J...)
-    _unsafe_setindex!(l, reshape(A, Val{N}), x, J...)
+    _unsafe_setindex!(l, _maybe_reshape(l, A, J), x, J...)
     A
 end
 
@@ -433,21 +462,26 @@ for (f, fmod, op) = ((:cummin, :_cummin!, :min), (:cummax, :_cummax!, :max))
     end
 
     @eval function ($f)(A::AbstractArray, axis::Integer)
+        axis > 0 || throw(ArgumentError("axis must be a positive integer"))
         res = similar(A)
-        if size(A, axis) < 1
+        axis > ndims(A) && return copy!(res, A)
+        inds = indices(A)
+        if isempty(inds[axis])
             return res
         end
-        R1 = CartesianRange(size(A)[1:axis-1])
-        R2 = CartesianRange(size(A)[axis+1:end])
+        R1 = CartesianRange(inds[1:axis-1])
+        R2 = CartesianRange(inds[axis+1:end])
         ($fmod)(res, A, R1, R2, axis)
     end
 
     @eval @noinline function ($fmod)(res, A::AbstractArray, R1::CartesianRange, R2::CartesianRange, axis::Integer)
+        inds = indices(A, axis)
+        i1 = first(inds)
         for I2 in R2
             for I1 in R1
-                res[I1, 1, I2] = A[I1, 1, I2]
+                res[I1, i1, I2] = A[I1, i1, I2]
             end
-            for i = 2:size(A, axis)
+            for i = i1+1:last(inds)
                 for I1 in R1
                     res[I1, i, I2] = ($op)(A[I1, i, I2], res[I1, i-1, I2])
                 end
@@ -468,18 +502,16 @@ cumsum!(B, A, axis::Integer) = cumop!(+, B, A, axis)
 cumprod!(B, A, axis::Integer) = cumop!(*, B, A, axis)
 
 function cumop!(op, B, A, axis::Integer)
-    if size(B, axis) < 1
-        return B
-    end
-    indices(B) == indices(A) || throw(DimensionMismatch("Shape of B must match A"))
-    if axis > ndims(A)
-        copy!(B, A)
-        return B
-    end
+    axis > 0 || throw(ArgumentError("axis must be a positive integer"))
+    inds_t = indices(A)
+    indices(B) == inds_t || throw(DimensionMismatch("shape of B must match A"))
+    axis > ndims(A) && return copy!(B, A)
+    isempty(inds_t[axis]) && return B
     if axis == 1
-        # We can accumulate to a temporary variable, which allows register usage and will be slightly faster
-        ind1 = indices(A,1)
-        @inbounds for I in CartesianRange(tail(indices(A)))
+        # We can accumulate to a temporary variable, which allows
+        # register usage and will be slightly faster
+        ind1 = inds_t[1]
+        @inbounds for I in CartesianRange(tail(inds_t))
             tmp = convert(eltype(B), A[first(ind1), I])
             B[first(ind1), I] = tmp
             for i_1 = first(ind1)+1:last(ind1)
@@ -490,7 +522,7 @@ function cumop!(op, B, A, axis::Integer)
     else
         R1 = CartesianRange(indices(A)[1:axis-1])   # not type-stable
         R2 = CartesianRange(indices(A)[axis+1:end])
-        _cumop!(op, B, A, R1, indices(A, axis), R2) # use function barrier
+        _cumop!(op, B, A, R1, inds_t[axis], R2) # use function barrier
     end
     return B
 end
@@ -534,7 +566,7 @@ end
 # contiguous multidimensional indexing: if the first dimension is a range,
 # we can get some performance from using copy_chunks!
 @inline function _unsafe_getindex!(X::BitArray, B::BitArray, I0::Union{UnitRange{Int},Colon})
-    copy_chunks!(X.chunks, 1, B.chunks, _first(I0, B, :), index_lengths(B, I0)[1])
+    copy_chunks!(X.chunks, 1, B.chunks, indexoffset(I0)+1, index_lengths(B, I0)[1])
     return X
 end
 
@@ -545,7 +577,7 @@ end
         $(Expr(:meta, :inline))
         @nexprs $N d->(I_d = I[d])
 
-        f0 = _first(I0, B, 1)
+        f0 = indexoffset(I0)+1
         l0 = size(X, 1)
 
         gap_lst_1 = 0
@@ -555,7 +587,7 @@ end
         @nexprs $N d->begin
             stride *= size(B, d)
             stride_lst_d = stride
-            ind += stride * (_first(I_d, B, d) - 1)
+            ind += stride * indexoffset(I_d)
             gap_lst_{d+1} *= stride
         end
 
@@ -604,7 +636,7 @@ end
     l0 = index_lengths(B, I0)[1]
     setindex_shape_check(X, l0)
     l0 == 0 && return B
-    f0 = _first(I0, B, :)
+    f0 = indexoffset(I0)+1
     copy_to_bitarray_chunks!(B.chunks, f0, X, 1, l0)
     return B
 end
@@ -614,7 +646,7 @@ end
     y = Bool(x)
     l0 = index_lengths(B, I0)[1]
     l0 == 0 && return B
-    f0 = _first(I0, B, :)
+    f0 = indexoffset(I0)+1
     fill_chunks!(B.chunks, y, f0, l0)
     return B
 end
@@ -630,7 +662,7 @@ end
         idxlens = @ncall $N index_lengths B I0 d->I[d]
         @ncall $N setindex_shape_check X idxlens[1] d->idxlens[d+1]
         isempty(X) && return B
-        f0 = _first(I0, B, 1)
+        f0 = indexoffset(I0)+1
         l0 = idxlens[1]
 
         gap_lst_1 = 0
@@ -640,7 +672,7 @@ end
         @nexprs $N d->begin
             stride *= size(B, d)
             stride_lst_d = stride
-            ind += stride * (_first(I[d], B, d) - 1)
+            ind += stride * indexoffset(I[d])
             gap_lst_{d+1} *= stride
         end
 
@@ -669,7 +701,7 @@ end
         y = Bool(x)
         idxlens = @ncall $N index_lengths B I0 d->I[d]
 
-        f0 = _first(I0, B, 1)
+        f0 = indexoffset(I0)+1
         l0 = idxlens[1]
         l0 == 0 && return B
         @nexprs $N d->(isempty(I[d]) && return B)
@@ -681,7 +713,7 @@ end
         @nexprs $N d->begin
             stride *= size(B, d)
             stride_lst_d = stride
-            ind += stride * (_first(I[d], B, d) - 1)
+            ind += stride * indexoffset(I[d])
             gap_lst_{d+1} *= stride
         end
 
@@ -874,3 +906,6 @@ If `dim` is specified, returns unique regions of the array `itr` along `dim`.
         @nref $N A d->d == dim ? sort!(uniquerows) : (indices(A, d))
     end
 end
+
+indexoffset(i) = first(i)-1
+indexoffset(::Colon) = 0
diff --git a/base/operators.jl b/base/operators.jl
index fd997396d3fd4..3f3a8fc95977c 100644
--- a/base/operators.jl
+++ b/base/operators.jl
@@ -2,6 +2,10 @@
 
 ## types ##
 
+typealias Dims{N} NTuple{N,Int}
+typealias DimsInteger{N} NTuple{N,Integer}
+typealias Indices{N} NTuple{N,AbstractUnitRange}
+
 const (<:) = issubtype
 
 supertype(T::DataType) = T.super
@@ -332,6 +336,8 @@ eltype(x) = eltype(typeof(x))
 
 # array shape rules
 
+promote_shape(::Tuple{}, ::Tuple{}) = ()
+
 function promote_shape(a::Tuple{Int,}, b::Tuple{Int,})
     if a[1] != b[1]
         throw(DimensionMismatch("dimensions must match"))
@@ -373,20 +379,24 @@ function promote_shape(a::Dims, b::Dims)
 end
 
 function promote_shape(a::AbstractArray, b::AbstractArray)
-    if ndims(a) < ndims(b)
+    promote_shape(indices(a), indices(b))
+end
+
+function promote_shape(a::Indices, b::Indices)
+    if length(a) < length(b)
         return promote_shape(b, a)
     end
-    for i=1:ndims(b)
-        if indices(a, i) != indices(b, i)
+    for i=1:length(b)
+        if a[i] != b[i]
             throw(DimensionMismatch("dimensions must match"))
         end
     end
-    for i=ndims(b)+1:ndims(a)
-        if indices(a, i) != 1:1
+    for i=length(b)+1:length(a)
+        if a[i] != 1:1
             throw(DimensionMismatch("dimensions must match"))
         end
     end
-    return indices(a)
+    return a
 end
 
 function throw_setindex_mismatch(X, I)
@@ -407,12 +417,12 @@ function setindex_shape_check(X::AbstractArray, I::Integer...)
     lj = length(I)
     i = j = 1
     while true
-        ii = size(X,i)
+        ii = length(indices(X,i))
         jj = I[j]
         if i == li || j == lj
             while i < li
                 i += 1
-                ii *= size(X,i)
+                ii *= length(indices(X,i))
             end
             while j < lj
                 j += 1
@@ -437,22 +447,22 @@ function setindex_shape_check(X::AbstractArray, I::Integer...)
 end
 
 setindex_shape_check(X::AbstractArray) =
-    (length(X)==1 || throw_setindex_mismatch(X,()))
+    (_length(X)==1 || throw_setindex_mismatch(X,()))
 
 setindex_shape_check(X::AbstractArray, i::Integer) =
-    (length(X)==i || throw_setindex_mismatch(X, (i,)))
+    (_length(X)==i || throw_setindex_mismatch(X, (i,)))
 
 setindex_shape_check{T}(X::AbstractArray{T,1}, i::Integer) =
-    (length(X)==i || throw_setindex_mismatch(X, (i,)))
+    (_length(X)==i || throw_setindex_mismatch(X, (i,)))
 
 setindex_shape_check{T}(X::AbstractArray{T,1}, i::Integer, j::Integer) =
-    (length(X)==i*j || throw_setindex_mismatch(X, (i,j)))
+    (_length(X)==i*j || throw_setindex_mismatch(X, (i,j)))
 
 function setindex_shape_check{T}(X::AbstractArray{T,2}, i::Integer, j::Integer)
     if length(X) != i*j
         throw_setindex_mismatch(X, (i,j))
     end
-    sx1 = size(X,1)
+    sx1 = length(indices(X,1))
     if !(i == 1 || i == sx1 || sx1 == 1)
         throw_setindex_mismatch(X, (i,j))
     end
diff --git a/base/permuteddimsarray.jl b/base/permuteddimsarray.jl
index 8ab7d60f26301..aa2944f7c4bb0 100644
--- a/base/permuteddimsarray.jl
+++ b/base/permuteddimsarray.jl
@@ -7,13 +7,12 @@ export permutedims
 # Some day we will want storage-order-aware iteration, so put perm in the parameters
 immutable PermutedDimsArray{T,N,perm,iperm,AA<:AbstractArray} <: AbstractArray{T,N}
     parent::AA
-    dims::NTuple{N,Int}
 
     function PermutedDimsArray(data::AA)
         (isa(perm, NTuple{N,Int}) && isa(iperm, NTuple{N,Int})) || error("perm and iperm must both be NTuple{$N,Int}")
         isperm(perm) || throw(ArgumentError(string(perm, " is not a valid permutation of dimensions 1:", N)))
         all(map(d->iperm[perm[d]]==d, 1:N)) || throw(ArgumentError(string(perm, " and ", iperm, " must be inverses")))
-        new(data, genperm(size(data), perm))
+        new(data)
     end
 end
 
@@ -24,7 +23,7 @@ function PermutedDimsArray{T,N}(data::AbstractArray{T,N}, perm)
 end
 
 Base.parent(A::PermutedDimsArray) = A.parent
-Base.size(A::PermutedDimsArray) = A.dims
+Base.size{T,N,perm}(A::PermutedDimsArray{T,N,perm})    = genperm(size(parent(A)),    perm)
 Base.indices{T,N,perm}(A::PermutedDimsArray{T,N,perm}) = genperm(indices(parent(A)), perm)
 
 @inline function Base.getindex{T,N,perm,iperm}(A::PermutedDimsArray{T,N,perm,iperm}, I::Vararg{Int,N})
diff --git a/base/range.jl b/base/range.jl
index 48a591438ebd3..62997616404a4 100644
--- a/base/range.jl
+++ b/base/range.jl
@@ -2,9 +2,6 @@
 
 ## 1-dimensional ranges ##
 
-typealias Dims{N} NTuple{N,Int}
-typealias DimsInteger{N} NTuple{N,Integer}
-
 abstract Range{T} <: AbstractArray{T,1}
 
 ## ordinal ranges
@@ -715,6 +712,7 @@ convert{T<:Real}(::Type{OneTo{T}}, r::OneTo) = OneTo{T}(r.stop)
 promote_rule{T1,UR<:AbstractUnitRange}(::Type{UnitRange{T1}}, ::Type{UR}) =
     UnitRange{promote_type(T1,eltype(UR))}
 convert{T<:Real}(::Type{UnitRange{T}}, r::AbstractUnitRange) = UnitRange{T}(first(r), last(r))
+convert(::Type{UnitRange}, r::AbstractUnitRange) = UnitRange(first(r), last(r))
 
 promote_rule{T1a,T1b,T2a,T2b}(::Type{StepRange{T1a,T1b}},::Type{StepRange{T2a,T2b}}) =
     StepRange{promote_type(T1a,T2a),promote_type(T1b,T2b)}
diff --git a/base/reduce.jl b/base/reduce.jl
index 679fa4b9bb7d8..4d26a67eee985 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -93,7 +93,7 @@ foldr(op, itr) = mapfoldr(identity, op, itr)
 
 ## reduce & mapreduce
 
-function mapreduce_impl(f, op, A::AbstractArray, ifirst::Int, ilast::Int, blksize::Int=pairwise_blocksize(f, op))
+function mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ilast::Integer, blksize::Int=pairwise_blocksize(f, op))
     if ifirst + blksize > ilast
         # sequential portion
         fx1 = r_promote(op, f(A[ifirst]))
@@ -141,23 +141,26 @@ mr_empty(f, op::typeof(|), T) = false
 _mapreduce(f, op, A::AbstractArray) = _mapreduce(f, op, linearindexing(A), A)
 
 function _mapreduce{T}(f, op, ::LinearFast, A::AbstractArray{T})
-    n = Int(length(A))
-    if n == 0
-        return mr_empty(f, op, T)
-    elseif n == 1
-        return r_promote(op, f(A[1]))
-    elseif n < 16
-        fx1 = r_promote(op, f(A[1]))
-        fx2 = r_promote(op, f(A[2]))
-        s = op(fx1, fx2)
-        i = 2
-        while i < n
-            @inbounds Ai = A[i+=1]
-            s = op(s, f(Ai))
+    inds = linearindices(A)
+    n = length(inds)
+    @inbounds begin
+        if n == 0
+            return mr_empty(f, op, T)
+        elseif n == 1
+            return r_promote(op, f(A[inds[1]]))
+        elseif n < 16
+            fx1 = r_promote(op, f(A[inds[1]]))
+            fx2 = r_promote(op, f(A[inds[2]]))
+            s = op(fx1, fx2)
+            i = inds[2]
+            while i < last(inds)
+                Ai = A[i+=1]
+                s = op(s, f(Ai))
+            end
+            return s
+        else
+            return mapreduce_impl(f, op, A, first(inds), last(inds))
         end
-        return s
-    else
-        return mapreduce_impl(f, op, A, 1, n)
     end
 end
 
diff --git a/base/reducedim.jl b/base/reducedim.jl
index 7f82a176cac98..f6fa07169a511 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -159,15 +159,13 @@ function check_reducedims(R, A)
     #   it will be size(A, 1) or size(A, 1) * size(A, 2).
     # - Otherwise, e.g. sum(A, 2) or sum(A, (1, 3)), it returns 0.
     #
-    ndims(R) <= ndims(A) || throw(DimensionMismatch("Cannot reduce $(ndims(A))-dimensional array to $(ndims(R)) dimensions"))
+    ndims(R) <= ndims(A) || throw(DimensionMismatch("cannot reduce $(ndims(A))-dimensional array to $(ndims(R)) dimensions"))
     lsiz = 1
     had_nonreduc = false
     for i = 1:ndims(A)
-        sRi = size(R, i)
-        sAi = size(A, i)
+        Ri, Ai = indices(R, i), indices(A, i)
+        sRi, sAi = length(Ri), length(Ai)
         if sRi == 1
-            first(indices(R, i)) == first(indices(A, i)) ||
-                throw(DimensionMismatch("Reduction along dimension $i must use lower indices"))
             if sAi > 1
                 if had_nonreduc
                     lsiz = 0  # to reduce along i, but some previous dimensions were non-reducing
@@ -176,8 +174,7 @@ function check_reducedims(R, A)
                 end
             end
         else
-            indices(R, i) == indices(A, i) ||
-                throw(DimensionMismatch("Reduction on array with indices $(indices(A)) with output with indices  $(indices(R))"))
+            Ri == Ai || throw(DimensionMismatch("reduction on array with indices $(indices(A)) with output with indices $(indices(R))"))
             had_nonreduc = true
         end
     end
@@ -187,11 +184,10 @@ end
 function _mapreducedim!{T,N}(f, op, R::AbstractArray, A::AbstractArray{T,N})
     lsiz = check_reducedims(R,A)
     isempty(A) && return R
-    sizA1 = size(A, 1)
 
     if has_fast_linear_indexing(A) && lsiz > 16
         # use mapreduce_impl, which is probably better tuned to achieve higher performance
-        nslices = div(length(A), lsiz)
+        nslices = div(_length(A), lsiz)
         ibase = first(linearindices(A))-1
         for i = 1:nslices
             @inbounds R[i] = op(R[i], mapreduce_impl(f, op, A, ibase+1, ibase+lsiz))
@@ -199,12 +195,13 @@ function _mapreducedim!{T,N}(f, op, R::AbstractArray, A::AbstractArray{T,N})
         end
         return R
     end
-    IRmax = dims_tail(map(last, indices(R)), A)
-    if size(R, 1) == 1 && sizA1 > 1
+    indsAt, indsRt = safe_tail(indices(A)), safe_tail(indices(R)) # handle d=1 manually
+    keep, Idefault = Broadcast.newindexer(indsAt, indsRt)
+    if reducedim1(R, A)
         # keep the accumulator as a local variable when reducing along the first dimension
-        i1 = first(indices(A, 1))
-        @inbounds for IA in CartesianRange(tail(indices(A)))
-            IR = min(IA, IRmax)
+        i1 = first(indices1(R))
+        @inbounds for IA in CartesianRange(indsAt)
+            IR = Broadcast.newindex(IA, keep, Idefault)
             r = R[i1,IR]
             @simd for i in indices(A, 1)
                 r = op(r, f(A[i, IA]))
@@ -212,8 +209,8 @@ function _mapreducedim!{T,N}(f, op, R::AbstractArray, A::AbstractArray{T,N})
             R[i1,IR] = r
         end
     else
-        @inbounds for IA in CartesianRange(tail(indices(A)))
-            IR = min(IA, IRmax)
+        @inbounds for IA in CartesianRange(indsAt)
+            IR = Broadcast.newindex(IA, keep, Idefault)
             @simd for i in indices(A, 1)
                 R[i,IR] = op(R[i,IR], f(A[i,IA]))
             end
@@ -273,18 +270,19 @@ end
 
 function findminmax!{T,N}(f, Rval, Rind, A::AbstractArray{T,N})
     (isempty(Rval) || isempty(A)) && return Rval, Rind
-    check_reducedims(Rval, A)
+    lsiz = check_reducedims(Rval, A)
     for i = 1:N
         indices(Rval, i) == indices(Rind, i) || throw(DimensionMismatch("Find-reduction: outputs must have the same indices"))
     end
     # If we're reducing along dimension 1, for efficiency we can make use of a temporary.
     # Otherwise, keep the result in Rval/Rind so that we traverse A in storage order.
-    IRmax = dims_tail(map(last, indices(Rval)), A)
+    indsAt, indsRt = safe_tail(indices(A)), safe_tail(indices(Rval))
+    keep, Idefault = Broadcast.newindexer(indsAt, indsRt)
     k = 0
-    if size(Rval, 1) < size(A, 1)
-        i1 = first(indices(A, 1))
-        @inbounds for IA in CartesianRange(tail(indices(A)))
-            IR = min(IRmax, IA)
+    if reducedim1(Rval, A)
+        i1 = first(indices1(Rval))
+        @inbounds for IA in CartesianRange(indsAt)
+            IR = Broadcast.newindex(IA, keep, Idefault)
             tmpRv = Rval[i1,IR]
             tmpRi = Rind[i1,IR]
             for i in indices(A,1)
@@ -299,8 +297,8 @@ function findminmax!{T,N}(f, Rval, Rind, A::AbstractArray{T,N})
             Rind[i1,IR] = tmpRi
         end
     else
-        @inbounds for IA in CartesianRange(tail(indices(A)))
-            IR = min(IRmax, IA)
+        @inbounds for IA in CartesianRange(indsAt)
+            IR = Broadcast.newindex(IA, keep, Idefault)
             for i in indices(A, 1)
                 k += 1
                 tmpAv = A[i,IA]
@@ -359,6 +357,4 @@ function findmax{T}(A::AbstractArray{T}, region)
             zeros(Int, reduced_dims0(A, region)), A)
 end
 
-dims_tail{T}(dims::Tuple{}, Aref::AbstractArray{T,0}) = CartesianIndex(())
-dims_tail{T,N}(dims::NTuple{N,Int}, Aref::AbstractArray{T,N}) = CartesianIndex(tail(dims))
-@inline dims_tail{T,M,N}(dims::NTuple{M,Int}, Aref::AbstractArray{T,N}) = dims_tail(tuple(dims..., 1), Aref)
+reducedim1(R, A) = length(indices1(R)) == 1
diff --git a/base/reshapedarray.jl b/base/reshapedarray.jl
index 053a554eb821a..d1a4b0961a5f6 100644
--- a/base/reshapedarray.jl
+++ b/base/reshapedarray.jl
@@ -52,7 +52,8 @@ end
 @pure rdims{N}(out::Tuple, inds::Tuple{Any, Vararg{Any}}, ::Type{Val{N}}) = rdims((out..., first(inds)), tail(inds), Val{N})
 
 function _reshape(parent::AbstractArray, dims::Dims)
-    prod(dims) == length(parent) || throw(DimensionMismatch("parent has $(length(parent)) elements, which is incompatible with size $dims"))
+    n = _length(parent)
+    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
     __reshape((parent, linearindexing(parent)), dims)
 end
 _reshape(R::ReshapedArray, dims::Dims) = _reshape(R.parent, dims)
diff --git a/base/serialize.jl b/base/serialize.jl
index 30c38312a1f2b..f0e6d938fba49 100644
--- a/base/serialize.jl
+++ b/base/serialize.jl
@@ -227,8 +227,8 @@ end
 trimmedsize(V) = index_lengths(V.parent, V.indexes...)
 
 function _trimmedsubarray{T,N,P,I,LD}(A, V::SubArray{T,N,P,I,LD}, newindexes)
-    LD && return SubArray{T,N,P,I,LD}(A, newindexes, size(V), Base.compute_offset1(A, 1, newindexes), 1)
-    SubArray{T,N,P,I,LD}(A, newindexes, size(V), 0, 0)
+    LD && return SubArray{T,N,P,I,LD}(A, newindexes, Base.compute_offset1(A, 1, newindexes), 1)
+    SubArray{T,N,P,I,LD}(A, newindexes, 0, 0)
 end
 _trimmedsubarray(A, V, newindexes, index::ViewIndex, indexes...) = _trimmedsubarray(A, V, (newindexes..., trimmedindex(V.parent, length(newindexes)+1, index)), indexes...)
 
diff --git a/base/show.jl b/base/show.jl
index fe9db2c308e44..51b04e057ae47 100644
--- a/base/show.jl
+++ b/base/show.jl
@@ -1286,7 +1286,7 @@ function alignment(io::IO, X::AbstractVecOrMat,
             break
         end
     end
-    if 1 < length(a) < size(X,2)
+    if 1 < length(a) < length(indices(X,2))
         while sum(map(sum,a)) + sep*length(a) >= cols_otherwise
             pop!(a)
         end
@@ -1384,8 +1384,9 @@ function print_matrix(io::IO, X::AbstractVecOrMat,
     postsp = ""
     @assert strwidth(hdots) == strwidth(ddots)
     sepsize = length(sep)
-    m, n = size(X,1), size(X,2)
-    rowsA, colsA = collect(indices(X,1)), collect(indices(X,2))
+    inds1, inds2 = indices(X,1), indices(X,2)
+    m, n = length(inds1), length(inds2)
+    rowsA, colsA = collect(inds1), collect(inds2)
     # To figure out alignments, only need to look at as many rows as could
     # fit down screen. If screen has at least as many rows as A, look at A.
     # If not, then we only need to look at the first and last chunks of A,
@@ -1488,19 +1489,19 @@ function show_nd(io::IO, a::AbstractArray, print_matrix, label_slices)
     if isempty(a)
         return
     end
-    tail = indices(a)[3:end]
+    tailinds = tail(tail(indices(a)))
     nd = ndims(a)-2
-    for I in CartesianRange(tail)
+    for I in CartesianRange(tailinds)
         idxs = I.I
         if limit
             for i = 1:nd
                 ii = idxs[i]
-                ind = tail[i]
+                ind = tailinds[i]
                 if length(ind) > 10
-                    if ii == ind[4] && all(d->idxs[d]==first(tail[d]),1:i-1)
+                    if ii == ind[4] && all(d->idxs[d]==first(tailinds[d]),1:i-1)
                         for j=i+1:nd
                             szj = size(a,j+2)
-                            indj = tail[j]
+                            indj = tailinds[j]
                             if szj>10 && first(indj)+2 < idxs[j] <= last(indj)-3
                                 @goto skip
                             end
@@ -1522,7 +1523,7 @@ function show_nd(io::IO, a::AbstractArray, print_matrix, label_slices)
         end
         slice = view(a, indices(a,1), indices(a,2), idxs...)
         print_matrix(io, slice)
-        print(io, idxs == map(last,tail) ? "" : "\n\n")
+        print(io, idxs == map(last,tailinds) ? "" : "\n\n")
         @label skip
     end
 end
@@ -1536,10 +1537,11 @@ function print_matrix_repr(io, X::AbstractArray)
     if compact && !haskey(io, :compact)
         io = IOContext(io, :compact => compact)
     end
-    nr, nc = size(X,1), size(X,2)
+    indr, indc = indices(X,1), indices(X,2)
+    nr, nc = length(indr), length(indc)
     rdots, cdots = false, false
-    rr1, rr2 = indices(X,1), 1:0
-    cr1, cr2 = indices(X,2), 1:0
+    rr1, rr2 = UnitRange{Int}(indr), 1:0
+    cr1, cr2 = UnitRange{Int}(indc), 1:0
     if limit
         if nr > 4
             rr1, rr2 = rr1[1:2], rr1[nr-1:nr]
@@ -1563,8 +1565,8 @@ function print_matrix_repr(io, X::AbstractArray)
                         show(io, el)
                     end
                 end
-                if last(cr) == last(indices(X,2))
-                    i < last(indices(X,1)) && print(io, "; ")
+                if last(cr) == last(indc)
+                    i < last(indr) && print(io, "; ")
                 elseif cdots
                     print(io, " \u2026 ")
                 end
diff --git a/base/sort.jl b/base/sort.jl
index 99900be146ef4..d388e55f6647c 100644
--- a/base/sort.jl
+++ b/base/sort.jl
@@ -482,16 +482,16 @@ function sort(A::AbstractArray, dim::Integer;
               order::Ordering=Forward,
               initialized::Bool=false)
     order = ord(lt,by,rev,order)
+    n = length(indices(A, dim))
     if dim != 1
         pdims = (dim, setdiff(1:ndims(A), dim)...)  # put the selected dimension first
-        Ap = permutedims(A, pdims)    # note Ap is an Array, no matter what A is
-        n = size(Ap, 1)
+        Ap = permutedims(A, pdims)
         Av = vec(Ap)
         sort_chunks!(Av, n, alg, order)
         ipermutedims(Ap, pdims)
     else
         Av = A[:]
-        sort_chunks!(Av, size(A,1), alg, order)
+        sort_chunks!(Av, n, alg, order)
         reshape(Av, indices(A))
     end
 end
diff --git a/base/statistics.jl b/base/statistics.jl
index 764a8fa70a394..05d9b6fe3b28e 100644
--- a/base/statistics.jl
+++ b/base/statistics.jl
@@ -96,10 +96,9 @@ function centralize_sumabs2!{S,T,N}(R::AbstractArray{S}, A::AbstractArray{T,N},
     lsiz = check_reducedims(R,A)
     isempty(R) || fill!(R, zero(S))
     isempty(A) && return R
-    sizA1 = size(A, 1)
 
     if has_fast_linear_indexing(A) && lsiz > 16
-        nslices = div(length(A), lsiz)
+        nslices = div(_length(A), lsiz)
         ibase = first(linearindices(A))-1
         for i = 1:nslices
             @inbounds R[i] = centralize_sumabs2(A, means[i], ibase+1, ibase+lsiz)
@@ -107,11 +106,12 @@ function centralize_sumabs2!{S,T,N}(R::AbstractArray{S}, A::AbstractArray{T,N},
         end
         return R
     end
-    IRmax = dims_tail(map(last, indices(R)), A)
-    if size(R, 1) == 1 && sizA1 > 1
-        i1 = first(indices(A, 1))
-        @inbounds for IA in CartesianRange(tail(indices(A)))
-            IR = min(IA, IRmax)
+    indsAt, indsRt = safe_tail(indices(A)), safe_tail(indices(R)) # handle d=1 manually
+    keep, Idefault = Broadcast.newindexer(indsAt, indsRt)
+    if reducedim1(R, A)
+        i1 = first(indices1(R))
+        @inbounds for IA in CartesianRange(indsAt)
+            IR = Broadcast.newindex(IA, keep, Idefault)
             r = R[i1,IR]
             m = means[i1,IR]
             @simd for i in indices(A, 1)
@@ -120,8 +120,8 @@ function centralize_sumabs2!{S,T,N}(R::AbstractArray{S}, A::AbstractArray{T,N},
             R[i1,IR] = r
         end
     else
-        @inbounds for IA in CartesianRange(tail(indices(A)))
-            IR = min(IA, IRmax)
+        @inbounds for IA in CartesianRange(indsAt)
+            IR = Broadcast.newindex(IA, keep, Idefault)
             @simd for i in indices(A, 1)
                 R[i,IR] += abs2(A[i,IA] - means[i,IR])
             end
diff --git a/base/subarray.jl b/base/subarray.jl
index 48446fcc60cb4..390b6fc5aa07e 100644
--- a/base/subarray.jl
+++ b/base/subarray.jl
@@ -8,25 +8,24 @@ abstract AbstractCartesianIndex{N} # This is a hacky forward declaration for Car
 immutable SubArray{T,N,P,I,L} <: AbstractArray{T,N}
     parent::P
     indexes::I
-    dims::NTuple{N,Int}
     offset1::Int       # for linear indexing and pointer, only valid when L==true
     stride1::Int       # used only for linear indexing
-    function SubArray(parent, indexes, dims, offset1, stride1)
+    function SubArray(parent, indexes, offset1, stride1)
         check_parent_index_match(parent, indexes)
-        new(parent, indexes, dims, offset1, stride1)
+        new(parent, indexes, offset1, stride1)
     end
 end
 # Compute the linear indexability of the indices, and combine it with the linear indexing of the parent
 function SubArray(parent::AbstractArray, indexes::Tuple, dims::Tuple)
-    SubArray(linearindexing(viewindexing(indexes), linearindexing(parent)), parent, indexes, map(dimlength, dims))
+    SubArray(linearindexing(viewindexing(indexes), linearindexing(parent)), parent, indexes, dims)
 end
-function SubArray{P, I, N}(::LinearSlow, parent::P, indexes::I, dims::NTuple{N, Int})
-    SubArray{eltype(P), N, P, I, false}(parent, indexes, dims, 0, 0)
+function SubArray{P, I, N}(::LinearSlow, parent::P, indexes::I, dims::NTuple{N})
+    SubArray{eltype(P), N, P, I, false}(parent, indexes, 0, 0)
 end
-function SubArray{P, I, N}(::LinearFast, parent::P, indexes::I, dims::NTuple{N, Int})
+function SubArray{P, I, N}(::LinearFast, parent::P, indexes::I, dims::NTuple{N})
     # Compute the stride and offset
     stride1 = compute_stride1(parent, indexes)
-    SubArray{eltype(P), N, P, I, true}(parent, indexes, dims, compute_offset1(parent, stride1, indexes), stride1)
+    SubArray{eltype(P), N, P, I, true}(parent, indexes, compute_offset1(parent, stride1, indexes), stride1)
 end
 
 check_parent_index_match{T,N}(parent::AbstractArray{T,N}, indexes::NTuple{N}) = nothing
@@ -47,12 +46,8 @@ viewindexing(I::Tuple{Vararg{Any}}) = LinearSlow()
 # Of course, all other array types are slow
 viewindexing(I::Tuple{AbstractArray, Vararg{Any}}) = LinearSlow()
 
-dimlength(r::Range) = length(r)
-dimlength(i::Integer) = Int(i)
-
 # Simple utilities
-size(V::SubArray) = V.dims
-length(V::SubArray) = prod(V.dims)
+size(V::SubArray) = (@_inline_meta; map(n->Int(unsafe_length(n)), indices(V)))
 
 similar(V::SubArray, T::Type, dims::Dims) = similar(V.parent, T, dims)
 
@@ -83,12 +78,12 @@ end
 function unsafe_view{T,N}(A::AbstractArray{T,N}, I::Vararg{ViewIndex,N})
     @_inline_meta
     J = to_indexes(I...)
-    SubArray(A, J, map(dimlength, index_shape(A, J...)))
+    SubArray(A, J, map(unsafe_length, index_shape(A, J...)))
 end
 function unsafe_view{T,N}(V::SubArray{T,N}, I::Vararg{ViewIndex,N})
     @_inline_meta
     idxs = reindex(V, V.indexes, to_indexes(I...))
-    SubArray(V.parent, idxs, map(dimlength, (index_shape(V.parent, idxs...))))
+    SubArray(V.parent, idxs, map(unsafe_length, (index_shape(V.parent, idxs...))))
 end
 
 # Re-indexing is the heart of a view, transforming A[i, j][x, y] to A[i[x], j[y]]
@@ -199,34 +194,25 @@ substrides(s, parent, dim, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("st
 
 stride(V::SubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]
 
-compute_stride1(parent, I::Tuple) = compute_stride1(1, parent, 1, I)
-compute_stride1(s, parent, dim, I::Tuple{}) = s
-compute_stride1(s, parent, dim, I::Tuple{DroppedScalar, Vararg{Any}}) =
-    (@_inline_meta; compute_stride1(s*size(parent, dim), parent, dim+1, tail(I)))
-compute_stride1(s, parent, dim, I::Tuple{Range, Vararg{Any}}) = s*step(I[1])
-compute_stride1(s, parent, dim, I::Tuple{Colon, Vararg{Any}}) = s
-compute_stride1(s, parent, dim, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("invalid strided index type $(typeof(I[1]))"))
+compute_stride1{N}(parent::AbstractArray, I::NTuple{N}) =
+    compute_stride1(1, fill_to_length(indices(parent), OneTo(1), Val{N}), I)
+compute_stride1(s, inds, I::Tuple{}) = s
+compute_stride1(s, inds, I::Tuple{DroppedScalar, Vararg{Any}}) =
+    (@_inline_meta; compute_stride1(s*unsafe_length(inds[1]), tail(inds), tail(I)))
+compute_stride1(s, inds, I::Tuple{Range, Vararg{Any}}) = s*step(I[1])
+compute_stride1(s, inds, I::Tuple{Colon, Vararg{Any}}) = s
+compute_stride1(s, inds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("invalid strided index type $(typeof(I[1]))"))
 
 iscontiguous(A::SubArray) = iscontiguous(typeof(A))
 iscontiguous{S<:SubArray}(::Type{S}) = false
 iscontiguous{F<:FastContiguousSubArray}(::Type{F}) = true
 
-# Fast linear SubArrays have their first index cached
-first_index(V::FastSubArray) = V.offset1 + V.stride1
-first_index(V::SubArray) = first_index(V.parent, V.indexes)
-function first_index(P::AbstractArray, indexes::Tuple)
-    f = first(linearindices(P))
-    s = 1
-    for i = 1:length(indexes)
-        f += (_first(indexes[i], P, i)-first(indices(P, i)))*s
-        s *= size(P, i)
-    end
-    f
+first_index(V::FastSubArray) = V.offset1 + V.stride1 # cached for fast linear SubArrays
+function first_index(V::SubArray)
+    P, I = parent(V), V.indexes
+    s1 = compute_stride1(P, I)
+    s1 + compute_offset1(P, s1, I)
 end
-_first(::Colon, P, ::Colon) = first(linearindices(P))
-_first(i, P, ::Colon) = first(i)
-_first(::Colon, P, d) = first(indices(P, d))
-_first(i, P, d) = first(i)
 
 # Computing the first index simply steps through the indices, accumulating the
 # sum of index each multiplied by the parent's stride.
@@ -238,17 +224,26 @@ compute_offset1(parent, stride1::Integer, I::Tuple) = (@_inline_meta; compute_of
 compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{Colon}, I::Tuple) = compute_linindex(parent, I) - stride1*first(indices(parent, dims[1]))  # index-preserving case
 compute_offset1(parent, stride1::Integer, dims, inds, I::Tuple) = compute_linindex(parent, I) - stride1  # linear indexing starts with 1
 
-compute_linindex(parent, I) = compute_linindex(1, 1, parent, 1, I)
-compute_linindex(f, s, parent, dim, I::Tuple{Real, Vararg{Any}}) =
-    (@_inline_meta; compute_linindex(f + (I[1]-first(indices(parent,dim)))*s, s*size(parent, dim), parent, dim+1, tail(I)))
+function compute_linindex{N}(parent, I::NTuple{N})
+    IP = fill_to_length(indices(parent), OneTo(1), Val{N})
+    compute_linindex(1, 1, IP, I)
+end
+function compute_linindex(f, s, IP::Tuple, I::Tuple{Real, Vararg{Any}})
+    @_inline_meta
+    Δi = I[1]-first(IP[1])
+    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
+end
 # Just splat out the cartesian indices and continue
-compute_linindex(f, s, parent, dim, I::Tuple{AbstractCartesianIndex, Vararg{Any}}) =
-    (@_inline_meta; compute_linindex(f, s, parent, dim, (I[1].I..., tail(I)...)))
-compute_linindex(f, s, parent, dim, I::Tuple{Colon, Vararg{Any}}) =
-    (@_inline_meta; compute_linindex(f, s*size(parent, dim), parent, dim+1, tail(I)))
-compute_linindex(f, s, parent, dim, I::Tuple{Any, Vararg{Any}}) =
-    (@_inline_meta; compute_linindex(f + (first(I[1])-first(indices(parent,dim)))*s, s*size(parent, dim), parent, dim+1, tail(I)))
-compute_linindex(f, s, parent, dim, I::Tuple{}) = f
+compute_linindex(f, s, IP::Tuple, I::Tuple{AbstractCartesianIndex, Vararg{Any}}) =
+    (@_inline_meta; compute_linindex(f, s, IP, (I[1].I..., tail(I)...)))
+compute_linindex(f, s, IP::Tuple, I::Tuple{Colon, Vararg{Any}}) =
+    (@_inline_meta; compute_linindex(f, s*unsafe_length(IP[1]), tail(IP), tail(I)))
+function compute_linindex(f, s, IP::Tuple, I::Tuple{Any, Vararg{Any}})
+    @_inline_meta
+    Δi = first(I[1])-first(IP[1])
+    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
+end
+compute_linindex(f, s, IP::Tuple, I::Tuple{}) = f
 
 find_extended_dims(I) = (@_inline_meta; _find_extended_dims((), (), 1, I...))
 _find_extended_dims(dims, inds, dim) = dims, inds
diff --git a/base/tuple.jl b/base/tuple.jl
index 08c245f837977..eabbce23e4f18 100644
--- a/base/tuple.jl
+++ b/base/tuple.jl
@@ -43,6 +43,10 @@ first(t::Tuple) = t[1]
 eltype(::Type{Tuple{}}) = Bottom
 eltype{T}(::Type{Tuple{Vararg{T}}}) = T
 
+# version of tail that doesn't throw on empty tuples (used in array indexing)
+safe_tail(t::Tuple) = tail(t)
+safe_tail(t::Tuple{}) = ()
+
 # front (the converse of tail: it skips the last entry)
 
 function front(t::Tuple)
diff --git a/doc/devdocs/julia.rst b/doc/devdocs/julia.rst
index 32c1c212caebb..141c9f4933b81 100644
--- a/doc/devdocs/julia.rst
+++ b/doc/devdocs/julia.rst
@@ -25,3 +25,4 @@
    promote-op
    boundscheck
    locks
+   offset-arrays
diff --git a/doc/devdocs/offset-arrays.rst b/doc/devdocs/offset-arrays.rst
new file mode 100644
index 0000000000000..5d6beb8563251
--- /dev/null
+++ b/doc/devdocs/offset-arrays.rst
@@ -0,0 +1,293 @@
+.. currentmodule:: Base
+
+.. _devdocs-offsetarrays:
+
+**************************
+Arrays with custom indices
+**************************
+
+Julia 0.5 adds experimental support for arrays with arbitrary indices.
+Conventionally, Julia's arrays are indexed starting at 1, whereas some
+other languages start numbering at 0, and yet others (e.g., Fortran)
+allow you to specify arbitrary starting indices.  While there is much
+merit in picking a standard (i.e., 1 for Julia), there are some
+algorithms which simplify considerably if you can index outside the
+range ``1:size(A,d)`` (and not just ``0:size(A,d)-1``, either). Such
+array types are expected to be supplied through packages.
+
+The purpose of this page is to address the question, "what do I
+have to do to support such arrays in my own code?"  First, let's
+address the simplest case: if you know that your code will never need
+to handle arrays with unconventional indexing, hopefully the answer is
+"nothing." Old code, on conventional arrays, should function
+essentially without alteration as long as it was using the exported
+interfaces of Julia.
+
+Generalizing existing code
+--------------------------
+
+As an overview, the steps are:
+
+- replace many uses of ``size`` with ``indices``
+- replace ``1:length(A)`` with ``linearindices(A)``, and ``length(A)`` with ``length(linearindices(A))``
+- replace explicit allocations like ``Array{Int}(size(B))`` with ``similar(Array{Int}, indices(B))``
+
+These are described in more detail below.
+
+Background
+~~~~~~~~~~
+
+Because unconventional indexing breaks deeply-held assumptions
+throughout the Julia ecosystem, early adopters running code that has
+not been updated are likely to experience errors.  The most
+frustrating bugs would be incorrect results or segfaults (total
+crashes of Julia).  For example, consider the following function::
+
+  function mycopy!(dest::AbstractVector, src::AbstractVector)
+      length(dest) == length(src) || throw(DimensionMismatch("vectors must match"))
+      # OK, now we're safe to use @inbounds, right? (not anymore!)
+      for i = 1:length(src)
+          @inbounds dest[i] = src[i]
+      end
+      dest
+  end
+
+This code implicitly assumes that vectors are indexed from 1.
+Previously that was a safe assumption, so this code was fine, but
+(depending on what types the user passes to this function) it may no
+longer be safe.  If this code continued to work when passed a vector
+with non-1 indices, it would either produce an incorrect answer or it
+would segfault.  (If you do get segfaults, to help locate the
+cause try running julia with the option ``--check-bounds=yes``.)
+
+To ensure that such errors are caught, in Julia 0.5 both ``length``
+and ``size`` **should** throw an error when passed an array with
+non-1 indexing.  This is designed to force users of such arrays to
+check the code, and inspect it for whether it needs to be generalized.
+
+Using ``indices`` for bounds checks and loop iteration
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``indices(A)`` (reminiscent of ``size(A)``) returns a tuple of
+``AbstractUnitRange`` objects, specifying the range of valid indices
+along each dimension of ``A``.  When ``A`` has unconventional
+indexing, the ranges may not start at 1.  If you just want the range
+for a particular dimension ``d``, there is ``indices(A, d)``.
+
+Base implements a custom range type, ``OneTo``, where ``OneTo(n)``
+means the same thing as ``1:n`` but in a form that guarantees (via the
+type system) that the lower index is 1.  For any new ``AbstractArray``
+type, this is the default returned by ``indices``, and it indicates
+that this array type uses "conventional" 1-based indexing.  Note that
+if you don't want to be bothered supporting arrays with non-1
+indexing, you can add the following line::
+
+   @assert all(x->isa(x, Base.OneTo), indices(A))
+
+at the top of any function.
+
+For bounds checking, note that there are dedicated functions
+``checkbounds`` and ``checkindex`` which can sometimes simplify such
+tests.
+
+Linear indexing (``linearindices``)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Some algorithms are most conveniently (or efficiently) written in
+terms of a single linear index, ``A[i]`` even if ``A`` is
+multi-dimensional.  In "true" linear indexing, the indices always
+range from ``1:length(A)``. However, this raises an ambiguity for
+one-dimensional arrays (a.k.a., ``AbstractVector``): does ``v[i]``
+mean linear indexing, or Cartesian indexing with the array's native
+indices?
+
+For this reason, if you want to use linear indexing in an algorithm,
+your best option is 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.
+
+In a sense, one can say that 1-dimensional arrays always use Cartesian
+indexing. 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.
+
+Using ``indices`` and ``linearindices``, here is one way you could
+rewrite ``mycopy!``::
+
+  function mycopy!(dest::AbstractVector, src::AbstractVector)
+      indices(dest) == indices(src) || throw(DimensionMismatch("vectors must match"))
+      for i in linearindices(src)
+          @inbounds dest[i] = src[i]
+      end
+      dest
+  end
+
+Allocating storage using generalizations of ``similar``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Storage is often allocated with ``Array{Int}(dims)`` or ``similar(A, args...)``.
+When the result needs to match the indices of some other array, this may not always suffice.
+The generic replacement for such patterns is to use
+``similar(storagetype, shape)``.  ``storagetype`` indicates the kind
+of underlying "conventional" behavior you'd like, e.g., ``Array{Int}``
+or ``BitArray`` or even ``dims->zeros(Float32, dims)`` (which would
+allocate an all-zeros array). ``shape`` is a tuple of ``Integer`` or
+``AbstractUnitRange`` values, specifying the indices
+that you want the result to use.
+
+Let's walk through a couple of explicit examples. First, if ``A`` has
+conventional indices, then ``similar(Array{Int}, indices(A))``
+would end up calling ``Array{Int}(size(A))``, and thus return an
+array.  If ``A`` is an ``AbstractArray`` type with unconventional
+indexing, then ``similar(Array{Int}, indices(A))`` should return
+something that "behaves like" an ``Array{Int}`` but with a shape
+(including indices) that matches ``A``.  (The most obvious
+implementation is to allocate an ``Array{Int}(size(A))`` and then
+"wrap" it in a type that shifts the indices.)
+
+Note also that ``similar(Array{Int}, (indices(A, 2),))`` would allocate
+an ``AbstractVector{Int}`` (i.e., 1-dimensional array) that matches the
+indices of the columns of ``A``.
+
+
+Deprecations
+~~~~~~~~~~~~
+
+In generalizing Julia's code base, at least one deprecation was
+unavoidable: earlier versions of Julia defined ``first(::Colon) = 1``,
+meaning that the first index along a dimension indexed by ``:`` is 1.
+This definition can no longer be justified, so it was deprecated.
+There is no provided replacement, because the proper replacement
+depends on what you are doing and might need to know more about the
+array. However, it appears that many uses of ``first(::Colon)`` are
+really about computing an index offset; when that is the case, a
+candidate replacement is::
+
+  indexoffset(r::AbstractVector) = first(r) - 1
+  indexoffset(::Colon) = 0
+
+In other words, while ``first(:)`` does not itself make sense, in
+general you can say that the offset associated with a colon-index is
+zero.
+
+Writing custom array types with non-1 indexing
+----------------------------------------------
+
+Most of the methods you'll need to define are standard for any
+``AbstractArray`` type, see :ref:`man-interfaces-abstractarray`.
+This page focuses on the steps needed to define unconventional
+indexing.
+
+Do **not** implement ``size`` or ``length``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Perhaps the majority of pre-existing code that uses ``size`` will not
+work properly for arrays with non-1 indices.  For that reason, it is
+much better to avoid implementing these methods, and use the resulting
+``MethodError`` to identify code that needs to be audited and perhaps
+generalized.
+
+Do **not** annotate bounds checks
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Julia 0.5 includes ``@boundscheck`` to annotate code that can be
+removed for callers that exploit ``@inbounds``. Initially, it seems
+far preferable to run with bounds checking always enabled (i.e., omit
+the ``@boundscheck`` annotation so the check always runs).
+
+Custom ``AbstractUnitRange`` types
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+If you're writing a non-1 indexed array type, you will want to
+specialize ``indices`` so it returns a ``UnitRange``, or (perhaps
+better) a custom ``AbstractUnitRange``.  The advantage of a custom
+type is that it "signals" the allocation type for functions like
+``similar``. If we're writing an array type for which indexing will
+start at 0, we likely want to begin by creating a new
+``AbstractUnitRange``, ``ZeroRange``, where ``ZeroRange(n)`` is equivalent
+to ``0:n-1``.
+
+In general, you should probably *not* export ``ZeroRange`` from your
+package: there may be other packages that implement their own
+``ZeroRange``, and having multiple distinct ``ZeroRange`` types is (perhaps
+counterintuitively) an advantage: ``ModuleA.ZeroRange`` indicates that
+``similar`` should create a ``ModuleA.ZeroArray``, whereas
+``ModuleB.ZeroRange`` indicates a ``ModuleB.ZeroArray`` type.  This
+design allows peaceful coexistence among many different custom array
+types.
+
+Note that the Julia package ``CustomUnitRanges.jl`` can sometimes be used to
+avoid the need to write your own ``ZeroRange`` type.
+
+Specializing ``indices``
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Once you have your ``AbstractUnitRange`` type, then specialize ``indices``::
+
+  Base.indices(A::ZeroArray) = map(n->ZeroRange(n), A.size)
+
+where here we imagine that ``ZeroArray`` has a field called ``size``
+(there would be other ways to implement this).
+
+In some cases, the fallback definition for ``indices(A, d)``::
+
+   indices{T,N}(A::AbstractArray{T,N}, d) = d <= N ? indices(A)[d] : OneTo(1)
+
+may not be what you want: you may need to specialize it to return
+something other than ``OneTo(1)`` when ``d > ndims(A)``.  Likewise, in
+``Base`` there is a dedicated function ``indices1`` which is
+equivalent to ``indices(A, 1)`` but which avoids checking (at
+runtime) whether ``ndims(A) > 0``. (This is purely a performance
+optimization.)  It is defined as::
+
+  indices1{T}(A::AbstractArray{T,0}) = OneTo(1)
+  indices1{T}(A::AbstractArray{T})   = indices(A)[1]
+
+If the first of these (the zero-dimensional case) is problematic for
+your custom array type, be sure to specialize it appropriately.
+
+Specializing ``similar``
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Given your custom ``ZeroRange`` type, then you should also add the following two specializations for ``similar``::
+
+  function Base.similar(A::AbstractArray, T::Type, shape::Tuple{ZeroRange,Vararg{ZeroRange}})
+      # body
+  end
+
+  function Base.similar(f::Union{Function,DataType}, shape::Tuple{ZeroRange,Vararg{ZeroRange}})
+      # body
+  end
+
+Both of these should allocate your custom array type.
+
+Specializing ``reshape``
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Optionally, define a method
+::
+
+   Base.reshape(A::AbstractArray, shape::Tuple{ZeroRange,Vararg{ZeroRange}}) = ...
+
+and you can ``reshape`` an array so that the result has custom indices.
+
+Summary
+-------
+
+Writing code that doesn't make assumptions about indexing requires a
+few extra abstractions, but hopefully the necessary changes are
+relatively straightforward.
+
+As a reminder, this support is still experimental. While much of
+Julia's base code has been updated to support unconventional indexing,
+without a doubt there are many omissions that will be discovered only
+through usage.  Moreover, at the time of this writing, most packages
+do not support unconventional indexing.  As a consequence, early
+adopters should be prepared to identify and/or fix bugs.  On the other
+hand, only through practical usage will it become clear whether this
+experimental feature should be retained in future versions of Julia;
+consequently, interested parties are encouraged to accept some
+ownership for putting it through its paces.
diff --git a/doc/manual/interfaces.rst b/doc/manual/interfaces.rst
index b0f0886da6c22..15465047b2402 100644
--- a/doc/manual/interfaces.rst
+++ b/doc/manual/interfaces.rst
@@ -166,8 +166,9 @@ Methods to implement
 :func:`similar(A, dims::NTuple{Int}) <similar>`                       ``similar(A, eltype(A), dims)``              Return a mutable array with the same element type and size `dims`
 :func:`similar(A, ::Type{S}, dims::NTuple{Int}) <similar>`            ``Array{S}(dims)``                           Return a mutable array with the specified element type and size
 **Non-traditional indices**                                           **Default definition**                       **Brief description**
-:func:`indices(A, d) <indices>`                                       ``OneTo(size(A, d))``                        Return the ``AbstractUnitRange`` of valid indices along dimension ``d``
-:func:`Base.similar(A, ::Type{S}, inds::NTuple{Ind}) <similar>`       ``similar(A, S, map(Base.dimlength, inds))`` Return a mutable array with the specified indices ``inds`` (see below)
+:func:`indices(A) <indices>`                                          ``map(OneTo, size(A))``                      Return the ``AbstractUnitRange`` of valid indices
+:func:`Base.similar(A, ::Type{S}, inds::NTuple{Ind}) <similar>`       ``similar(A, S, Base.to_shape(inds))``       Return a mutable array with the specified indices ``inds`` (see below)
+:func:`Base.similar(T::Union{Type,Function}, inds) <similar>`         ``T(Base.to_shape(inds))``                   Return an array similar to ``T`` with the specified indices ``inds`` (see below)
 ===================================================================== ============================================ =======================================================================================
 
 If a type is defined as a subtype of ``AbstractArray``, it inherits a very large set of rich behaviors including iteration and multidimensional indexing built on top of single-element access.  See the :ref:`arrays manual page <man-arrays>` and :ref:`standard library section <stdlib-arrays>` for more supported methods.
@@ -292,6 +293,4 @@ If you are defining an array type that allows non-traditional indexing
 ``indices``.  You should also specialize ``similar`` so that the
 ``dims`` argument (ordinarily a ``Dims`` size-tuple) can accept
 ``AbstractUnitRange`` objects, perhaps range-types ``Ind`` of your own
-design.  For example, if indexing always starts with 0 for your
-arrays, you likely want to define a ``ZeroTo`` range type.  Otherwise,
-you can use standard ``UnitRange``.
+design.  For more information, see :ref:`devdocs-offsetarrays`.
diff --git a/test/arrayops.jl b/test/arrayops.jl
index 741fac3fb5816..00e4ca331153c 100644
--- a/test/arrayops.jl
+++ b/test/arrayops.jl
@@ -1218,27 +1218,37 @@ a = ones(5,0)
 b = view(a, :, :)
 @test mdsum(b) == 0
 
-a = copy(reshape(1:60, 3, 4, 5))
-@test a[CartesianIndex{3}(2,3,4)] == 44
-a[CartesianIndex{3}(2,3,3)] = -1
-@test a[CartesianIndex{3}(2,3,3)] == -1
-@test a[2,CartesianIndex{2}(3,4)] == 44
-a[1,CartesianIndex{2}(3,4)] = -2
-@test a[1,CartesianIndex{2}(3,4)] == -2
-@test a[CartesianIndex{1}(2),3,CartesianIndex{1}(4)] == 44
-a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] = -3
-@test a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] == -3
-@test a[:, :, CartesianIndex((1,))] == a[:,:,1]
-@test a[CartesianIndex((1,)), [1,2], :] == a[1,[1,2],:]
-@test a[CartesianIndex((2,)), 3:4, :] == a[2,3:4,:]
-
-a = view(zeros(3, 4, 5), :, :, :)
-a[CartesianIndex{3}(2,3,3)] = -1
-@test a[CartesianIndex{3}(2,3,3)] == -1
-a[1,CartesianIndex{2}(3,4)] = -2
-@test a[1,CartesianIndex{2}(3,4)] == -2
-a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] = -3
-@test a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] == -3
+for a in (copy(reshape(1:60, 3, 4, 5)),
+          view(copy(reshape(1:60, 3, 4, 5)), 1:3, :, :))
+    @test a[CartesianIndex{3}(2,3,4)] == 44
+    a[CartesianIndex{3}(2,3,3)] = -1
+    @test a[CartesianIndex{3}(2,3,3)] == -1
+    @test a[2,CartesianIndex{2}(3,4)] == 44
+    a[1,CartesianIndex{2}(3,4)] = -2
+    @test a[1,CartesianIndex{2}(3,4)] == -2
+    @test a[CartesianIndex{1}(2),3,CartesianIndex{1}(4)] == 44
+    a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] = -3
+    @test a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] == -3
+    @test a[:, :, CartesianIndex((1,))] == a[:,:,1]
+    @test a[CartesianIndex((1,)), [1,2], :] == a[1,[1,2],:]
+    @test a[CartesianIndex((2,)), 3:4, :] == a[2,3:4,:]
+    @test a[[CartesianIndex(1,3),CartesianIndex(2,4)],3:3] == reshape([a[1,3,3]; a[2,4,3]], 2, 1)
+    @test_throws BoundsError a[[CartesianIndex(1,5),CartesianIndex(2,4)],3:3]
+    @test_throws BoundsError a[1:4, [CartesianIndex(1,3),CartesianIndex(2,4)]]
+end
+
+for a in (view(zeros(3, 4, 5), :, :, :),
+          view(zeros(3, 4, 5), 1:3, :, :))
+    a[CartesianIndex{3}(2,3,3)] = -1
+    @test a[CartesianIndex{3}(2,3,3)] == -1
+    a[1,CartesianIndex{2}(3,4)] = -2
+    @test a[1,CartesianIndex{2}(3,4)] == -2
+    a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] = -3
+    @test a[CartesianIndex{1}(2),3,CartesianIndex{1}(3)] == -3
+    a[[CartesianIndex(1,3),CartesianIndex(2,4)],3:3] = -4
+    @test a[1,3,3] == -4
+    @test a[2,4,3] == -4
+end
 
 I1 = CartesianIndex((2,3,0))
 I2 = CartesianIndex((-1,5,2))
diff --git a/test/broadcast.jl b/test/broadcast.jl
index 2668e21e68e05..2d6adca165282 100644
--- a/test/broadcast.jl
+++ b/test/broadcast.jl
@@ -42,13 +42,14 @@ check_broadcast_shape((-1:1, 6:9), (1, 6:9))
 check_broadcast_shape((-1:1, 6:9), 1)
 
 ci(x) = CartesianIndex(x)
-@test @inferred(newindex(ci((2,2)), (true, true)))   == ci((2,2))
-@test @inferred(newindex(ci((2,2)), (true, false)))  == ci((2,1))
-@test @inferred(newindex(ci((2,2)), (false, true)))  == ci((1,2))
-@test @inferred(newindex(ci((2,2)), (false, false))) == ci((1,1))
-@test @inferred(newindex(ci((2,2)), (true,)))   == ci((2,))
-@test @inferred(newindex(ci((2,2)), (false,))) == ci((1,))
-@test @inferred(newindex(ci((2,2)), ())) == 1
+@test @inferred(newindex(ci((2,2)), (true, true), (-1,-1)))   == ci((2,2))
+@test @inferred(newindex(ci((2,2)), (true, false), (-1,-1)))  == ci((2,-1))
+@test @inferred(newindex(ci((2,2)), (false, true), (-1,-1)))  == ci((-1,2))
+@test @inferred(newindex(ci((2,2)), (false, false), (-1,-1))) == ci((-1,-1))
+@test @inferred(newindex(ci((2,2)), (true,), (-1,-1)))   == ci((2,))
+@test @inferred(newindex(ci((2,2)), (true,), (-1,)))   == ci((2,))
+@test @inferred(newindex(ci((2,2)), (false,), (-1,))) == ci((-1,))
+@test @inferred(newindex(ci((2,2)), (), ())) == ci(())
 
 end
 
diff --git a/test/offsetarray.jl b/test/offsetarray.jl
index 4a91bf6a29f26..4e83a3299226d 100644
--- a/test/offsetarray.jl
+++ b/test/offsetarray.jl
@@ -7,7 +7,7 @@
 
 module OAs
 
-using Base: DimOrInd, DimsOrInds, Indices, LinearSlow, LinearFast
+using Base: Indices, LinearSlow, LinearFast, tail
 
 export OffsetArray
 
@@ -20,7 +20,7 @@ typealias OffsetVector{T,AA<:AbstractArray} OffsetArray{T,1,AA}
 OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) = OffsetArray{T,N,typeof(A)}(A, offsets)
 OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) = OffsetArray(A, offsets)
 
-(::Type{OffsetArray{T,N}}){T,N}(inds::Indices{N}) = OffsetArray{T,N,Array{T,N}}(Array{T,N}(map(Base.dimlength, inds)), map(indsoffset, inds))
+(::Type{OffsetArray{T,N}}){T,N}(inds::Indices{N}) = OffsetArray{T,N,Array{T,N}}(Array{T,N}(map(length, inds)), map(indsoffset, inds))
 (::Type{OffsetArray{T}}){T,N}(inds::Indices{N}) = OffsetArray{T,N}(inds)
 
 Base.linearindexing{T<:OffsetArray}(::Type{T}) = Base.linearindexing(parenttype(T))
@@ -28,59 +28,62 @@ parenttype{T,N,AA}(::Type{OffsetArray{T,N,AA}}) = AA
 parenttype(A::OffsetArray) = parenttype(typeof(A))
 
 Base.parent(A::OffsetArray) = A.parent
-Base.size(A::OffsetArray) = size(parent(A))
+
+errmsg(A) = error("size not supported for arrays with indices $(indices(A)); see http://docs.julialang.org/en/latest/devdocs/offset-arrays/")
+Base.size(A::OffsetArray) = errmsg(A)
+Base.size(A::OffsetArray, d) = errmsg(A)
+Base.eachindex(::LinearSlow, A::OffsetArray) = CartesianRange(indices(A))
+Base.eachindex(::LinearFast, A::OffsetVector) = indices(A, 1)
 Base.summary(A::OffsetArray) = string(typeof(A))*" with indices "*string(indices(A))
 
 # Implementations of indices and indices1. Since bounds-checking is
 # performance-critical and relies on indices, these are usually worth
 # optimizing thoroughly.
-@inline Base.indices(A::OffsetArray, d) = 1 <= d <= length(A.offsets) ? (o = A.offsets[d]; (1+o:size(parent(A),d)+o)) : (1:1)
-@inline Base.indices(A::OffsetArray) = _indices((), A)  # would rather use ntuple, but see #15276
-@inline _indices{T,N}(out::NTuple{N}, A::OffsetArray{T,N}) = out
-@inline _indices(out, A::OffsetArray) = (d = length(out)+1; o = A.offsets[d]; _indices((out..., (1+o:size(parent(A),d)+o)), A))
-# By optimizing indices1 we can avoid a branch on the dim-check
-Base.indices1{T}(A::OffsetArray{T,0}) = 1:1
-@inline Base.indices1(A::OffsetArray) = (o = A.offsets[1]; 1+o:size(parent(A),1)+o)
+@inline Base.indices(A::OffsetArray, d) = 1 <= d <= length(A.offsets) ? indices(parent(A))[d] + A.offsets[d] : (1:1)
+@inline Base.indices(A::OffsetArray) = _indices(indices(parent(A)), A.offsets)  # would rather use ntuple, but see #15276
+@inline _indices(inds, offsets) = (inds[1]+offsets[1], _indices(tail(inds), tail(offsets))...)
+_indices(::Tuple{}, ::Tuple{}) = ()
+Base.indices1{T}(A::OffsetArray{T,0}) = 1:1  # we only need to specialize this one
 
 function Base.similar(A::OffsetArray, T::Type, dims::Dims)
     B = similar(parent(A), T, dims)
 end
 function Base.similar(A::AbstractArray, T::Type, inds::Tuple{UnitRange,Vararg{UnitRange}})
-    B = similar(A, T, map(Base.dimlength, inds))
+    B = similar(A, T, map(length, inds))
     OffsetArray(B, map(indsoffset, inds))
 end
 
-Base.similar(f::Union{Function,Type}, shape::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(f(map(Base.dimlength, shape)), map(indsoffset, shape))
+Base.similar(f::Union{Function,DataType}, shape::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(f(map(length, shape)), map(indsoffset, shape))
 
-Base.reshape(A::AbstractArray, inds::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(reshape(A, map(Base.dimlength, inds)), map(indsoffset, inds))
+Base.reshape(A::AbstractArray, inds::Tuple{UnitRange,Vararg{UnitRange}}) = OffsetArray(reshape(A, map(length, inds)), map(indsoffset, inds))
 
 @inline function Base.getindex{T,N}(A::OffsetArray{T,N}, I::Vararg{Int,N})
-    @boundscheck checkbounds(A, I...)
+    checkbounds(A, I...)
     @inbounds ret = parent(A)[offset(A.offsets, I)...]
     ret
 end
 @inline function Base._getindex(::LinearFast, A::OffsetVector, i::Int)
-    @boundscheck checkbounds(A, i)
+    checkbounds(A, i)
     @inbounds ret = parent(A)[offset(A.offsets, (i,))[1]]
     ret
 end
 @inline function Base._getindex(::LinearFast, A::OffsetArray, i::Int)
-    @boundscheck checkbounds(A, i)
+    checkbounds(A, i)
     @inbounds ret = parent(A)[i]
     ret
 end
 @inline function Base.setindex!{T,N}(A::OffsetArray{T,N}, val, I::Vararg{Int,N})
-    @boundscheck checkbounds(A, I...)
+    checkbounds(A, I...)
     @inbounds parent(A)[offset(A.offsets, I)...] = val
     val
 end
 @inline function Base._setindex!(::LinearFast, A::OffsetVector, val, i::Int)
-    @boundscheck checkbounds(A, i)
+    checkbounds(A, i)
     @inbounds parent(A)[offset(A.offsets, (i,))[1]] = val
     val
 end
 @inline function Base._setindex!(::LinearFast, A::OffsetArray, val, i::Int)
-    @boundscheck checkbounds(A, i)
+    checkbounds(A, i)
     @inbounds parent(A)[i] = val
     val
 end
@@ -101,8 +104,12 @@ let
 # Basics
 A0 = [1 3; 2 4]
 A = OffsetArray(A0, (-1,2))                   # LinearFast
-S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))  # LinearSlow
+S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))   # LinearSlow
 @test indices(A) == indices(S) == (0:1, 3:4)
+@test_throws ErrorException size(A)
+@test_throws ErrorException size(A, 1)
+
+# Scalar indexing
 @test A[0,3] == A[1] == S[0,3] == S[1] == 1
 @test A[1,3] == A[2] == S[1,3] == S[2] == 2
 @test A[0,4] == A[3] == S[0,4] == S[3] == 3
@@ -204,7 +211,6 @@ cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,10^3), (10,-9))) # neithe
 # Similar
 B = similar(A, Float32)
 @test isa(B, OffsetArray{Float32,2})
-@test size(B) == size(A)
 @test indices(B) === indices(A)
 B = similar(A, (3,4))
 @test isa(B, Array{Int,2})
@@ -307,10 +313,10 @@ cumsum!(C, A, 1)
 @test parent(cumsum(A, 1)) == cumsum(parent(A), 1)
 cumsum!(C, A, 2)
 @test parent(C) == cumsum(parent(A), 2)
-R = similar(A, (-2:-2, 6:9))
+R = similar(A, (1:1, 6:9))
 maximum!(R, A)
 @test parent(R) == maximum(parent(A), 1)
-R = similar(A, (-2:1, 6:6))
+R = similar(A, (-2:1, 1:1))
 maximum!(R, A)
 @test parent(R) == maximum(parent(A), 2)
 amin, iamin = findmin(A)