diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0a286d1..83e6a26 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -15,7 +15,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1.5' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
         os: [ubuntu-latest, windows-latest, macOS-latest] 
         arch:
diff --git a/Project.toml b/Project.toml
index 8f5990b..3864e38 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,7 @@ repo = "https://github.com/JuliaArrays/ShiftedArrays.jl.git"
 version = "2.0.0"
 
 [compat]
-julia = "1"
+julia = "1.5"
 
 [extras]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
diff --git a/docs/src/api.md b/docs/src/api.md
index d2bd672..fdb59d3 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -36,4 +36,6 @@ default
 ```@docs
 ShiftedArrays.padded_tuple
 ShiftedArrays.ft_center_diff
+ShiftedArrays.materialize_checkerboard!
+ShiftedArrays.refine_view
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 99044e4..43c5207 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -15,11 +15,11 @@ julia> v = reshape(1:16, 4, 4)
  4  8  12  16
 
 julia> s = ShiftedArray(v, (2, 0))
-4×4 ShiftedArray{Int64, Missing, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 ShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{2, 0}, Missing}:
   missing   missing    missing    missing
   missing   missing    missing    missing
  1         5          9         13
- 2         6         10         14 
+ 2         6         10         14
 ```
 
 The parent Array as well as the amount of shifting can be recovered with `parent` and `shifts` respectively.
@@ -53,7 +53,7 @@ If you pass an integer, it will shift in the first dimension:
 
 ```julia
 julia> ShiftedArray(v, 1)
-4×4 ShiftedArray{Int64, Missing, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 ShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{1, 0}, Missing}:
   missing   missing    missing    missing
  1         5          9         13
  2         6         10         14
@@ -64,7 +64,7 @@ A custom default value (other than `missing`) can be provided with the `default`
 
 ```julia
 julia> ShiftedArray([1.2, 3.1, 4.5], 1, default = NaN)
-3-element ShiftedVector{Float64, Float64, Vector{Float64}}:
+3-element ShiftedVector{Float64, Vector{Float64}, Tuple{1}, Val{NaN}}:
  NaN
    1.2
    3.1
@@ -76,7 +76,7 @@ Accessing indexes outside the `ShiftedArray` give a `BoundsError`, even if the s
 
 ```julia
 julia> ShiftedArray([1, 2, 3], 1)[4]
-ERROR: BoundsError: attempt to access 3-element ShiftedVector{Int64, Missing, Vector{Int64}} at index [4]
+ERROR: BoundsError: attempt to access 3-element ShiftedVector{Int64, Vector{Int64}, Tuple{1}, Missing} at index [4]
 ```
 
 ## Shifting the data
@@ -87,11 +87,11 @@ Using the `ShiftedArray` type, this package provides two operations for lazily s
 julia> v = [1, 3, 5, 4];
 
 julia> ShiftedArrays.lag(v)
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
+4-element ShiftedVector{Int64, Vector{Int64}, Tuple{1}, Missing}:
   missing
  1
  3
- 5       
+ 5     
 
 julia> v .- ShiftedArrays.lag(v) # compute difference from previous element without unnecessary allocations
 4-element Vector{Union{Missing, Int64}}:
@@ -101,7 +101,7 @@ julia> v .- ShiftedArrays.lag(v) # compute difference from previous element with
  -1       
 
 julia> s = ShiftedArrays.lag(v, 2) # shift by more than one element
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
+4-element ShiftedVector{Int64, Vector{Int64}, Tuple{2}, Missing}:
   missing
   missing
  1
@@ -114,7 +114,7 @@ julia> s = ShiftedArrays.lag(v, 2) # shift by more than one element
 julia> v = [1, 3, 5, 4];
 
 julia> ShiftedArrays.lead(v)
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
+4-element ShiftedVector{Int64, Vector{Int64}, Tuple{-1}, Missing}:
  3
  5
  4
@@ -134,7 +134,7 @@ Our implementation of `circshift` relies on them to avoid copying:
 julia> w = reshape(1:16, 4, 4);
 
 julia> s = ShiftedArrays.circshift(w, (1, -1))
-4×4 CircShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 CircShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{1, 3}}:
  8  12  16  4
  5   9  13  1
  6  10  14  2
diff --git a/src/ShiftedArrays.jl b/src/ShiftedArrays.jl
index 20b6f82..91c6118 100644
--- a/src/ShiftedArrays.jl
+++ b/src/ShiftedArrays.jl
@@ -3,6 +3,7 @@ module ShiftedArrays
 import Base: checkbounds, getindex, setindex!, parent, size, axes
 export ShiftedArray, ShiftedVector, shifts, default
 export CircShiftedArray, CircShiftedVector
+export CircShift
 
 include("shiftedarray.jl")
 include("circshiftedarray.jl")
diff --git a/src/circshift.jl b/src/circshift.jl
index f845ca4..4637f6b 100644
--- a/src/circshift.jl
+++ b/src/circshift.jl
@@ -12,7 +12,7 @@ remaining dimensions is assumed to be `0`.
 julia> v = [1, 3, 5, 4];
 
 julia> ShiftedArrays.circshift(v, 1)
-4-element CircShiftedVector{Int64, Vector{Int64}}:
+4-element ShiftedArray{Int64, 1, Vector{Int64}, Tuple{1}, CircShift}:
  4
  1
  3
@@ -21,7 +21,7 @@ julia> ShiftedArrays.circshift(v, 1)
 julia> w = reshape(1:16, 4, 4);
 
 julia> ShiftedArrays.circshift(w, (1, -1))
-4×4 CircShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 CircShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Tuple{1, 3}}:
  8  12  16  4
  5   9  13  1
  6  10  14  2
diff --git a/src/circshiftedarray.jl b/src/circshiftedarray.jl
index 1fa0bf2..689c69c 100644
--- a/src/circshiftedarray.jl
+++ b/src/circshiftedarray.jl
@@ -3,7 +3,7 @@
 
 Custom `AbstractArray` object to store an `AbstractArray` `parent` circularly shifted
 by `shifts` steps (where `shifts` is a `Tuple` with one `shift` value per dimension of `parent`).
-Use `copy` to collect the values of a `CircShiftedArray` into a normal `Array`.
+Use `copy` or `collect` to collect the values of a `ShiftedArray` into a normal `Array`.
 
 !!! note
     `shift` is modified with a modulo operation and does not store the passed value
@@ -19,7 +19,7 @@ Use `copy` to collect the values of a `CircShiftedArray` into a normal `Array`.
 julia> v = [1, 3, 5, 4];
 
 julia> s = CircShiftedArray(v, (1,))
-4-element CircShiftedVector{Int64, Vector{Int64}}:
+4-element ShiftedArray{Int64, 1, Vector{Int64}, CircShift}:
  4
  1
  3
@@ -33,57 +33,51 @@ julia> copy(s)
  5
 ```
 """
-struct CircShiftedArray{T, N, S<:AbstractArray} <: AbstractArray{T, N}
-    parent::S
-    # the field `shifts` stores the circular shifts modulo the size of the parent array
-    shifts::NTuple{N, Int}
-    function CircShiftedArray(p::AbstractArray{T, N}, n = ()) where {T, N}
-        shifts = map(mod, padded_tuple(p, n), size(p))
-        return new{T, N, typeof(p)}(p, shifts)
+const CircShiftedArray{T, N, A<:AbstractArray} = ShiftedArray{T, N, A, CircShift} 
+CircShiftedArray(p::AbstractArray, n=()) = ShiftedArray(p, map(mod, padded_tuple(p, n), size(p)); default=CircShift())
+function CircShiftedArray(p::ShiftedArray, n=()) 
+    ns = map(mod, padded_tuple(p, n) .+ shifts(p), size(p))
+    if all(ns.==0)
+        return p.parent
+    else
+        return ShiftedArray(p.parent, ns; default=CircShift())
     end
 end
 
-function CircShiftedArray(c::CircShiftedArray, n = ())
-    shifts = map(+, ShiftedArrays.shifts(c), padded_tuple(c, n))
-    return CircShiftedArray(parent(c), shifts)
-end
-
 """
-    CircShiftedVector{T, S<:AbstractArray}
+    CircShiftedVector{T, A<:AbstractArray}
 
-Shorthand for `CircShiftedArray{T, 1, S}`.
+Shorthand for `CircShiftedArray{T, 1, A}`.
 """
-const CircShiftedVector{T, S<:AbstractArray} = CircShiftedArray{T, 1, S}
+const CircShiftedVector{T, A<:AbstractArray} = ShiftedVector{T, A, CircShift}
 
 CircShiftedVector(v::AbstractVector, n = ()) = CircShiftedArray(v, n)
 
-size(s::CircShiftedArray) = size(parent(s))
-axes(s::CircShiftedArray) = axes(parent(s))
+has_circ_type(a::CircShiftedArray) = true
 
-@inline function bringwithin(ind_with_offset::Int, ranges::AbstractUnitRange)
-    return ifelse(ind_with_offset < first(ranges), ind_with_offset + length(ranges), ind_with_offset)
+
+# mod1 avoids first subtracting one and then adding one
+@inline function Base.getindex(csa::CircShiftedArray{T,N,A}, i::Vararg{Int,N}) where {T,N,A} 
+    getindex(csa.parent, (mod1(i[j]-shifts(csa)[j], size(csa.parent, j)) for j in 1:N)...)
 end
 
-@inline function getindex(s::CircShiftedArray{T, N}, x::Vararg{Int, N}) where {T, N}
-    @boundscheck checkbounds(s, x...)
-    v, ind = parent(s), offset(shifts(s), x)
-    i = map(bringwithin, ind, axes(s))
-    return @inbounds v[i...]
+@inline function Base.setindex!(csa::CircShiftedArray{T,N,A}, v, i::Int) where {T,N,A}
+    # @show "si circ"
+    setindex!(csa.parent, v, i)
 end
 
-@inline function setindex!(s::CircShiftedArray{T, N}, el, x::Vararg{Int, N}) where {T, N}
-    @boundscheck checkbounds(s, x...)
-    v, ind = parent(s), offset(shifts(s), x)
-    i = map(bringwithin, ind, axes(s))
-    @inbounds v[i...] = el
-    return s
+@inline function Base.setindex!(csa::CircShiftedArray{T,N,A}, v, i::Vararg{Int,N}) where {T,N,A}
+    setindex!(csa.parent, v, (mod1(i[j]-shifts(csa)[j], size(csa.parent, j)) for j in 1:N)...)
+    csa
 end
 
-parent(s::CircShiftedArray) = s.parent
+# for speed reasons use the optimized version in Base for actually perfoming the circshift in this case:
+Base.collect(csa::CircShiftedArray{T,N,A}) where {T,N,A} = Base.circshift(csa.parent, shifts(csa))
 
-"""
-    shifts(s::CircShiftedArray)
-
-Return amount by which `s` is shifted compared to `parent(s)`.
-"""
-shifts(s::CircShiftedArray) = s.shifts
+# this is not really fully in place, but the only way to emulate the reverse! function
+function Base.reverse!(csa::CircShiftedArray; dims=:)
+    tmp = Base.reverse(csa.parent; dims=dims)
+    # keep the old shift but compensate by an appropriate circshift
+    Base.circshift!(csa.parent, tmp, -2 .* shifts(csa))
+    return csa
+end
diff --git a/src/lag.jl b/src/lag.jl
index 6dde9f3..35459c5 100644
--- a/src/lag.jl
+++ b/src/lag.jl
@@ -13,7 +13,7 @@ remaining dimensions is assumed to be `0`.
 julia> v = [1, 3, 5, 4];
 
 julia> ShiftedArrays.lag(v)
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
+4-element ShiftedVector{Int64, Vector{Int64}, Missing}:
   missing
  1
  3
@@ -23,7 +23,7 @@ julia> w = 1:2:9
 1:2:9
 
 julia> s = ShiftedArrays.lag(w, 2)
-5-element ShiftedVector{Int64, Missing, StepRange{Int64, Int64}}:
+5-element ShiftedVector{Int64, StepRange{Int64, Int64}, Missing}:
   missing
   missing
  1
@@ -41,7 +41,7 @@ julia> copy(s)
 julia> v = reshape(1:16, 4, 4);
 
 julia> s = ShiftedArrays.lag(v, (0, 2))
-4×4 ShiftedArray{Int64, Missing, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 ShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Missing}:
  missing  missing  1  5
  missing  missing  2  6
  missing  missing  3  7
@@ -67,7 +67,7 @@ remaining dimensions is assumed to be `0`.
 julia> v = [1, 3, 5, 4];
 
 julia> ShiftedArrays.lead(v)
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
+4-element ShiftedVector{Int64, Vector{Int64}, Missing}:
  3
  5
  4
@@ -77,7 +77,7 @@ julia> w = 1:2:9
 1:2:9
 
 julia> s = ShiftedArrays.lead(w, 2)
-5-element ShiftedVector{Int64, Missing, StepRange{Int64, Int64}}:
+5-element ShiftedVector{Int64, StepRange{Int64, Int64}, Missing}:
  5
  7
  9
@@ -95,7 +95,7 @@ julia> copy(s)
 julia> v = reshape(1:16, 4, 4);
 
 julia> s = ShiftedArrays.lead(v, (0, 2))
-4×4 ShiftedArray{Int64, Missing, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 ShiftedArray{Int64, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Missing}:
   9  13  missing  missing
  10  14  missing  missing
  11  15  missing  missing
diff --git a/src/shiftedarray.jl b/src/shiftedarray.jl
index b36c12c..5746bde 100644
--- a/src/shiftedarray.jl
+++ b/src/shiftedarray.jl
@@ -1,25 +1,8 @@
-"""
-    padded_tuple(v::AbstractVector, s)
-
-Internal function used to compute shifts. Return a `Tuple` with as many element
-as the dimensions of `v`. The first `length(s)` entries are filled with values
-from `s`, the remaining entries are `0`. `s` should be an integer, in which case
-`length(s) == 1`, or a container of integers with keys `1:length(s)`.
+using Base
+using Core.Compiler
 
-# Examples
-
-```jldoctest padded_tuple
-julia> ShiftedArrays.padded_tuple(rand(10, 10), 3)
-(3, 0)
-
-julia> ShiftedArrays.padded_tuple(rand(10, 10), (4,))
-(4, 0)
-
-julia> ShiftedArrays.padded_tuple(rand(10, 10), (1, 5))
-(1, 5)
-```
-"""
-padded_tuple(v::AbstractArray, s) = ntuple(i -> i ≤ length(s) ? s[i] : 0, ndims(v))
+# just a type to indicate that this is a ShiftedArray rather than the ShiftedArray 
+struct CircShift end
 
 """
     ShiftedArray(parent::AbstractArray, shifts, default)
@@ -42,16 +25,7 @@ The recommended constructor is `ShiftedArray(parent, shifts; default = missing)`
 julia> v = [1, 3, 5, 4];
 
 julia> s = ShiftedArray(v, (1,))
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
-  missing
- 1
- 3
- 5
-
-julia> v = [1, 3, 5, 4];
-
-julia> s = ShiftedArray(v, (1,))
-4-element ShiftedVector{Int64, Missing, Vector{Int64}}:
+4-element ShiftedVector{Union{Missing, Int64}, Vector{Int64}, Missing}:
   missing
  1
  3
@@ -67,7 +41,7 @@ julia> copy(s)
 julia> v = reshape(1:16, 4, 4);
 
 julia> s = ShiftedArray(v, (0, 2))
-4×4 ShiftedArray{Int64, Missing, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}:
+4×4 ShiftedArray{Union{Missing, Int64}, 2, Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, Missing}:
  missing  missing  1  5
  missing  missing  2  6
  missing  missing  3  7
@@ -77,44 +51,172 @@ julia> shifts(s)
 (0, 2)
 ```
 """
-struct ShiftedArray{T, M, N, S<:AbstractArray} <: AbstractArray{Union{T, M}, N}
-    parent::S
+struct ShiftedArray{T, N, A<:AbstractArray, R} <: AbstractArray{T, N} 
+    parent::A
     shifts::NTuple{N, Int}
-    default::M
-end
 
-# low-level private constructor to handle type parameters
-function shiftedarray(v::AbstractArray{T, N}, shifts, default::M) where {T, N, M}
-    return ShiftedArray{T, M, N, typeof(v)}(v, padded_tuple(v, shifts), default)
+    function ShiftedArray(p::AbstractArray{Tb,N}, n=(); default=missing) where {Tb,N}
+        return new{shifted_array_base_type(Tb, default), N, typeof(p), to_default_type(default)}(p, padded_tuple(p, n))
+    end
+    # if a ShiftedArray is wrapped in a ShiftedArray, only a single CSA results ONLY if the default does not change! 
+    function ShiftedArray(p::ShiftedArray{Tb,N,A,R}, n=(); default=default(p)) where {Tb,N,A,R}
+        myshifts = padded_tuple(p, n)
+        if isa(default,R)
+            return new{shifted_array_base_type(Tb, default),N,A, to_default_type(default)}(p.parent, myshifts.+ shifts(p))
+        else
+            return new{eltype(p),N,typeof(p), to_default_type(default)}(p, myshifts)
+        end
+
+    end
 end
 
-function ShiftedArray(v::AbstractArray, n = (); default = ShiftedArrays.default(v))
-    return if v isa ShiftedArray && default === ShiftedArrays.default(v)
-        shifts = map(+, ShiftedArrays.shifts(v), padded_tuple(v, n))
-        shiftedarray(parent(v), shifts, default)
-    else
-        shiftedarray(v, n, default)
+shifted_array_base_type(::Type{T}, default::R) where {T,R} = Union{T,R}
+shifted_array_base_type(::Type{T}, default::CircShift) where {T} = T
+
+to_default_type(default::Val)=typeof(default)
+to_default_type(default::Type)=default
+to_default_type(default::Missing) = Missing
+to_default_type(default::Nothing) = Nothing
+to_default_type(default::CircShift) = CircShift
+to_default_type(default)=typeof(Val(default))
+
+
+"""
+    default(s::ShiftedArray)
+
+Return default value.
+"""
+default(s::ShiftedArray{T,N,A,R}) where {T,N,A,R} = default(R)
+default(::AbstractArray) = missing
+
+function default(R1,R2)
+    if R1 != R2
+        error("propagating multiple arrays in one expression which contains only shifted arrays all need the same default.")
     end
+    R1
 end
+# default(Datatype) should NOT be called.
+# default(::Type{T}) where T = T
+default(::Missing) = missing
+default(::Type{Missing}) = missing
+default(::Nothing) = nothing
+default(::Type{Nothing}) = nothing
+default(::Type{Val{N}}) where N = N
+# default(::T) where T = T
+default(::Type{ShiftedArray{T,N,A,R}}) where {T,N,A,R} = default(R)
+default(R::Type{CircShift}) = CircShift()
+default(R1, R2::Type{CircShift}) = R1()
+default(R1::Type{CircShift}, R2) = R2()
+default(R1::Type{CircShift}, R2::Type{CircShift}) = CircShift()
+# default applied to a Broadcast calculates the default value
+function default(bc::Base.Broadcast.Broadcasted)
+    bcd = replace_by_default_broadcast(bc)
+    #@show bcd
+    dv = collect(bcd)
+    #@show dv
+    if prod(size(dv)) != 1
+        error("wrong size of dv")
+    end
+    return dv[1]
+end
+
+# find out, whether any of the arguments is a circ-shifted array
+has_circ_type(a) = false
+has_circ_type(bc::Base.Broadcast.Broadcasted) = any(has_circ_type.(bc.args))
+
+# low-level private constructor to handle type parameters
+function shiftedarray(v::AbstractArray{T, N}, shifts, r=default(v)) where {T, N}
+    return ShiftedArray(v, padded_tuple(v, shifts); default=r)
+end
+
+"""
+    padded_tuple(v::AbstractVector, s)
+
+Internal function used to compute shifts. Return a `Tuple` with as many element
+as the dimensions of `v`. The first `length(s)` entries are filled with values
+from `s`, the remaining entries are `0`. `s` should be an integer, in which case
+`length(s) == 1`, or a container of integers with keys `1:length(s)`.
+
+# Examples
+
+```jldoctest padded_tuple
+julia> ShiftedArrays.padded_tuple(rand(10, 10), 3)
+(3, 0)
+
+julia> ShiftedArrays.padded_tuple(rand(10, 10), (4,))
+(4, 0)
+
+julia> ShiftedArrays.padded_tuple(rand(10, 10), (1, 5))
+(1, 5)
+```
+"""
+padded_tuple(v::AbstractArray, s) = ntuple(i -> i ≤ length(s) ? s[i] : 0, ndims(v))
+
+"""
+    shifts(s::ShiftedArray)
+
+Return amount by which `s` is shifted compared to `parent(s)`.
+"""
+@inline shifts(arr::ShiftedArray{T,N,A,R}) where {T,N,A,R} = arr.shifts
+
+# to_tuple(S::Type{T}) where {T<:Tuple}= tuple(S.parameters...)
+# shifts(arr::Type{ShiftedArray}) = arr.shifts
+
 
 """
     ShiftedVector{T, S<:AbstractArray}
 
-Shorthand for `ShiftedArray{T, 1, S}`.
+Shorthand for `ShiftedArray{T, 1, A, M}`.
 """
-const ShiftedVector{T, M, S<:AbstractArray} = ShiftedArray{T, M, 1, S}
+const ShiftedVector{T, N, A<:AbstractArray, M} = ShiftedArray{T, 1, A, M}
 
 function ShiftedVector(v::AbstractVector, n = (); default = ShiftedArrays.default(v))
     return ShiftedArray(v, n; default = default)
 end
 
-size(s::ShiftedArray) = size(parent(s))
-axes(s::ShiftedArray) = axes(parent(s))
-
 # Computing a shifted index (subtracting the offset)
 offset(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) where {N} = map(-, inds, offsets)
 
-@inline function getindex(s::ShiftedArray{<:Any, <:Any, N}, x::Vararg{Int, N}) where {N}
+# we keep this for compatability reasons
+@inline function bringwithin(ind_with_offset::Int, ranges::AbstractUnitRange)
+    return ifelse(ind_with_offset < first(ranges), ind_with_offset + length(ranges), ind_with_offset)
+end
+
+# wraps shifts into the range 0...N-1
+wrapshift(shift::NTuple, dims::NTuple) = ntuple(i -> mod(shift[i], dims[i]), length(dims))
+# wraps indices into the range 1...N
+wrapids(shift::NTuple, dims::NTuple) = ntuple(i -> mod1(shift[i], dims[i]), length(dims))
+invert_rng(s, sz) = wrapshift(sz .- s, sz)
+
+# define a new broadcast style. Stores dimensions (N), Shift(S) and default type (R)
+struct ShiftedArrayStyle{N} <: Base.Broadcast.AbstractArrayStyle{N} end
+
+# convenient constructor
+ShiftedArrayStyle{N}(::Val{M}) where {N,M} = ShiftedArrayStyle{max(N,M)}()
+
+# make it known to the system
+function Base.Broadcast.BroadcastStyle(::Type{T}) where (T<: ShiftedArray)
+    #@show "Style SA"
+    #@show shifts(T)
+     ShiftedArrayStyle{ndims(T)}()
+end
+
+# make subarrays (views) of ShiftedArray also broadcast inthe ShiftedArray style:
+Base.Broadcast.BroadcastStyle(::Type{SubArray{T,N,P,I,L}}) where {T,N,P<:ShiftedArray,I,L} = ShiftedArrayStyle{ndims(P)}()
+# ShiftedArray and default Array broadcast to ShiftedArray with max dimensions between the two
+Base.Broadcast.BroadcastStyle(::ShiftedArrayStyle{N}, ::Base.Broadcast.DefaultArrayStyle{M}) where {N,M} = ShiftedArrayStyle{max(N,M)}() 
+Base.Broadcast.BroadcastStyle(::Base.Broadcast.DefaultArrayStyle{M}, ::ShiftedArrayStyle{N}) where {M,N} = ShiftedArrayStyle{max(N,M)}() 
+Base.Broadcast.BroadcastStyle(::ShiftedArrayStyle{N}, ::Base.Broadcast.AbstractArrayStyle{M}) where {N,M} = ShiftedArrayStyle{max(N,M)}() 
+Base.Broadcast.BroadcastStyle(::Base.Broadcast.AbstractArrayStyle{M}, ::ShiftedArrayStyle{N}) where {M,N} = ShiftedArrayStyle{max(N,M)}() 
+
+@inline Base.size(csa::ShiftedArray) = size(csa.parent)
+@inline Base.size(csa::ShiftedArray, d::Int) = size(csa.parent, d)
+@inline Base.axes(csa::ShiftedArray) = axes(csa.parent)
+@inline Base.IndexStyle(::Type{<:ShiftedArray}) = IndexLinear()
+@inline Base.parent(csa::ShiftedArray) = csa.parent
+
+@inline function Base.getindex(s::ShiftedArray{<:Any, N, <:Any, <:Any}, x::Vararg{Int, N}) where {N}
+    #@show "gi shifted"
     @boundscheck checkbounds(s, x...)
     v, i = parent(s), offset(shifts(s), x)
     return if checkbounds(Bool, v, i...)
@@ -124,20 +226,461 @@ offset(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) where {N} = map(-, inds, off
     end
 end
 
-parent(s::ShiftedArray) = s.parent
+# linear indexing ignores the shifts
+@inline function Base.getindex(csa::ShiftedArray{T,N,A,R}, i::Int) where {T,N,A,R} 
+    getindex(csa.parent, i)
+end
+
+@inline function Base.setindex!(csa::ShiftedArray{T,N,A,R}, v, i::Int) where {T,N,A,R}
+    # note that we simply use the cyclic method, since the missing values are simply ignored
+    setindex!(csa.parent, v, i)
+end
+
+@inline function Base.setindex!(csa::ShiftedArray{T,N,A,R}, v, i::Vararg{Int,N}) where {T,N,A,R}
+    # note that we simply use the cyclic method, since the missing values are simply ignored
+    setindex!(csa.parent, v, (mod1(i[j]-shifts(csa)[j], size(csa.parent, j)) for j in 1:N)...)
+    csa
+end
+
+# setting a value which corresponds to the border type is ignored
+@inline function Base.setindex!(csa::ShiftedArray{T,N,A,R}, v::R, i::Vararg{Int,N}) where {T,N,A,R}
+    csa
+end
+
+# These apply for broadcasted assignment operations, if the shifts are identical
+@inline function Base.Broadcast.materialize!(dest::ShiftedArray{T,N,A,R}, src::ShiftedArray) where {T,N,A,R} 
+    if shifts(dest) != shifts(src)
+        error("Copyiing into a ShiftedArray of different shift is disallowed. Use the same shift or an ordinary array.")
+    end
+    Base.Broadcast.materialize!(dest.parent, src.parent)
+end
+
+function Base.Broadcast.copyto!(dest::AbstractArray, bc::Base.Broadcast.Broadcasted{ShiftedArrayStyle{N}}) where {N}
+     Base.Broadcast.materialize!(dest, bc)
+     return dest
+end
+
+# This copy operation is performed when an expression (in bc) needs to be evaluated
+function Base.Broadcast.copy(bc::Base.Broadcast.Broadcasted{ShiftedArrayStyle{N}}) where {N}
+    if only_shifted(bc)
+        ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
+        # since there are only shifted arrays we can use the default
+        if has_circ_type(bc)
+            res = similar(bc, ElType, size(bc); shifts = shifts(bc), default = CircShift)
+        else
+            # calculate the default here
+            res = similar(bc, ElType, size(bc); shifts = shifts(bc), default = default(bc))
+        end
+    else
+        ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
+        bcr = remove_sa_broadcast(bc)
+        res = similar(bcr, ElType, size(bcr))
+    end
+    return copyto!(res, bc)
+end
+
+# remove all the (circ-)shift part if all shifts are the same (or constants)
+# @inline function materialize!(dest::ShiftedArray{T, N, A, R}, bc::Base.Broadcast.Broadcasted{ShiftedArrays.ShiftedArrayStyle{N, S}}) where {T, N, A, S, N, S, R}
+@inline function Base.Broadcast.materialize!(dest::ShiftedArray{T, N, A, R}, bc::Base.Broadcast.Broadcasted{ShiftedArrayStyle{N}}) where {T, N, A, R}
+    #     # get all not-shifted arrays and apply the materialize operations piecewise using array views
+    return materialize_checkerboard!(dest.parent, bc, Tuple(1:N), shifts(dest), true; array_type=R)
+end
+
+# we cannot specialize the Broadcast style here, since the rhs may not contain a ShiftedArray and still wants to be assigned
+@inline function Base.Broadcast.materialize!(dest::ShiftedArray{T,N,A,R}, bc::Base.Broadcast.Broadcasted{BT}) where {T,N,A,R,BT}
+    materialize_checkerboard!(dest.parent, bc, Tuple(1:N), shifts(dest), true; array_type=R)
+    return dest
+end
+
+# for ambiguous conflict resolution
+@inline function Base.Broadcast.materialize!(dest::ShiftedArray{T,N,A,R}, bc::Base.Broadcast.Broadcasted{ShiftedArrayStyle{N2}}) where {T,N,A,R,N2}
+    if shifts(dest) != S2
+        error("assigment of ShiftedArray to another ShiftedArray needs to have the same shift.")
+    end
+    materialize_checkerboard!(dest.parent, bc, Tuple(1:N), wrapshift(size(dest) .- shifts(dest), size(dest)), true; default=R)
+    return dest
+end
+
+@inline function Base.Broadcast.materialize!(dest::AbstractArray, bc::Base.Broadcast.Broadcasted{ShiftedArrayStyle{N}}) where {N}
+    materialize_checkerboard!(dest, bc, Tuple(1:N), nothing, false; array_type=Nothing)
+    return dest
+end
+
+# This collect function needs to be defined to prevent the linear indexing to take over and just copy the raw (not shifted) data.
+function Base.collect(src::ShiftedArray{T,N,A,R}) where {T,N,A,R}
+    src = if (isa(src.parent,ShiftedArray))
+        ShiftedArray(collect(src.parent), shifts(src); default=default(src))
+    else
+        src
+    end
+    bc = Base.Broadcast.broadcasted(identity, src)
+    dest = similar(src.parent, eltype(src))
+    materialize_checkerboard!(dest, bc, Tuple(1:N), shifts(src), false; array_type=R)
+    return dest
+end
+
+# needs to generate both ranges as both appear in mixed broadcasting expressions
+function generate_shift_ranges(dest, myshift)
+    shift_rng_1 = ntuple((d)->firstindex(dest,d):firstindex(dest,d)+myshift[d]-1, ndims(dest))
+    shift_rng_2 = ntuple((d)->firstindex(dest,d)+myshift[d]:lastindex(dest,d), ndims(dest))
+    noshift_rng_1 = ntuple((d)->lastindex(dest,d)-myshift[d]+1:lastindex(dest,d), ndims(dest))
+    noshift_rng_2 = ntuple((d)->firstindex(dest,d):lastindex(dest,d)-myshift[d], ndims(dest))
+    return ((shift_rng_1, shift_rng_2), (noshift_rng_1, noshift_rng_2))
+end
 
 """
-    shifts(s::ShiftedArray)
+    materialize_checkerboard!(dest, bc, dims, myshift, dest_is_cs_array=false; array_type=CircShift) 
+
+this function subdivides the array into tiles, which each needs to be processed individually via calls to `materialize!`.
+
+|--------|
+| a| b   |
+|--|-----|---|
+| c| dD  | C |
+|--+-----|---|
+   | B   | A |
+   |---------|
 
-Return amount by which `s` is shifted compared to `parent(s)`.
 """
-shifts(s::ShiftedArray) = s.shifts
+function materialize_checkerboard!(dest, bc, dims, myshift, dest_is_cs_array=false; array_type=CircShift) 
+    dest = refine_view(dest)
+
+    if isnothing(myshift)
+        myshift = shifts(bc)
+    end
+    # gets Tuples of Tuples of 1D ranges (low and high) for each dimension
+    cs_rngs, ns_rngs = generate_shift_ranges(dest, wrapshift(size(dest) .- myshift, size(dest)))
+    nonflipped = nonflipped_source(myshift)
+    # use the broadcast, with the ShiftedArrays which are not CircShiftedArrays replaced by the default values 
+    bcd = replace_by_default_broadcast(bc)
+    for n in CartesianIndices(ntuple((x)->2, ndims(dest)))
+        cs_rng = Tuple(cs_rngs[n[d]][d] for d=1:ndims(dest))
+        ns_rng = Tuple(ns_rngs[n[d]][d] for d=1:ndims(dest))
+        dst_rng = ifelse(dest_is_cs_array, cs_rng, ns_rng)
+        dst_rng = refine_shift_rng(dest, dst_rng)
+        dst_view = @view dest[dst_rng...]
+
+        bc1 = split_array_broadcast(bc, ns_rng, cs_rng)
+        if (prod(size(dst_view)) > 0)
+            if (Tuple(n) == nonflipped) # array_type <: CircShift || 
+                @inbounds Base.Broadcast.materialize!(dst_view, bc1)
+            else
+                # check if there is any need to calculate the other quadrants.
+                if (! isa(dest, ShiftedArray) || default(dest) == CircShift())
+                    bc2 = split_array_broadcast(bcd, ns_rng, cs_rng)
+                    @inbounds Base.Broadcast.materialize!(dst_view, bc2)
+                end
+            end
+        end
+    end
+end
+function materialize_checkerboard!(dest::ShiftedArray{T,N,A,R}, bc, dims, myshift, dest_is_cs_array=(R<:CircShift); array_type=R) where {T,N,A,R}
+    @show "materialize_checkerboard SA"
+    materialize_checkerboard!(dest, bc, dims, myshift, dest_is_cs_array; array_type=array_type) 
+end
+
+# identifies which quadrant corresponds to the original (non-flipped quadrant)
+function nonflipped_source(myshift)
+    (myshift.<=0) .+ 1
+end
+
+# some code which determines whether all arrays are shifted
+@inline only_shifted(bc)  = true
+@inline only_shifted(bc::Number)  = true
+@inline only_shifted(bc::AbstractArray)  = false
+@inline only_shifted(bc::ShiftedArray)  = true
+@inline only_shifted(bc::Base.Broadcast.Broadcasted) = all(only_shifted.(bc.args))
+@inline only_shifted(bc::Base.Broadcast.Extruded) = only_shifted(bc.x)
+
+# These functions remove the ShiftArray in a broadcast and replace each by a view into the original array 
+@inline split_array_broadcast(bc, noshift_rng, shift_rng) = bc
+@inline split_array_broadcast(bc::Number, noshift_rng, shift_rng) = bc
+@inline split_array_broadcast(bc::AbstractArray, noshift_rng, shift_rng) = @view bc[noshift_rng...]
+@inline split_array_broadcast(bc::ShiftedArray, noshift_rng, shift_rng) = @view bc.parent[shift_rng...]
+@inline split_array_broadcast(bc::ShiftedArray{T,N,A,R}, noshift_rng, shift_rng)  where {T,N,A,R} =  @view bc.parent[noshift_rng...]
+@inline function split_array_broadcast(v::SubArray{T,N,P,I,L}, noshift_rng, shift_rng) where {T,N,P<:ShiftedArray,I,L}    
+    new_cs = refine_view(v)
+    new_shift_rng = refine_shift_rng(v, shift_rng)
+    res = split_array_broadcast(new_cs, noshift_rng, new_shift_rng)
+    return res
+end
+
+function split_array_broadcast(bc::Base.Broadcast.Broadcasted, noshift_rng, shift_rng)
+    # Ref below protects the argument from broadcasting
+    bc_modified = split_array_broadcast.(bc.args, Ref(noshift_rng), Ref(shift_rng))
+    res = Base.Broadcast.broadcasted(bc.f, bc_modified...)
+    return res
+end
+
+# These function remove the ShiftedArray properties from the Broadcast chain
+@inline remove_sa_broadcast(bc::Number) = bc
+@inline remove_sa_broadcast(bc::AbstractArray) = bc
+@inline remove_sa_broadcast(bc::ShiftedArray) = bc.parent
+function remove_sa_broadcast(bc::Base.Broadcast.Broadcasted)
+    # Ref below protects the argument from broadcasting
+    bc_modified = remove_sa_broadcast.(bc.args)
+    res = Base.Broadcast.broadcasted(bc.f, bc_modified...)
+    return res
+end
+
+# leave all circshifted arrays in place but replace the shifted array by a single default number
+@inline replace_by_default_broadcast(arg) = arg
+@inline replace_by_default_broadcast(sa::ShiftedArray{T,N,A,R}) where {T,N,A,R} = ifelse(R<:CircShift, sa, default(sa))
+function replace_by_default_broadcast(bc::Base.Broadcast.Broadcasted)
+    # Ref below protects the argument from broadcasting
+    bc_modified = replace_by_default_broadcast.(bc.args)
+    res = Base.Broadcast.broadcasted(bc.f, bc_modified...)
+    return res
+end
+
+# leave all circshifted arrays in place but replace the shifted array by a single default number
+shifts(anything) = missing
+shifts(::AbstractArray{T,N}) where {T,N} = missing
+function shifts(bc::Base.Broadcast.Broadcasted)
+    # Ref below protects the argument from broadcasting
+    all_shifts = shifts.(bc.args)
+    first_shift = coalesce(all_shifts...)
+    if !all(skipmissing(all_shifts) .== Ref(first_shift))
+        error("Shifts of arrays during broadcast differ. Please use `collect`.")
+    end
+    return first_shift
+end
+
+@inline function refine_shift_rng(v::SubArray{T,N,P,I,L}, shift_rng) where {T,N,P,I,L}    
+    new_shift_rng = ntuple((d)-> ifelse(isa(v.indices[d], Base.Slice), shift_rng[d], Base.Colon()), ndims(v.parent))
+    return new_shift_rng
+end
+@inline refine_shift_rng(v, shift_rng) = shift_rng
 
 """
-    default(s::ShiftedArray)
+    function refine_view(v::SubArray{T,N,P,I,L}, shift_rng)
 
-Return default value.
+returns a refined view of a SubArray of a ShiftedArray as a ShiftedArray, if necessary. Otherwise just the original array.
+find out, if the range of this view crosses any boundary of the parent ShiftedArray
+by calculating the new indices
+if, so though an error. find the full slices, which can stay a (circ-)shifted array withs shifts
 """
-default(s::ShiftedArray) = s.default
+function refine_view(v::SubArray{T,N,P,I,L}) where {T,N,P<:ShiftedArray,I,L}
+    myshift = shifts(v.parent)
+    sz = size(v.parent)
+    # find out, if the range of this view crosses any boundary of the parent ShiftedArray
+    # by calculating the new indices
+    # if, so though an error.
+    # find the full slices, which can stay a circ shifted array withs shifts
+    sub_rngs = ntuple((d)-> !isa(v.indices[d], Base.Slice), ndims(v.parent))
+
+    # in the line below one should better use "begin" instead of "1" but this is not supported by early Julia versions.
+    new_ids_begin = wrapids(ntuple((d)-> v.indices[d][1] .- myshift[d], ndims(v.parent)), sz)
+    new_ids_end = wrapids(ntuple((d)-> v.indices[d][end] .- myshift[d], ndims(v.parent)), sz)
+    if any(sub_rngs .& (new_ids_end .< new_ids_begin))
+        error("a view of a shifted array is not allowed to cross boarders of the original array. Do not use a view here.")
+        # potentially this can be remedied, once there is a decent CatViews implementation
+    end
+    new_rngs = ntuple((d)-> ifelse(isa(v.indices[d],Base.Slice), v.indices[d], new_ids_begin[d]:new_ids_end[d]), ndims(v.parent))
+    new_shift = ntuple((d)-> ifelse(isa(v.indices[d],Base.Slice), 0, myshift[d]), ndims(v.parent))
+    new_cs = ShiftedArray((@view v.parent.parent[new_rngs...]), new_shift)
+    return new_cs
+end
+
+refine_view(csa::AbstractArray) = csa
+
+function Base.reverse(csa::ShiftedArray; dims=:)
+    rev = Base.reverse(csa.parent; dims=dims)
+    ShiftedArray(rev, .-shifts(csa), default=default(csa))
+end
+
+# these array isequal and == functions are defined to be compatible with the previous definition of equality (equal values only)
+function Base.isequal(csa::ShiftedArray, arr::AbstractArray)
+    if isequal(Ref(csa), Ref(arr))
+        return true
+    end
+    return all(isequal.(csa, arr))
+end
+Base.isequal(arr::AbstractArray, csa::ShiftedArray) = isequal(csa, arr)
+function Base.isequal(arr::ShiftedArray, csa::ShiftedArray)
+    if isequal(Ref(csa), Ref(arr))
+        return true
+    end
+    if isequal(default(arr), CircShift()) && isequal(default(csa), CircShift()) && shifts(arr)==shifts(csa)
+        return isequal(arr.parent, csa.parent)
+    end
+    return all(isequal.(csa, arr))
+end
+
+# cases for views of ShiftedArray should still be added at some point.
+
+function Base. ==(csa::ShiftedArray, arr::AbstractArray) 
+    if isequal(Ref(csa), Ref(arr))
+        return true
+    end
+    return all(.==(csa, arr))
+end
+Base. ==(arr::AbstractArray, csa::ShiftedArray) = (csa==arr)
+
+function Base. ==(arr::ShiftedArray, csa::ShiftedArray)
+    if isequal(Ref(csa), Ref(arr))
+        return true
+    end
+    if default(arr)==CircShift() && default(csa) == CircShift() && shifts(arr)==shifts(csa)
+        return arr.parent == csa.parent
+    end
+    return all(.==(csa, arr))
+end
+
+function Base.isapprox(csa::ShiftedArray, arr::AbstractArray; kwargs...) 
+    if isequal(Ref(csa), Ref(arr))
+        return true
+    end
+    return isapprox(collect(csa), arr; kwargs...)
+end
+Base.isapprox(arr::AbstractArray, csa::ShiftedArray; kwargs...) = isapprox(csa, arr; kwargs...)
+
+function Base.isapprox(arr::ShiftedArray, csa::ShiftedArray; kwargs...)
+    if isequal(Ref(csa), Ref(arr))
+        return true
+    end
+    if default(arr)==CircShift() && default(csa) == CircShift() && shifts(arr)==shifts(csa)
+        return isapprox(arr.parent, csa.parent; kwargs...)
+    end
+    return isapprox(collect(csa), arr; kwargs...)
+end
+
+# Base.isequal(arr::AbstractArray, csa::ShiftedArray) = isequal(csa, arr)
+# Base.isequal(csa1::ShiftedArray, csa2::ShiftedArray)  = invoke(isequal, Tuple{ShiftedArray, AbstractArray}, csa1, csa2)
+# Base. ==(csa::ShiftedArray, arr::AbstractArray) = isequal(csa, arr)
+# Base. ==(arr::AbstractArray, csa::ShiftedArray) = isequal(csa, arr)
+# Base. ==(csa1::ShiftedArray, csa2::ShiftedArray) = isequal(csa1,csa2)
+
+function Base.copy(arr::ShiftedArray)
+    collect(arr)
+end
+
+Base.eltype(arr::ShiftedArray{T,N,A,R})  where {T,N,A,R<:CircShift} = eltype(parent(arr))
+
+# The code below supports the default-types for ShiftedArrays interacting with constants, which remain a ShiftedArray. 
+propagate_default(f,a,b) = f(a,b)
+propagate_default(f, a::Missing, b::AbstractArray) = missing
+propagate_default(f, a::AbstractArray, b::Missing) = missing
+propagate_default(f, a::Nothing, b::AbstractArray) = nothing
+propagate_default(f, a::AbstractArray, b::Nothing) = nothing
+#propagate_default(f, a::Base.Broadcast.Broadcasted, b) = Base.promote_op(Base.broadcast, typeof(f), typeof(a), typeof(b))
+#propagate_default(f, a, ::Base.Broadcast.Broadcasted) = a
+propagate_default(f, a::CircShift,b) = CircShift()
+propagate_default(f, a, b::CircShift) = CircShift()
+propagate_default(f, a::CircShift, b::CircShift) = CircShift()
+
+function Base.broadcasted(f::Function, a::ShiftedArray, b::Number)
+    if isa(b, Base.Broadcast.Broadcasted)
+        a = collect(a)
+        b = reshape(collect(b), size(b))
+        return Base.broadcasted(f, a, b)
+    else
+        new_default = propagate_default(f,default(a), b)
+        a = ShiftedArray(a.parent, shifts(a), default=new_default)
+    end
+    return invoke(Base.broadcasted, Tuple{typeof(f), AbstractArray, typeof(b)}, f, a, b)
+end
+
+function Base.broadcasted(f::Function, b::Number, a::ShiftedArray)
+    if isa(b, Base.Broadcast.Broadcasted)
+        a = collect(a)
+        b = collect(b)
+    else
+        new_default = propagate_default(f, b, default(a))
+        a = ShiftedArray(a.parent, shifts(a), default=new_default)
+    end
+    return invoke(Base.broadcasted, Tuple{typeof(f), typeof(b), AbstractArray}, f, b, a)
+end
+
+function Base.broadcasted(f::Function, a::ShiftedArray, b::ShiftedArray)
+    new_default = propagate_default(f,default(a), default(b))
+    if shifts(a) != shifts(b)
+        error("shifts ($(shifts(a)) and $(shifts(b))) of both arrays need to be equal.")
+    end
+    raw = f.(a.parent, b.parent)
+    res = ShiftedArray(raw, shifts(a), default=new_default)
+    return res
+end
+
+function Base.similar(arr::ShiftedArray, eltp::Type, dims::Tuple{Int64, Vararg{Int64, N}}) where {N}
+    na = similar(arr.parent, eltp, dims)
+    # the results-type depends on whether the result size is the same or not. Same size can remain ShiftedArray.
+    # its important that a shifted array is the result for reductions on CircShiftedArray, since only then the broadcasting works
+    # Since similar cannot infer the shift information when a sub-range ws created e.g. sv[1:3] 
+    return na # ifelse(size(arr)==dims, na, ShiftedArray(na, shifts(arr); default=default(arr)))
+end
+
+# specialized version for CircShiftArray which removes the Union Type
+function Base.similar(arr::ShiftedArray, ::Type{Union{CircShift,T}}, dims::Tuple{Int64, Vararg{Int64, N}}) where {T,N}
+    na = similar(arr.parent, T, dims)
+    return na 
+end
+
+Base.similar(arr::ShiftedArray, eltp::Type) = similar(arr, eltp, size(arr))
+Base.similar(arr::ShiftedArray, dims::Tuple{Int64, Vararg{Int64, N}}) where {N} = similar(arr, eltype(arr), dims)
+Base.similar(arr::ShiftedArray) = similar(arr, eltype(arr), size(arr))
+
+# this should be used for sum, max etc.
+# function Base.reduce(op, itr::ShiftedArray; init=init)
+#     @show "reduce"
+#     reduce(op, itr.parent; init=init)
+# end
+
+# The order of non-broadcast-reduction does not matter, but ONLY, if the dims keyword is not used.
+function Base.mapreduce(mapf, op, a::ShiftedArray; dims=:, kw...)
+    res = mapreduce(mapf, op, a.parent; dims=dims, kw...)
+    if !isa(dims, Colon)
+        return ShiftedArray(res, shifts(a), default=default(a))
+    end
+    return res
+end
+
+Base.any(a::ShiftedArray; dims=:) = any(a.parent; dims)
+Base.all(a::ShiftedArray; dims=:) = all(a.parent; dims)
+
+# function Base.mapreducedim!(f, op, B::ShiftedArray, A)
+#     @show "reducedim! 1"
+#     Base.mapreducedim!(f, op, B.parent, A)
+# end
+
+# function Base.mapreducedim!(f, op, B::ShiftedArray, A::ShiftedArray)
+#     @show "reducedim! 2"
+#     Base.mapreducedim!(f, op, B.parent, A.parent)
+# end
+
+# This one is called for reduce operations to allocate like similar, but always a ShiftedArray
+function Base.reducedim_initarray(arr::ShiftedArray, region, init, r::Type{R}) where R
+    @show "reducedim_initarray"
+    res = invoke(Base.reducedim_initarray, Tuple{AbstractArray, typeof(region), typeof(init), typeof(r)}, arr, region,init,r)
+    return ShiftedArray(res, shifts(arr), default=default(arr))
+end
+
+function Base.similar(bc::Base.Broadcast.Broadcasted{ShiftedArrayStyle{N},Ax,F,Args}, et::ET, dims::Any; shifts=shifts(bc), default=default(bc)) where {N,ET,Ax,F,Args}
+    # remove the ShiftedArrayStyle from broadcast to call the original "similar" function 
+    bc_tmp = remove_sa_broadcast(bc) #Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{N}}(bc.f, bc.args, bc.axes)
+    res = invoke(Base.Broadcast.similar, Tuple{typeof(bc_tmp),ET,Any}, bc_tmp, et, dims)
+    # if all broadcast members are shifted, also return a shifted version
+    if only_shifted(bc)
+        return ShiftedArray(res, shifts; default=default)
+    else
+        return res
+    end
+end
+
+function Base.show(io::IO, mm::MIME"text/plain", cs::ShiftedArray) 
+    # a bit of a hack to determine whether the datatype is CuArray without needing a CUDA.jl dependence
+    if startswith(string(typeof(cs.parent)),"CuArray")
+        # this is needed such that the show method does not throw an error when individual element access is used
+        buffer = IOBuffer()
+        ioc = IOContext(buffer, io)
+        # unfortunately we have to collect here
+        show(ioc, mm, collect(cs))
+        buffer = String(take!(buffer))
+        lines = split(buffer,"\n")
+        lines[1] = string(typeof(cs)) * ":"
+        print(io, join(lines,"\n"))
+    else
+        invoke(Base.show, Tuple{IO, typeof(mm), AbstractArray}, io, mm, cs) 
+    end
+end
 
-default(::AbstractArray) = missing
\ No newline at end of file
diff --git a/test/benchmark.jl b/test/benchmark.jl
new file mode 100644
index 0000000..afc025b
--- /dev/null
+++ b/test/benchmark.jl
@@ -0,0 +1,40 @@
+using ShiftedArrays
+using BenchmarkTools
+sz = (1024,1024)
+
+v = rand(sz...);
+useCuda = true
+if useCuda
+    using CUDA
+    CUDA.allowscalar(false);
+    v = CuArray(v)
+
+    macro mytime(expr)
+        return :( @btime CUDA.@sync $expr)
+    end    
+else
+    macro mytime(expr)
+        return :( @btime $expr)
+    end
+end
+sh = (335, 444)
+
+sv = ShiftedArray(v, sh)
+cv = CircShiftedArray(v, sh)
+
+
+# timings stated for Dell Laptop XPS 15 (i7 11800) on Windows 10
+@mytime q = $sv .+ 5.0; # bc version: 1.48 ms, CuArray bc: 0.099 ms, old version: 2.73 ms
+res = sv .+ 5.0
+
+@mytime res .= $sv .+ 5.0 # bc version: 0.18 ms, CuArray bc: 0.097 ms, old version: 0.37 ms
+
+sv = ShiftedArray(v, sh, default=0.0)
+resn = copy(v)
+@mytime $resn .= $sv .+ 5.0; # bc version: 0.28 ms, CuArray bc: 0.118 ms, old version: 0.42 ms
+
+@mytime $res .= $sv .+ 5.0 .* $v .*$sv; # bc version: 0.727 ms, CuArray bc: 0.264 ms, old version: 1.65 ms
+
+svi = ShiftedArrays.ifftshift(v)
+@mytime $resn .= $svi .+ 5.0 .* $v .*$svi; # bc version: 0.53 ms, CuArray bc: 0.324 ms, old version: 3.98 ms
+@mytime $resn .= ShiftedArrays.fftshift($svi .+ 5.0 .* $v .*$svi); # bc version: 2.41 ms,  CuArray bc: 0.443 ms, old version: 3.98 ms
diff --git a/test/runtests.jl b/test/runtests.jl
index 4415ed5..125d22f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,31 +5,33 @@ using AbstractFFTs
     v = [1, 3, 5, 4]
     @test all(v .== ShiftedVector(v))
     sv = ShiftedVector(v, -1)
-    @test isequal(sv, ShiftedVector(v, (-1,)))
+    @test isequal(sv, ShiftedVector(v, (-1,)))   
     @test length(sv) == 4
-    @test all(sv[1:3] .== [3, 5, 4])
+    # causes individual getindex:
+    @test all(sv[1:3] .== [3, 5, 4]) 
     @test ismissing(sv[4])
-    diff = v .- sv
+    diff = v .- sv;
     @test isequal(diff, [-2, -2, 1, missing])
     @test shifts(sv) == (-1,)
-    svneg = ShiftedVector(v, -1, default = -100)
+    svneg = ShiftedVector(v, -1, default = -100);
     @test default(svneg) == -100
     @test copy(svneg) == coalesce.(sv, -100)
-    @test isequal(sv[1:3], Union{Int64, Missing}[3, 5, 4])
-    svnest = ShiftedVector(ShiftedVector(v, 1), 2)
-    sv = ShiftedVector(v, 3)
+    # causes individual getindex:
+    @test isequal(sv[1:3], Union{Int64, Missing}[3, 5, 4]) 
+    svnest = ShiftedVector(ShiftedVector(v, 1), 2);
+    sv = ShiftedVector(v, 3);
     @test sv === svnest
-    sv = ShiftedVector(v, 2, default = nothing)
-    sv1 = ShiftedVector(sv, 1)
-    sv2 = ShiftedVector(sv, 1, default = 0)
+    sv = ShiftedVector(v, 2, default = nothing);
+    sv1 = ShiftedVector(sv, 1);
+    sv2 = ShiftedVector(sv, 1, default = 0);
     @test isequal(collect(sv1), [nothing, nothing, nothing, 1])
     @test isequal(collect(sv2), [0, nothing, nothing, 1])
 end
 
 @testset "ShiftedArray" begin
-    v = reshape(1:16, 4, 4)
+    v = reshape(1:16, 4, 4);
     @test all(v .== ShiftedArray(v))
-    sv = ShiftedArray(v, (-2, 0))
+    sv = ShiftedArray(v, (-2, 0));
     @test length(sv) == 16
     @test sv[1, 3] == 11
     @test ismissing(sv[3, 3])
@@ -37,21 +39,23 @@ end
     @test isequal(sv, ShiftedArray(v, -2))
     @test isequal(@inferred(ShiftedArray(v, (2,))), @inferred(ShiftedArray(v, 2)))
     @test isequal(@inferred(ShiftedArray(v)), @inferred(ShiftedArray(v, (0, 0))))
-    s = ShiftedArray(v, (0, -2))
+    @test isequal(ShiftedArray(v, (2,)), ShiftedArray(v, 2))
+    @test isequal(ShiftedArray(v), ShiftedArray(v, (0, 0)))
+    s = ShiftedArray(v, (0, -2));
     @test isequal(collect(s), [ 9 13 missing missing;
                                10 14 missing missing;
                                11 15 missing missing;
                                12 16 missing missing])
-    sneg = ShiftedArray(v, (0, -2), default = -100)
+    sneg = ShiftedArray(v, (0, -2), default = -100);
     @test all(sneg .== coalesce.(s, default(sneg)))
     @test checkbounds(Bool, sv, 2, 2)
     @test !checkbounds(Bool, sv, 123, 123)
-    svnest = ShiftedArray(ShiftedArray(v, (1, 1)), 2)
-    sv = ShiftedArray(v, (3, 1))
+    svnest = ShiftedArray(ShiftedArray(v, (1, 1)), 2);
+    sv = ShiftedArray(v, (3, 1));
     @test sv === svnest
-    sv = ShiftedArray(v, 2, default = nothing)
-    sv1 = ShiftedArray(sv, (1, 1))
-    sv2 = ShiftedArray(sv, (1, 1), default = 0)
+    sv = ShiftedArray(v, 2, default = nothing);
+    sv1 = ShiftedArray(sv, (1, 1));
+    sv2 = ShiftedArray(sv, (1, 1), default = 0);
     @test isequal(collect(sv1), [nothing   nothing   nothing   nothing
                                  nothing   nothing   nothing   nothing
                                  nothing   nothing   nothing   nothing
@@ -63,7 +67,7 @@ end
 end
 
 @testset "padded_tuple" begin
-    v = rand(2, 2)
+    v = rand(2, 2);
     @test (1, 0) == @inferred ShiftedArrays.padded_tuple(v, 1)
     @test (0, 0) == @inferred ShiftedArrays.padded_tuple(v, ())
     @test (3, 0) == @inferred ShiftedArrays.padded_tuple(v, (3,))
@@ -100,7 +104,7 @@ end
     sv[3] = 12 
     @test collect(sv) == [3, 0, 12, 1]
     @test v == [1, 3, 0, 12]
-    @test sv === setindex!(sv, 12, 3) 
+    @test sv[3] === setindex!(sv, 12, 3)[3] # this was changed to adhere to the default return of setindex! which is an array and not a value
     @test checkbounds(Bool, sv, 2)
     @test !checkbounds(Bool, sv, 123)
     sv = CircShiftedArray(v, 3)
@@ -118,6 +122,8 @@ end
     @test isequal(sv, CircShiftedArray(v, -2))
     @test isequal(@inferred(CircShiftedArray(v, 2)), @inferred(CircShiftedArray(v, (2,))))
     @test isequal(@inferred(CircShiftedArray(v)), @inferred(CircShiftedArray(v, (0, 0))))
+    @test isequal(CircShiftedArray(v, 2), CircShiftedArray(v, (2,)))
+    @test isequal(CircShiftedArray(v), CircShiftedArray(v, (0, 0)))
     s = CircShiftedArray(v, (0, 2))
     @test isequal(collect(s), [ 9 13 1 5;
                                10 14 2 6;
@@ -128,6 +134,73 @@ end
     @test sv === svnest
 end
 
+@testset "ShiftedArray broadcast" begin
+    # Some tests for the broadcasting functionality
+    function test_broadcast(expr, src, sv)
+        ca = circshift(src, sv);
+        csa = CircShiftedArray(src, sv);
+        @test eltype(ca) == eltype(csa)
+        @test ca == csa
+        # approx is needed since the summing order is slightly different in both cases
+        @test sum(ca) ≈ sum(csa)
+        for d = 1:ndims(src)
+            # approx causes individual getindex:
+            @test sum(ca, dims=d) ≈ sum(csa, dims=d)
+        end
+        res = expr.(ca)
+        res_c = expr.(csa)
+        @test res_c == res
+        for d = 1:ndims(src)
+            # approx causes individual getindex:
+            @test sum(expr.(ca), dims=d) ≈ sum(expr.(csa), dims=d)
+        end
+        res_c = expr.(csa) .+ src
+        @test res_c == res .+ src
+        for d = 1:ndims(src)
+            @test sum(expr.(csa) .+ src, dims=d) ≈ sum(res .+ src, dims=d)
+        end
+        res_c = expr.(csa)
+        res1 = copy(src)
+        res1 .= expr.(ca)
+        @test res_c == res1 
+        res1 .= expr.(ca) .+ src
+        @test collect(res_c) .+ src == collect(res1)
+        res2 = CircShiftedArray(res1, sv)
+        res2 .= expr.(ca)
+        @test res2 == res_c
+        @test collect(res2) == collect(res_c)
+        @test sum(res2) == sum(res_c)
+        # test some shifted arrays
+    end
+    v = reshape(1:16, 4, 4)
+    test_broadcast(x->x+1,v,(2,-1))
+    sv = CircShiftedArray(v,(3,2))
+    @test collect(sv)[1:2,1:2] == sv[1:2,1:2]
+    @test sv[1:2,1:2] == @view sv[1:2,1:2]
+    v = rand(Int,3,4,5)
+    test_broadcast(x->x+1,v,(2,0,3))
+    v = rand(Int,6,4,8)
+    # test whether views are handled correctly
+    test_broadcast(x->x+1,(@view v[1:2,1:2,4:6]),(2,0,3))
+    v = rand(ComplexF32,5,4,3)
+    test_broadcast(x->x+2+x*x,v,(2,-1,3))
+
+    v = reshape(1:16, 4, 4);
+    sv = ShiftedArray(v, (-2, 1));
+    bcv = sv .+ v
+    sv2 = ShiftedArray(v, (-1, 1));
+    # mixing a shifted an not-shifted vector yields no shifted vector
+    @test !isa(bcv, ShiftedArray)
+    @test isa(sv .+ sv, ShiftedArray)
+
+    # this test is disabled, since Julia 1.5 has trouble with the syntax in @test_throws
+    #@test_throws "shifts ((-2, 1) and (-1, 1)) of both arrays need to be equal." (isa(sv .+ sv2, ShiftedArray))
+
+    ref = [missing 8 16 24; missing 10 18 26; missing missing missing missing; missing missing missing missing]
+    @test isequal(bcv, ref)
+    @test coalesce.(sv .+ v, 100) == coalesce.(ref, 100)
+end
+
 @testset "circshift" begin
     v = reshape(1:16, 4, 4)
     @test all(circshift(v, (1, -1)) .== ShiftedArrays.circshift(v, (1, -1)))
@@ -138,6 +211,17 @@ end
     @test sv === svnest
 end
 
+if VERSION >= v"1.6"
+    @testset "reverse" begin
+        v = collect(reshape(1:16, 4, 4))
+        cs = ShiftedArrays.circshift(v, (1,1))
+        q = reverse(cs)
+        reverse!(cs)
+        @test collect(q) == collect(cs)
+        @test q[1:1,1:1] == [11;;]
+    end
+end
+
 @testset "fftshift and ifftshift" begin
     function test_fftshift(x, dims=1:ndims(x))
         @test fftshift(x, dims) == ShiftedArrays.fftshift(x, dims)
@@ -159,6 +243,10 @@ end
 
     @test (2, 2, 0) == ShiftedArrays.ft_center_diff((4, 5, 6), (1, 2)) # Fourier center is at (2, 3, 0)
     @test (2, 2, 3) == ShiftedArrays.ft_center_diff((4, 5, 6), (1, 2, 3)) # Fourier center is at (2, 3, 4)
+    # ensure that circ-shifts with zero are converted back to ordinary arrays
+    @test isa(ShiftedArrays.ifftshift(ShiftedArrays.fftshift(randn(10,11,12))), Array)
+    @test isa(ShiftedArrays.fftshift(ShiftedArrays.ifftshift(randn(11,12,10))), Array)
+
 end
 
 @testset "laglead" begin