# Performance challenge: can you write a faster sum?

**URL:** <https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456>\
**Category:** Performance\
**Tags:** simd\
**Created:** [July 3, 2025, 3:47pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456 "2025-07-03T15:47:21Z")\
**Posts on this page:** 13\
**Page:** 2

<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:** [July 7, 2025, 11:02pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/21 "2025-07-07T23:02:51Z")

</div>

> [@mbauman](#):
>
> ```julia
> julia> using BenchmarkTools
> 
> julia> A = rand(10000); B = @view A[begin:2:end];
> 
> julia> @btime sum($A);
> 946.640 ns (0 allocations: 0 bytes)
> 
> julia> @btime sum($B);
> 4.149 μs (0 allocations: 0 bytes)
> 
> ```

Results are of course platform-dependent; zen5 seems to prefer `@turbo`:

```julia
julia> @btime reduce_mb_3(+, $A); @btime reduce_mb_3(+, $B)
  440.970 ns (0 allocations: 0 bytes)
  1.053 μs (0 allocations: 0 bytes)
2475.6003f0

julia> @btime sum_turbo($A); @btime sum_turbo($B)
  58.546 ns (0 allocations: 0 bytes)
  748.756 ns (0 allocations: 0 bytes)
2475.6f0

julia> versioninfo()
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: 32 × AMD Ryzen 9 9950X 16-Core Processor

```

---

<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:** [July 7, 2025, 11:07pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/22 "2025-07-07T23:07:04Z")

</div>

> [@mbauman](#):
>
> That’s not possible for a general `@view A[1, :]` sort of access (think a reduction with `dims=2`)

LV seems to use the exact same order of operations here, at least at this size:

```julia
julia> sum(reshape(A,(100,100)),dims=2) == vreduce(+,reshape(A,(100,100)),dims=2)
true

julia> @benchmark sum(reshape($A,(100,100)),dims=2)
BenchmarkTools.Trial: 8993 samples with 198 evaluations per sample.
 Range (min … max): 443.495 ns … 10.497 μs ┊ GC (min … max): 0.00% … 94.60%
 Time (median): 482.333 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 557.728 ns ± 736.524 ns ┊ GC (mean ± σ): 12.81% ± 9.23%

     ▁█▂                                                         
  ▂▃▃███▆▆▆▅▅▅▅▅▅▅▅▆▅▅▅▅▄▅▅▄▄▄▄▅▅▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▁▁▁▁▁▁▂▁▁▁▂ ▃
  443 ns Histogram: frequency by time 619 ns <

 Memory estimate: 544 bytes, allocs estimate: 3.

julia> @benchmark vreduce(+,reshape($A,(100,100)),dims=2)
BenchmarkTools.Trial: 6451 samples with 767 evaluations per sample.
 Range (min … max): 127.382 ns … 3.704 μs ┊ GC (min … max): 0.00% … 95.52%
 Time (median): 141.150 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 200.908 ns ± 375.622 ns ┊ GC (mean ± σ): 25.97% ± 13.02%

  █                                                              
  █▄▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▂▂▁▁▁▁▂▂▂▂▂▂▂▂▂▂ ▂
  127 ns Histogram: frequency by time 2.81 μs <

 Memory estimate: 544 bytes, allocs estimate: 3.

```

Note that for `dims=2`, if your array is column major, you should SIMD across the rows.  
LV should just be using the naive summation order above, which is definitely not what you want for accuracy.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 8, 2025, 3:34pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/23 "2025-07-08T15:34:27Z")

</div>

Nice! I should’ve thought to check LV’s code. I can see that the generic `vreduce((x,y)->x+y, A)` uses a well-chosen [N-accumulators re-ordering strategy](https://github.com/JuliaSIMD/LoopVectorization.jl/blob/675b7b9b89426ed798dae5edb31b1537e2c480d9/src/simdfunctionals/mapreduce.jl#L86-L94), but it’s not as clear to me what the ordering and associativity is for the `::typeof(+)` specialization — those are defined with the `@turbo` macros over the loop and I’ve not tried to detangle that macro (or the generated code) further.

On my architecture for a `Float64` dense array, though, LV is behaving similarly (albeit not bit-exact) to my `reduce_mb_3`, which uses an interleaved 4x4 pattern — that is, it uses 4 accumulators, each summing four values at a time before reducing with the accumulator. That also seems to be doing quite well in general across multiple datatypes and layouts — including other iterators that are outside the scope of LV.

---

<div class="post-metadata">

**Author:** ![simsurace](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simsurace/32/30216_2.png) [@simsurace](https://discourse.julialang.org/u/simsurace)\
**Post date:** [July 8, 2025, 8:07pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/24 "2025-07-08T20:07:14Z")

</div>

I‘ve been wondering if pairwise reduction intent could be encoded via an iterator or a transducer so it could be wrapped around other iterators. But maybe that would make it harder for the compiler to optimize…

---

<div class="post-metadata">

**Author:** ![Gesee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gesee/32/217235_2.png) [@Gesee](https://discourse.julialang.org/u/Gesee)\
**Post date:** [July 8, 2025, 8:28pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/25 "2025-07-08T20:28:07Z")

</div>

It’s normal for sum `B` to be slower than sum `A`. `B` isn’t contiguous in memory, when iterating on him, we need to jump from an element `i` to an element `i+2` instead of `i+1` making it less faster an disabling advanced optimization like `@simd`. But if you find a way to make that negligeable, I would be glad.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [July 8, 2025, 10:47pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/26 "2025-07-08T22:47:01Z")

</div>

> [@simsurace](#):
>
> I‘ve been wondering if pairwise reduction intent could be encoded via an iterator or a transducer so it could be wrapped around other iterators. But maybe that would make it harder for the compiler to optimize…

Good question! I tried this with TransducersNext.jl and it worked pretty great:

```julia
using Base.Cartesian: @nexprs
using TransducersNext
using TransducersNext: Executor, complete, next, combine, @unroll
struct BlockedCommutativeEx <: Executor end

function TransducersNext. __fold__ (rf::RF, init::T, A, ::BlockedCommutativeEx) where {RF, T}
    op(x, y) = next(rf, x, y)
    inds = eachindex(A)
    i1, iN = firstindex(inds), lastindex(inds)
    n = length(inds)
    v = init
    # statically peel off the first iteration for type stability reasons if we don't have a numeric `init`
    @unroll 1 for batch in 0:(n>>4)-1
        i = i1 + batch*16
        @nexprs 16 N-> a_N = @inbounds A[inds[i+(N-1)]]
        @nexprs 2 N-> v_N = op(
            op(op(op(op(op(op(op(init,
                                 a_{(N-1)*8+1}),
                              a_{(N-1)*8+2}),
                           a_{(N-1)*8+3}),
                        a_{(N-1)*8+4}),
                     a_{(N-1)*8+5}),
                  a_{(N-1)*8+6}),
               a_{(N-1)*8+7}),
            a_{(N-1)*8+8})
        v = combine(rf, v, combine(rf, v_1, v_2))
    end
    i = i1 + (n>>4)*16 - 1
    i == iN && return v
    for i in i+1:iN
        ai = @inbounds A[inds[i]]
        v = op(v, ai)
    end
    return complete(rf, v)
end

```

```julia-repl
julia> @btime fold(+, $A; executor=BlockedCommutativeEx())
  1.146 μs (0 allocations: 0 bytes)
5003.618592521285

julia> @btime fold(+, $B; executor=BlockedCommutativeEx())
  767.919 ns (0 allocations: 0 bytes)
2509.5895307786136

```

Basically the same timings as with @mbauman’s function for `A`, but surprisingly it does better than Matt’s version for `B`:

```julia-repl
julia> @btime reduce_mb_4(+, $A)
  1.144 μs (0 allocations: 0 bytes)
5003.618592521285

julia> @btime reduce_mb_4(+, $B)
  1.110 μs (0 allocations: 0 bytes)
2509.5895307786136

```

What’s especially cool about doing this with TransducersNext is that you can do stuff like

```julia
executor = ChunkedEx(BlockedCommutativeEx; chunksize=1024)

```

or whatever to do further outer-chunking of the algorithm in a composable way, with minimal performance loss:

```julia
julia> @btime fold(+, $B; executor=ChunkedEx(BlockedCommutativeEx(); chunksize=1024), init=0)
  869.175 ns (0 allocations: 0 bytes)
2507.819648819429

```

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 9, 2025, 1:57pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/27 "2025-07-09T13:57:11Z")

</div>

That sure is a funny way of spelling `reduce` 🙂

I think of transducers as _changing the reduction operator_, not the “executor.” The executor is a well-structured reduction loop — exactly what I’m trying to define — and can equivalently be swapped with the function call itself. Transducers themselves aren’t (easily) able to rearrange their execution. I wouldn’t expect any difference of behavior, unless you have infrastructure to unwrap SubArrays (with a transducer!) or some such in some indirection before you get to ` __fold__ `. Also: are you missing an element there? Or did you use different data?

I do really like how you expose the ability to do something like `ChunkedEx(BlockedCommutativeEx(); chunksize=1024)`, that’s very cool, and I can immediately imagine exactly how that works.

I really think Transducers have a bit of the Monad problem. It’s a great and straightforward idea, wrapped up in an abstraction and language that makes it seem more complicated than it is… and then once you start shifting compute around your potential design space becomes intractably big.

Transducers fundamentally rely upon the _exact particulars_ of the executor; that’s something we’ve not well-defined for `Base.reduce` and is at the core of what I’m trying to do.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [July 9, 2025, 2:24pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/28 "2025-07-09T14:24:59Z")

</div>

> [@mbauman](#):
>
> That sure is a funny way of spelling `reduce` 🙂

`reduce` was already taken 🙂

> [@mbauman](#):
>
> I think of transducers as _changing the reduction operator_, not the “executor.” The executor is a well-structured reduction loop — exactly what I’m trying to define — and can equivalently be swapped with the function call itself. Transducers themselves aren’t (easily) able to rearrange their execution. I wouldn’t expect any difference of behavior, unless you have infrastructure to unwrap SubArrays (with a transducer!) or some such in some indirection before you get to ` __fold__ `.
> 
> […]
> 
> I do really like how you expose the ability to do something like `ChunkedEx(BlockedCommutativeEx(); chunksize=1024)`, that’s very cool, and I can immediately imagine exactly how that works.

Yeah, this is one of the main things that TransducersNext.jl is trying to tackle (and maybe it needs a rename). Yes, a transducer is only about the reduction operator. However, we _also_ have a concept of “executors” which express how the data itself is partitioned and fed into the reducing operator, and I want to bring this more into the foreground.

The idea is that I want to make it so that people writing Executors don’t need to worry about re-implementing transducers, and people writing custom transducers shouldn’t have to worry about re-implementing all the executors.

So in TransducersNext.jl, I have

- `SequentialEx` (this basically does `foldl`)
- `SIMDEx` for making a simple `@simd` loop
- `ChunkedEx` for chunking
- `ThreadEx` for multithreading partitioning
- `DistributedEx` for distributed partitioning
- `KernelAbstractionsEx` under development
- Maybe soon a `BlockedCommutativeEx` under development

All but sequential, SIMD, and BlockedCommutative take an inner executor so you can compose them.

> [@mbauman](#):
>
> Also: are you missing an element there? Or did you use different data?

Yes, there was a bug. Our messages crossed and I updated my comment right before you sent yours, sorry about that.

> [@mbauman](#):
>
> I really think Transducers have a bit of the Monad problem. It’s a great and straightforward idea, wrapped up in an abstraction and language that makes it seem more complicated than it is… and then once you start shifting compute around your potential design space becomes intractably big.

Totally agreed. I really want to find a way of simplifying this stuff down and making it less scary. I’ve had some success, but there’s a LONG way to go and it’ll require some very hard design and pedagogical decisions.

> [@mbauman](#):
>
> Transducers fundamentally rely upon the _exact particulars_ of the executor; that’s something we’ve not well-defined for `Base.reduce` and is at the core of what I’m trying to do.

Actually I disagree with this, they should be agnostic about the particulars of the executor other than the usual stuff (like certain executors requiring associativity or commutativity).

For instance, here’s how the `BlockedCommutativeEx` works totally fine (and quick!) with the `Map` transducer (and is faster than `mapreduce`!):

```julia
julia> @btime fold(+, Map(sin), $B; executor=BlockedCommutativeEx())
  12.533 μs (0 allocations: 0 bytes)
2306.4133570501153

julia> @btime mapreduce(sin, +, $B)
  21.440 μs (0 allocations: 0 bytes)
2306.413357050116

```

Here’s using it with `Map` and `Filter` (though with no perf benefit here vs traditional executors):

```julia
julia> C = rand(1:10, 10000);

julia> @btime fold(+, Filter(iseven) ⨟ Map(sin), $C; executor=SIMDEx())
  31.288 μs (0 allocations: 0 bytes)
274.85264475617697

julia> @btime mapreduce(sin, +, Iterators.filter(iseven, $C))
  31.268 μs (0 allocations: 0 bytes)
274.85264475617697

```

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [July 9, 2025, 2:32pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/29 "2025-07-09T14:32:13Z")

</div>

> [@mbauman](#):
>
> I wouldn’t expect any difference of behavior, unless you have infrastructure to unwrap SubArrays (with a transducer!) or some such in some indirection before you get to ` __fold__ `

Just to expand on this, if one is working with a transducer, it’s very important to properly handle when you use the innermost reducing operator (e.g. `+`) vs when you’re using the transformed operator (e.g. `Map(sin)'(+)`). One needs to account for this whenever doing re-associations, and that’s the source of the differences between my version and yours. Otherwise you end up applying `sin` multiple times when combining sub-reductions!

So there are subtle differences, but the general shape of the algorithm is the same.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 9, 2025, 2:39pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/30 "2025-07-09T14:39:41Z")

</div>

> [@Mason](#):
>
> Actually I disagree with this, they should be agnostic about the particulars of the executor other than the usual stuff (like certain executors requiring associativity or commutativity).

I’m actually not thinking about associativity here, but rather the things like:

- How `init` is used (or not!)
- How the empty- and one-element cases behave
- If `op` always receives the accumulator on the LHS
- If `op` might receive two elements directly, or if it always takes an accumulator (or init) and an element
- If `op` might receive two accumulators (or if it’s a potentially separate `combine` operation)
- If `op` might be speculatively executed or memoized
- … and yes, exactly as you just wrote “it’s very important to properly handle when you use the innermost reducing operator”

You can’t get any of those assumptions wrong or leave them to chance as you design your transformed `op`.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [July 9, 2025, 2:45pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/31 "2025-07-09T14:45:55Z")

</div>

Ah I see what you mean. Yeah, because Transducers / TransducersNext run into these issues way more ‘catastrophically’ than regular `reduce`, we take a much more rigid set of rules here that all executors must follow or chaos ensues:

- `init` is _always_ used on every (sub)reduction.
- one-element case is done with `init`. zero element case is an error if it hits the default `init`.
- `rf` (the transformed reducing operator e.g. `Map(sin)'(+)`) always takes the accumulator on the LHS. The inner-most reducing operator (e.g. `+`) can take two accumulators\*\*
- We don’t currently have a policy on speculative execution or memoization, but we probably should.

* * *

Edit: I should mention that we have an overloabable interface for how two accumulators are combined: `combine(rf, l, r)`. By default, this strips out all the transducers out of `rf` and goes down to the inner-most reducing operator, but each transducer in the chain and the innermost reducing operator are both allowed to intervene and change what `combine` means for them.

On the other hand, it’s `next(rf, acc, elem)` that always gets an accumulator on the left and an element of the iterator in on right.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 9, 2025, 5:03pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/32 "2025-07-09T17:03:23Z")

</div>

So, yeah, to bring this back to my original intent here: my overarching goal is to bring _some sort_ of a “better” _generic_ associativity (and _maybe_ a reordering for some ops if necessary) and rules to `Base.reduce` in a way that can move [WIP: The great pairwise reduction refactor by mbauman · Pull Request #58418 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/58418) forward.

My biggest question is if there _exists_ a decent-enough generic “SIMD-oblivious algorithm” (which isn’t really a thing to my knowledge) that explicitly encodes an associativity that allows for _sufficient_ use of _some_ SIMD instructions (if possible, on most platforms) to avoid major performance regressions from the status-quo `@simd` loops. There are many subtle reasons that have all added up to my distaste for `@simd` for something as fundamental as `reduce`. Things that can change its results include (but are not limited to): `--bounds-check=yes`, inlineability, type instabilities, views, transposes (and different `dims` thereof), generator vs. comprehension, iterator style, fma fusion, etc., etc., etc., etc.

The wonderful part is that batched pairwise and inner SIMD-like associativities tend to improve performance on modern CPUs _even if they don’t use SIMD instructions_ thanks to their reduced data dependencies. If `op` itself is much more expensive than a single instruction or two, then that should dominate any sort of overhead we incur in making this more complicated than the straight loop. And by explicitly batching up the extraction (getindex/iterate) of N elements at a time, N applications of `f`, and then the applications of `op`, we can often _recover_ the use of SIMD instructions in cases that wouldn’t have used them with `@simd` in the first place.

If we can get the pairwise reduction refactor in, it’s going to cause some churn (because folks depend on its internals, bad usages of `init`, exact values, etc). So all the more reason to make sure its new general behaviors are working as best as they can be. Maybe it makes sense to even explicitly document its associativity — numpy documents theirs.

---

<div class="post-metadata">

**Author:** ![epilliat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/epilliat/32/219946_2.png) [@epilliat](https://discourse.julialang.org/u/epilliat)\
**Post date:** [December 11, 2025, 10:03pm UTC](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456/33 "2025-12-11T22:03:55Z")

</div>

Hello, I just tried your sum from the PR, and I’ve generalized your ideas using `@generated` functions: implementing accumulators (for your `Base.jl` and `reduce_mb_4` versions), and incorporating tree reduction within the loop itself (`reduce_mb_3`).

Here is my attempt:

```julia
macro tree_reduce(op, prefix, N, M, stride)
    prefix_str = string(prefix)

    function build_tree(indices)
        count = length(indices)
        if count == 1
            return Symbol(prefix_str, :_, indices[1])
        elseif count == 2
            return :($op($(Symbol(prefix_str, :_, indices[1])),
                $(Symbol(prefix_str, :_, indices[2]))))
        else
            mid = count ÷ 2
            left = build_tree(indices[1:mid])
            right = build_tree(indices[mid+1:end])
            return :($op($left, $right))
        end
    end

    # Generate strided indices: N, N+stride, N+2*stride, ..., N+(M-1)*stride
    indices = [N + i * stride for i in 0:M-1]

    return esc(build_tree(indices))
end

@generated function reduce_abstracted(op, A, ::Val{Ntot}, ::Val{Nacc}) where {Ntot,Nacc}
    Pacc = trailing_zeros(Nacc)
    Ptot = trailing_zeros(Ntot)
    quote
        f = identity
        inds = eachindex(A)
        i1, iN = firstindex(inds), lastindex(inds)
        n = length(inds)
        @nexprs $Ntot N -> a_N = @inbounds A[inds[i1+(N-1)]]
        @nexprs $Nacc N -> v_N = @tree_reduce(op, a, N, $(cld(Ntot, Nacc)), $Nacc)

        for batch in 1:(n>>$Ptot)-1
            i = 1 + (batch << $Ptot)
            @nexprs $Ntot N -> begin
                a_N = @inbounds A[inds[i+(N-1)]]
            end
            @nexprs $Nacc N -> v_N = op(v_N, @tree_reduce(op, a, N, $(cld(Ntot, Nacc)), $Nacc))

        end
        v = @tree_reduce(op, v, 1, $Nacc, 1)
        i = i1 + (n >> $Ptot) * $Ntot - 1
        i == iN && return v
        for i in i+1:iN
            ai = @inbounds A[inds[i]]
            v = op(v, f(ai))
        end
        return v
    end
end

```

The Base.sum from your PR correspond to Nacc=8 (8 accumulators) and Ntot=8 (no intern tree reduction). You mb\_3 function correspond to Ntot=16 and Nacc=4.

Here are my benchmarks:

```julia
n = 2^15
a = rand(Float64, n)
for Ntot in (4, 8, 16, 32)
    for Nacc in (1, 2, 4, 8)
        if Nacc <= Ntot
            println("n=$n, a[1:n], view(a, 1:2:n) (above, below) Ntot=$Ntot, Nacc=$Nacc")
            @btime reduce_abstracted(+, $a, $(Val(Ntot)), $(Val(Nacc)))
            @btime reduce_abstracted(+, $(view(a, 1:2:n)), $(Val(Ntot)), $(Val(Nacc)))
        end
    end
end

```

```repl
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=4, Nacc=1
  3.294 μs (0 allocations: 0 bytes)
  2.475 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=4, Nacc=2
  3.405 μs (0 allocations: 0 bytes)
  2.486 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=4, Nacc=4
  3.649 μs (0 allocations: 0 bytes)
  2.449 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=8, Nacc=1
  4.808 μs (0 allocations: 0 bytes)
  2.130 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=8, Nacc=2
  2.396 μs (0 allocations: 0 bytes)
  2.181 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=8, Nacc=4
  2.147 μs (0 allocations: 0 bytes)
  4.372 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=8, Nacc=8
  3.427 μs (0 allocations: 0 bytes)
  2.225 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=16, Nacc=1
  4.323 μs (0 allocations: 0 bytes)
  2.227 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=16, Nacc=2
  2.061 μs (0 allocations: 0 bytes)
  2.123 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=16, Nacc=4
  1.739 μs (0 allocations: 0 bytes)
  4.345 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=16, Nacc=8
  2.468 μs (0 allocations: 0 bytes)
  4.318 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=32, Nacc=1
  4.236 μs (0 allocations: 0 bytes)
  2.214 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=32, Nacc=2
  2.270 μs (0 allocations: 0 bytes)
  2.369 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=32, Nacc=4
  1.348 μs (0 allocations: 0 bytes)
  3.978 μs (0 allocations: 0 bytes)
n=32768, a[1:n], view(a, 1:2:n) (above, below) Ntot=32, Nacc=8
  1.651 μs (0 allocations: 0 bytes)
  3.959 μs (0 allocations: 0 bytes)

```

For comparison, the Base.sum from your PR on Float64:

```julia-auto
@btime Base.sum($a)
@btime Base.sum($(view(a, 1:2:n)))

```

```repl
3.459 μs (0 allocations: 0 bytes)
16295.798570780385

  4.040 μs (0 allocations: 0 bytes)
8128.148967556634

```

This is different from what I got because of the additionnal pairwise reduction I guess.  
Here is what the current Base.sum doing on Float64:

```julia-auto
2.488 μs (0 allocations: 0 bytes)
16410.273204764926

 6.539 μs (0 allocations: 0 bytes)
8170.225306474282

```

and just for completeness, the results from @simd:

```julia-auto
a=rand(Float32, 2^15)
function simdsum(a::AbstractArray{T}) where T
    s = T(0)
    @simd for x in a
        s += x
    end
    return s
end
@btime simdsum($a)
@btime simdsum($(view(a, 1:2:2^15)))

```

```repl
1.058 μs (0 allocations: 0 bytes)
16345.918f0

6.549 μs (0 allocations: 0 bytes)
8140.3384f0

```

On my Dell machine, the current `sum` implementation underperforms for non-view floating-point arrays. I’m curious whether you observe similar benchmark trends on your system:

- **Val(16), Val(2)**: Good all-around performance across array types
- **Val(32), Val(4)**: Optimal performance specifically for contiguous bitstype arrays

I’ve found that `@simd` performs poorly on views and larger element types, but on my machine it’s **unbeatable for contiguous bitstype arrays** —with one exception: very small arrays benefit more from vectorized loads followed by reduction.

For context, I started a discussion on this topic a few days ago (before discovering this thread):

> [@Imroving Base.mapreduce by 20% on vectors](https://discourse.julialang.org/t/imroving-base-mapreduce-by-20-on-vectors/134478):
>
> Hello everyone, I’ve been exploring reduction operations on GPUs recently, and this led me to investigate CPU performance as well. I discovered that we can improve mapreduce in Base while maintaining full precision for floating-point operations. My analysis suggests there may be an alignment issue with SIMD in the current Base implementation, which leads to suboptimal performance for common computational patterns. Additionally, I’ve implemented specialized functions for small arrays (size \< 32…

I also documented my initial attempt to outperform `mapreduce` on contiguous bitstype arrays using `@simd`:  
[https://epilliat.github.io/coding/notebooks/mapreduce\_cpu\_perf.html](https://epilliat.github.io/coding/notebooks/mapreduce_cpu_perf.html)

A key finding: **alignment is critical** for strong `@simd` performance, at least on my hardware. The current `Base.sum` implementation doesn’t properly initialize accumulators for 1024-element blocks, which significantly degrades performance.

[Previous page](https://discourse.julialang.org/t/performance-challenge-can-you-write-a-faster-sum/130456.md?page=1)