Ah okay, nevermind, I can reproduce the discrepancy you’re seeing by introducing a contracted index:

```
julia> let
N, M, P = 2,4, 3
A = rand(N, M)
B = rand(M, P)
C = rand(N)
@tullio out1[n, p] := A[n, m] * B[m, p] + C[n]
@tullio out2[n, p] := A[n, m] * B[m, p] + C[n] * ones(P)[p]
@tensor out3[n, p] := A[n, m] * B[m, p] + C[n] * ones(P)[p]
@info "" out1 out2 out3
end
┌ Info:
│ out1 =
│ 2×3 Matrix{Float64}:
│ 2.68322 2.8603 2.62903
│ 4.80266 4.25346 4.36135
│ out2 =
│ 2×3 Matrix{Float64}:
│ 2.68322 2.8603 2.62903
│ 4.80266 4.25346 4.36135
│ out3 =
│ 2×3 Matrix{Float64}:
│ 1.31961 1.49669 1.26541
└ 2.12742 1.57822 1.68612
```

It seems like what’s happening here is that instead of calculating

```
map(Iterators.product(1:N, 1:P)) do (n, p)
sum(1:M) do m
A[n,m] * B[m,p]
end + C[n]
end
```

Tullio is instead calculating

```
map(Iterators.product(1:N, 1:P)) do (n, p)
sum(1:M) do m
A[n,m] * B[m,p] + C[n]
end
end
```

That’s rather unexpected, I didn’t know it did this. cc @mcabbott

I guess this behaviour makes some sense, since Tullio is meant to represent general reductions, but there should maybe be some warnings about this in the README.