Search code examples
recursionviewjuliavectorization

Fill array recursively in Julia using dot-notation and views


I am using Julia 1.9.0 and I’m trying to fill an array using views and dot notation where each entry depends on the preceding one(s). Say, I want the first 10 entries of the Fibonacci series in my array, so I initialize the first two entries and want the rest computed without a for-loop.

My approach was:

v = zeros(Int64, 10)
v[1:2] .= 1

@views v[3:10] .= v[1:8] .+ v[2:9]
v

Output:

10-element Vector{Int64}:
 1
 1
 2
 1
 0
 0
 0
 0
 0
 0

This doesn’t work, though, because apparently v doesn’t get updated before the next value is calculated. From reading the docs, I couldn’t find out why, because as far as I understood no arrays are copied and no temporary array is created. This should be equivalent to:

v = zeros(Int64, 10)
v[1:2] .= 1

@views a, b, c = v[1:8], v[2:9], v[3:10]
for i in 1:8
  c[i] = a[i] + b[i]
end
v

Output:

10-element Vector{Int64}:
  1
  1
  2
  3
  5
  8
 13
 21
 34
 55

Which is exactly the result I expected to get. Why doesn’t the first version work and is it possible to solve this without the for-loop?


Solution

  • Broadcasting does alias detection. Therefore actually your operation copies data, because you have aliases.

    To convince yourself that this is what happens check this code:

    julia> v = zeros(Int64, 10^7);
    
    julia> v[1:2] .= 1;
    
    julia> @allocated v[3:end] .= @views v[1:end-2] .+ v[2:end-1]
    

    This check is done using the Base.mightalias function.

    Dealiasing is done to protect RHS to be affected by mutation of the LHS. So, on a basic level, the answer is that with broadcasting you cannot implement recursion.

    However, you could create a custom type that would always return false in Base.mightalias aliasing check in which case what you ask for would work.

    Here is an example:

    julia> using SentinelArrays
    
    julia> v = zeros(Int64, 10^7);
    
    julia> v[1:2] .= 1;
    
    julia> @allocated v[3:end] .= ChainedVector([@view v[1:end-2]]) .+ ChainedVector([@view v[2:end-1]])
    800
    
    julia> v
    10000000-element Vector{Int64}:
                        1
                        1
                        2
                        3
                        5
                        8
                       13
                       21
                       34
                       55
                       89
                      144
                      233
                      377
                      610
                      987
                        ⋮