`mapreduce` uses `@simd` which can cause different results depending on optimizations

[LilithHafner]

This produced incorrect answers in a real-world workflow. This is a simplified MWE:

$ julia --startup=no
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.11.5 (2025-04-14)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> x = Float32[429.17752, 429.2375, 429.29248, 430.21002, 430.995, 431.16498, 431.025, 430.995, 430.995, 430.995, 430.77252, 430.575, 430.255, 430.0, 430.375, 430.535];

julia> f = identity
identity (generic function with 1 method)

julia> sum(x -> identity(x), x)
6886.6f0

julia> sum(x -> f(x), x)
6886.6006f0

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (aarch64-linux-gnu)
  CPU: 64 × Neoverse-N1
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, neoverse-n1)
Threads: 64 default, 0 interactive, 32 GC (on 64 virtual cores)
Environment:
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_CONDAPKG_EXE = /home/ubuntu/temple/.venv/bin/python
  JULIA_PROJECT = @.
  JULIA_NUM_THREADS = auto

[oscardssmith]

as discussed on triage, this is a result of LLVM reassocation that is explicitly made legal by Base.add_sum. As such, math-mode=ieee causes this to return the same result. I think the only thing to potentially do here is to document the functions in Base which may return different results on different architecutres, and people that always want consistent results between systems need to use math-mode=ieee

[LilithHafner]

This is also different results based on different inference results, not just different systems. e.g. you could get one result with a production implementation and then add an @show statement which hurts inlining and inference and consequently get a different result. Then the presence of @show (or other spooky action at a distance that impacts compiler heuristics) would change the returned value.

[mbauman]

The element type thing is different, and Edit: oh both are type instabilities — I think this is something we should definitely fix; and #58418 addresses it. The fact that we can commute + but not other operators is a lot trickier. We could get consistent results if nothing ever commuted, with a ~1.5-2x perf regression.

[Seelengrab]

Are we sure that reassociation is the root cause? What about f being a non-const global, could that prevent inlining (which would lead to reassocation being blocked)?

[KristofferC]

Isn't all of this completely expected by the use of @simd in

julia/base/reduce.jl

Lines 251 to 265 in 4f1e471

	@noinline function mapreduce_impl(f, op, A::AbstractArrayOrBroadcasted,
	ifirst::Integer, ilast::Integer, blksize::Int)
	if ifirst == ilast
	@inbounds a1 = A[ifirst]
	return mapreduce_first(f, op, a1)
	elseif ilast - ifirst < blksize
	# sequential portion
	@inbounds a1 = A[ifirst]
	@inbounds a2 = A[ifirst+1]
	v = op(f(a1), f(a2))
	@simd for i = ifirst + 2 : ilast
	@inbounds ai = A[i]
	v = op(v, f(ai))
	end
	return v

?

[oscardssmith]

yes.

[KristofferC]

So is there anything to do here?

[oscardssmith]

I don't think so. (other than docs)

[mbauman]

I don't think so. (other than docs)

I'm not convinced. I don't think we should be using @simd here. I think we can get close enough to peak performance without it to make this worthwhile.