In the following code

```
import Base: ^
struct cccbbb
v::Matrix{Float64}
end
^(cb::cccbbb, p::Int) = cb.v^p
Base.Broadcast.broadcastable(cb::cccbbb) = Ref(cb)
n = 5
r = LinRange(0.0, 10.0, 5)
cb = cccbbb(rand(n,n))
tmp = cb.^2
@show typeof(cb.^2), typeof(cb.^2 ./r), typeof(tmp./r)
```

the result is

```
(typeof(cb .^ 2), typeof(cb .^ 2 ./ r), typeof(tmp ./ r)) = (Matrix{Float64}, Vector{Matrix{Float64}}, Matrix{Float64})
```

The type of (cb .^ 2 ./ r) is not as expected. Is this a bug?

Yeah, that’s a little confusing, but I think the result is correct. Here’s a slightly simplified version of your example:

```
import Base: ^
struct A
m::Matrix{Float64}
end
^(a::A, p::Int) = a.m ^ p
Base.broadcastable(a::A) = Ref(a)
m = [1 2
3 4]
a = A(m)
```

```
julia> m ^ 2
2×2 Array{Int64,2}:
7 10
15 22
julia> m .^ 2
2×2 Array{Int64,2}:
1 4
9 16
julia> m .^ 2 ./ [1, 2]
2×2 Array{Float64,2}:
1.0 4.0
4.5 8.0
julia> a .^ 2 ./ [1, 2]
2-element Array{Array{Float64,2},1}:
[7.0 10.0; 15.0 22.0]
[3.5 5.0; 7.5 11.0]
```

The `Base.broadcastable(a::A) = Ref(a)`

line tells Julia to treat an object of type `A`

as a scalar for broadcasting. So

```
a .^ 2 ./ [1, 2]
```

is effectively equivalent to

```
[(a^2)/1, (a^2)/2]
```

And they way you’ve defined `^`

on your new type means that it is doing a matrix power (`m * m`

), *not* an element by element power.

Can you elaborate on what output you expected and/or what output you are trying to achieve?

3 Likes

As in your example,

```
julia> b = a.^2
2×2 Matrix{Float64}:
7.0 10.0
15.0 22.0
```

and

```
julia> b ./ [ 1, 2]
2×2 Matrix{Float64}:
7.0 10.0
7.5 11.0
```

The result is the Matrix{Float64} type as expected. But the

```
a .^ 2 ./ [1, 2]
```

type is

```
2-element Array{Array{Float64,2},1}
```

which is not expected.

Can the `Base.broadcastable(a::A) = Ref(a)`

just affect the operator `^`

?

What we’re trying to tell you is that `b ./ [1, 2]`

and `a .^ 2 ./ [1, 2]`

are two distinct computations. An important feature of the dot syntax is that vectorized functions are fused into a single loop. In the second case, you’re thinking Julia would first broadcast exponentiation to obtain the matrix, but as @CameronBieganek noted in his reply, the operations are fused into `[(a^2)/1, (a^2/2)]`

. This is the expected behavior, since element-wise exponentiation *is* different from exponentiation after all (take, matrices, for example). The behavior here is more confusing than usual because you’re wrapping a matrix, but also tell Julia to treat your struct as a scalar.

If you want the same result with the first case, try `a ^ 2 ./ [1,2]`

instead.

3 Likes