# Arithmetic performance of expression

**URL:** <https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196>\
**Category:** Performance\
**Created:** [October 3, 2022, 11:31pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196 "2022-10-03T23:31:10Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jojo/32/42799_2.png) [@Jojo](https://discourse.julialang.org/u/Jojo)\
**Post date:** [October 3, 2022, 11:31pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/1 "2022-10-03T23:31:10Z")

</div>

Hello!,

Recently I’ve noticed that Julia performs better when I don’t associate the terms in a expression im using in a for loop. Why is (1), faster than (2)? According to @btime (2) is 0.3 seconds slower than (1) after I double loop in a big mesh. This is probably basic but I’m trying to get better at producing high performance code and the benchmarks go against my intuition.

(1) `(x/11.0)*(a+b+c)+(d+e)*y/11.0+f*z/11.0`  
(2) `(x*(a+b+c)+(d+e)*y+f*z)/11.0`

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 4, 2022, 1:22am UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/2 "2022-10-04T01:22:15Z")

</div>

Have a reproducible example?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 4, 2022, 4:45am UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/3 "2022-10-04T04:45:01Z")

</div>

yeah. 2 really should be faster here. it would be interesting to see the benchmark.

---

<div class="post-metadata">

**Author:** ![Jojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jojo/32/42799_2.png) [@Jojo](https://discourse.julialang.org/u/Jojo)\
**Post date:** [October 4, 2022, 6:27pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/4 "2022-10-04T18:27:34Z")

</div>

@Oscar_Smith, @Elrod Couldn’t set up a minimal working example yesterday but here’s an example where I consistenly get higher btimes for 2

````julia
function foo1()
   dt=0.125
   C=0.3
   C2=C^2
   Nt=40000
   Nx=12000
   m=rand(Nx,Nt+1)
for j = 4:Nt
   @tturbo for i=3:Nx-2
      m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]) 
   end
end
end

function foo2()
   dt=0.125
   C=0.3
   C2=C^2
   Nt=40000
   Nx=12000
   m=rand(Nx,Nt+1)
for j = 4:Nt
   @tturbo for i=3:Nx-2
      m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0 
   end
  end
end

@btime foo1()
@btime foo2()
```
````

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [October 4, 2022, 6:56pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/5 "2022-10-04T18:56:05Z")

</div>

I removed `@tturbo` and I got

```julia
  3.721 s (2 allocations: 3.58 GiB)
  3.708 s (2 allocations: 3.58 GiB)

```

Also, at least for me half of the time is spent in allocating the big arrays

---

<div class="post-metadata">

**Author:** ![Jojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jojo/32/42799_2.png) [@Jojo](https://discourse.julialang.org/u/Jojo)\
**Post date:** [October 4, 2022, 7:04pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/6 "2022-10-04T19:04:37Z")

</div>

That’s interesting, it might be related to how @tturbo works. I also get better times for foo1() when removing the macro.

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [October 4, 2022, 7:11pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/7 "2022-10-04T19:11:42Z")

</div>

I’d also suggest benchmarking something like

```julia
using BenchmarkTools

function foo1!(m, dt, C, C2, Nt, Nx)
    for j = 4:Nt
        for i=3:Nx-2
            m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]) 
        end
    end
end

function foo2!(m, dt, C, C2, Nt, Nx)
    for j = 4:Nt
        for i=3:Nx-2
            m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0 
        end
    end
end

@btime foo1() setup = (
    dt=0.125;
    C=0.3;
    C2=C^2;
    Nt=40000;
    Nx=12000;
    m=rand(Nx,Nt+1);
)

@btime foo2() setup = (
    dt=0.125;
    C=0.3;
    C2=C^2;
    Nt=40000;
    Nx=12000;
    m=rand(Nx,Nt+1);
)

```

to remove allocating the big arrays from the run, they create lots of noise.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 4, 2022, 7:11pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/8 "2022-10-04T19:11:53Z")

</div>

Comparing `@tturbo` and `@inbounds @fastmath` versions:

```julia
julia> GC.gc(); @time foo1()
  1.922473 seconds (2 allocations: 3.576 GiB, 0.12% gc time)

julia> GC.gc(); @time foo1_fm()
  2.182295 seconds (2 allocations: 3.576 GiB, 0.04% gc time)

julia> GC.gc(); @time foo2()
  1.920486 seconds (2 allocations: 3.576 GiB, 0.05% gc time)

julia> GC.gc(); @time foo2_fm()
  2.186049 seconds (2 allocations: 3.576 GiB, 0.04% gc time)

```

I don’t see any difference between `foo1` and `foo2`.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [October 4, 2022, 7:13pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/9 "2022-10-04T19:13:49Z")

</div>

> [@Oscar\_Smith](#):
>
> yeah. 2 really should be faster here. it would be interesting to see the benchmark.

Why? Out of curiosity. Fewer divisions?

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 4, 2022, 7:21pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/10 "2022-10-04T19:21:12Z")

</div>

> **Modification of @giordano 's code**
>
> ```julia
> julia> function foo1!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1! (generic function with 1 method)
> 
> julia> function foo1_fm!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @inbounds @fastmath for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1_fm! (generic function with 1 method)
> 
> julia> function foo2!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0
> end
> end
> end
> foo2! (generic function with 1 method)
> 
> julia> function foo2_fm!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @inbounds @fastmath for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0
> end
> end
> end
> foo2_fm! (generic function with 1 method)
> 
> julia> using LoopVectorization
> 
> julia> function foo1!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @turbo for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1! (generic function with 1 method)
> 
> julia> function foo1!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1! (generic function with 1 method)
> 
> julia> function foo1_fm!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @inbounds @fastmath for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1_fm! (generic function with 1 method)
> 
> julia> function foo2!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0
> end
> end
> end
> foo2! (generic function with 1 method)
> 
> julia> function foo2_fm!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @inbounds @fastmath for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0
> end
> end
> end
> foo2_fm! (generic function with 1 method)
> 
> julia> using LoopVectorization
> 
> julia> function foo1_turbo!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @turbo for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1_turbo! (generic function with 1 method)
> 
> julia> function foo2_turbo!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @turbo for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0
> end
> end
> end
> foo2_turbo! (generic function with 1 method)
> 
> julia> function foo1_inbounds!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @inbounds for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3])/11.0-(12.0/11.0)*dt^2+(C2/11.0)*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j])
> end
> end
> end
> foo1_inbounds! (generic function with 1 method)
> 
> julia> function foo2_inbounds!(m, dt, C, C2, Nt, Nx)
> for j = 4:Nt
> @inbounds for i=3:Nx-2
> m[i,j+1] = (20.0*m[i,j]-6.0*m[i,j-1]-4.0*m[i,j-2]+m[i,j-3]-12.0*dt^2+C2*(16.0*(m[i-1,j]+m[i+1,j])-30.0*m[i,j]-m[i+2,j]-m[i-2,j]))/11.0
> end
> end
> end
> foo2_inbounds! (generic function with 1 method)
> 
> ```

results

```julia
julia> @btime foo1!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  1.586 s (0 allocations: 0 bytes)

julia> @btime foo2!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  1.767 s (0 allocations: 0 bytes)

julia> @btime foo1_inbounds!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  445.422 ms (0 allocations: 0 bytes)

julia> @btime foo2_inbounds!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  445.125 ms (0 allocations: 0 bytes)

julia> @btime foo1_fm!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  388.491 ms (0 allocations: 0 bytes)

julia> @btime foo2_fm!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  387.870 ms (0 allocations: 0 bytes)

julia> @btime foo1_turbo!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  360.542 ms (0 allocations: 0 bytes)

julia> @btime foo2_turbo!(m, dt, C, C2, Nt, Nx) setup = (
           dt=0.125;
           C=0.3;
           C2=C^2;
           Nt=40000;
           Nx=12000;
           m=rand(Nx,Nt+1);
       )
  355.190 ms (0 allocations: 0 bytes)

```

I don’t see much of a difference, except without any annotations

---

<div class="post-metadata">

**Author:** ![Jojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jojo/32/42799_2.png) [@Jojo](https://discourse.julialang.org/u/Jojo)\
**Post date:** [October 4, 2022, 7:36pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/12 "2022-10-04T19:36:41Z")

</div>

@giordano @Elrod I’ve run the same code you posted in my computer and somehow I get the following times for the turbo versions. It’s not such a huge difference but I’m puzzled about how the foo2\_turbo performs a little worst in my computer. Could it be related to the fact that I’m using an M1 processor?

(Aside question, while tweaking my code I noticed that it performs better when I don’t use the \* multiplication symbol just as you did, is that normal?)

```julia
296.276 ms (0 allocations: 0 bytes)
319.983 ms (0 allocations: 0 bytes)
```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 4, 2022, 8:01pm UTC](https://discourse.julialang.org/t/arithmetic-performance-of-expression/88196/13 "2022-10-04T20:01:56Z")

</div>

Using much smaller sizes

```julia
julia> Nt = 400; Nx = 120;

julia> @btime foo1!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  165.522 μs (0 allocations: 0 bytes)

julia> @btime foo2!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  184.462 μs (0 allocations: 0 bytes)

julia> @btime foo1_inbounds!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  42.759 μs (0 allocations: 0 bytes)

julia> @btime foo2_inbounds!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  45.094 μs (0 allocations: 0 bytes)

julia> @btime foo1_fm!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  27.266 μs (0 allocations: 0 bytes)

julia> @btime foo2_fm!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  26.341 μs (0 allocations: 0 bytes)

julia> @btime foo1_turbo!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  18.326 μs (0 allocations: 0 bytes)

julia> @btime foo2_turbo!(m, dt, C, C2, Nt, Nx) setup = (dt=0.125; C=0.3; C2=C^2; Nt=$Nt; Nx=$Nx; m=rand(Nx,Nt+1);)
  19.072 μs (0 allocations: 0 bytes)

```

FWIW, assembly of the `@turbo` versions, foo1:

```nohighlight
.LBB0_3: # %L84
                                        # Parent Loop BB0_2 Depth=1
                                        # => This Inner Loop Header: Depth=2
        vmovupd zmm10, zmmword ptr [rbp + 8*r15 - 152]
        vmovupd zmm11, zmmword ptr [rbp + 8*r15 - 144]
        vmovupd zmm12, zmmword ptr [rbp + 8*r15 - 128]
        vmovupd zmm13, zmmword ptr [rbp + 8*r15 - 88]
        vmovupd zmm14, zmmword ptr [rbp + 8*r15 - 24]
        vaddpd zmm10, zmm10, zmmword ptr [rbp + 8*r15 - 136]
        vaddpd zmm13, zmm13, zmmword ptr [rbp + 8*r15 - 72]
        vaddpd zmm14, zmm14, zmmword ptr [rbp + 8*r15 - 8]
        vmovupd zmm15, zmmword ptr [rbp + 8*r15 - 80]
        vmovupd zmm16, zmmword ptr [rbp + 8*r15 - 16]
        vmovupd zmm17, zmmword ptr [rbp + 8*r15 - 64]
        vmovupd zmm18, zmmword ptr [rbp + 8*r15]
        vaddpd zmm12, zmm12, zmmword ptr [rbp + 8*r15 - 160]
        vaddpd zmm17, zmm17, zmmword ptr [rbp + 8*r15 - 96]
        vaddpd zmm18, zmm18, zmmword ptr [rbp + 8*r15 - 32]
        vfmadd231pd zmm12, zmm11, zmm2 # zmm12 = (zmm11 * zmm2) + zmm12
        vfmadd231pd zmm17, zmm15, zmm2 # zmm17 = (zmm15 * zmm2) + zmm17
        vfmadd231pd zmm18, zmm16, zmm2 # zmm18 = (zmm16 * zmm2) + zmm18
        vfnmadd231pd zmm12, zmm3, zmm10 # zmm12 = -(zmm3 * zmm10) + zmm12
        vfnmadd231pd zmm17, zmm3, zmm13 # zmm17 = -(zmm3 * zmm13) + zmm17
        vfnmadd231pd zmm18, zmm3, zmm14 # zmm18 = -(zmm3 * zmm14) + zmm18
        vmovupd zmm10, zmmword ptr [rcx + 8*r15 - 128]
        vmovupd zmm13, zmmword ptr [rcx + 8*r15 - 64]
        vmovupd zmm14, zmmword ptr [rcx + 8*r15]
        vfmadd213pd zmm10, zmm4, zmmword ptr [r13 + 8*r15 + 16] # zmm10 = (zmm4 * zmm10) + mem
        vfmadd213pd zmm13, zmm4, zmmword ptr [r13 + 8*r15 + 80] # zmm13 = (zmm4 * zmm13) + mem
        vfmadd213pd zmm14, zmm4, zmmword ptr [r13 + 8*r15 + 144] # zmm14 = (zmm4 * zmm14) + mem
        vfmadd231pd zmm10, zmm5, zmmword ptr [rbx + 8*r15 - 128] # zmm10 = (zmm5 * mem) + zmm10
        vfmadd231pd zmm13, zmm5, zmmword ptr [rbx + 8*r15 - 64] # zmm13 = (zmm5 * mem) + zmm13
        vfmadd231pd zmm14, zmm5, zmmword ptr [rbx + 8*r15] # zmm14 = (zmm5 * mem) + zmm14
        vfmadd231pd zmm10, zmm6, zmm11 # zmm10 = (zmm6 * zmm11) + zmm10
        vfmadd231pd zmm13, zmm6, zmm15 # zmm13 = (zmm6 * zmm15) + zmm13
        vfmadd231pd zmm14, zmm6, zmm16 # zmm14 = (zmm6 * zmm16) + zmm14
        vmulpd zmm10, zmm10, zmm9
        vmulpd zmm11, zmm13, zmm7
        vmulpd zmm13, zmm14, zmm7
        vfmsub231pd zmm10, zmm0, zmm12 # zmm10 = (zmm0 * zmm12) - zmm10
        vfmsub231pd zmm11, zmm0, zmm17 # zmm11 = (zmm0 * zmm17) - zmm11
        vfmsub231pd zmm13, zmm0, zmm18 # zmm13 = (zmm0 * zmm18) - zmm13
        vfmsub231pd zmm10, zmm1, zmm8 # zmm10 = (zmm1 * zmm8) - zmm10
        vfmsub231pd zmm11, zmm1, zmm8 # zmm11 = (zmm1 * zmm8) - zmm11
        vfmsub231pd zmm13, zmm1, zmm8 # zmm13 = (zmm1 * zmm8) - zmm13
        vmovupd zmmword ptr [rax + 8*r15 - 128], zmm10
        vmovupd zmmword ptr [rax + 8*r15 - 64], zmm11
        vmovupd zmmword ptr [rax + 8*r15], zmm13
        add r15, 24
        cmp r15, rsi
        jle .LBB0_3

```

foo2:

```nohighlight
.LBB0_3: # %L82
                                        # Parent Loop BB0_2 Depth=1
                                        # => This Inner Loop Header: Depth=2
        vmovupd zmm12, zmmword ptr [rcx + 8*r15 - 80]
        vmovupd zmm13, zmmword ptr [rcx + 8*r15 - 16]
        vmovupd zmm14, zmmword ptr [rcx + 8*r15 - 152]
        vmovupd zmm15, zmmword ptr [rcx + 8*r15 - 144]
        vmovupd zmm16, zmmword ptr [rcx + 8*r15 - 128]
        vmovupd zmm17, zmmword ptr [rcx + 8*r15 - 88]
        vmovupd zmm18, zmmword ptr [rcx + 8*r15 - 24]
        vaddpd zmm14, zmm14, zmmword ptr [rcx + 8*r15 - 136]
        vaddpd zmm17, zmm17, zmmword ptr [rcx + 8*r15 - 72]
        vaddpd zmm18, zmm18, zmmword ptr [rcx + 8*r15 - 8]
        vmovupd zmm19, zmmword ptr [rcx + 8*r15 - 64]
        vmovupd zmm20, zmmword ptr [rcx + 8*r15]
        vaddpd zmm16, zmm16, zmmword ptr [rcx + 8*r15 - 160]
        vaddpd zmm19, zmm19, zmmword ptr [rcx + 8*r15 - 96]
        vaddpd zmm20, zmm20, zmmword ptr [rcx + 8*r15 - 32]
        vfmadd231pd zmm16, zmm15, zmm2 # zmm16 = (zmm15 * zmm2) + zmm16
        vfmadd231pd zmm19, zmm12, zmm2 # zmm19 = (zmm12 * zmm2) + zmm19
        vfmadd231pd zmm20, zmm13, zmm2 # zmm20 = (zmm13 * zmm2) + zmm20
        vfnmadd231pd zmm16, zmm3, zmm14 # zmm16 = -(zmm3 * zmm14) + zmm16
        vfnmadd231pd zmm19, zmm3, zmm17 # zmm19 = -(zmm3 * zmm17) + zmm19
        vfnmadd231pd zmm20, zmm3, zmm18 # zmm20 = -(zmm3 * zmm18) + zmm20
        vfmsub213pd zmm16, zmm0, zmmword ptr [r13 + 8*r15 + 16] # zmm16 = (zmm0 * zmm16) - mem
        vfmsub213pd zmm19, zmm0, zmmword ptr [r13 + 8*r15 + 80] # zmm19 = (zmm0 * zmm19) - mem
        vfmsub213pd zmm20, zmm0, zmmword ptr [r13 + 8*r15 + 144] # zmm20 = (zmm0 * zmm20) - mem
        vfnmadd231pd zmm16, zmm1, zmm4 # zmm16 = -(zmm1 * zmm4) + zmm16
        vfnmadd231pd zmm19, zmm1, zmm4 # zmm19 = -(zmm1 * zmm4) + zmm19
        vfnmadd231pd zmm20, zmm1, zmm4 # zmm20 = -(zmm1 * zmm4) + zmm20
        vfnmadd231pd zmm16, zmm5, zmmword ptr [rbp + 8*r15 - 128] # zmm16 = -(zmm5 * mem) + zmm16
        vfnmadd231pd zmm19, zmm5, zmmword ptr [rbp + 8*r15 - 64] # zmm19 = -(zmm5 * mem) + zmm19
        vfnmadd231pd zmm20, zmm5, zmmword ptr [rbp + 8*r15] # zmm20 = -(zmm5 * mem) + zmm20
        vfnmadd231pd zmm16, zmm6, zmmword ptr [rbx + 8*r15 - 128] # zmm16 = -(zmm6 * mem) + zmm16
        vfnmadd231pd zmm19, zmm6, zmmword ptr [rbx + 8*r15 - 64] # zmm19 = -(zmm6 * mem) + zmm19
        vfnmadd231pd zmm20, zmm6, zmmword ptr [rbx + 8*r15] # zmm20 = -(zmm6 * mem) + zmm20
        vfmsub231pd zmm16, zmm7, zmm15 # zmm16 = (zmm7 * zmm15) - zmm16
        vmulpd zmm14, zmm16, zmm11
        vfmsub231pd zmm19, zmm7, zmm12 # zmm19 = (zmm7 * zmm12) - zmm19
        vfmsub231pd zmm20, zmm7, zmm13 # zmm20 = (zmm7 * zmm13) - zmm20
        vmulpd zmm12, zmm19, zmm9
        vmulpd zmm13, zmm20, zmm9
        vmovupd zmmword ptr [rax + 8*r15 - 128], zmm14
        vmovupd zmmword ptr [rax + 8*r15 - 64], zmm12
        vmovupd zmmword ptr [rax + 8*r15], zmm13
        add r15, 24
        cmp r15, rsi
        jle .LBB0_3

```

Assembly doesn’t look strikingly different. No divisions, it’s been replaced with multiplication by the inverse.  
Both versions require 10 arithmetic instructions, a mix of adds, fmas, and multiplies.

I have an M1 at home, so I can take a look later.