SIMD Math in @simd Loops

AppleAccelerate.SIMDMath exposes scalar math functions that LLVM replaces with Apple's SIMD math routines when they appear inside a loop that vectorises.

This is the non-array counterpart to Array Operations. Where AppleAccelerate.log!(out, X) needs a whole array in and a whole array out, SIMDMath.log(x) is an ordinary scalar call you can put in the middle of a loop.

When to use it

SIMDMath is for the cases where the array API does not apply:

• strided or otherwise non-contiguous access;
• values computed on the fly rather than read out of an array;
• loops with control flow, or over iterables that are not dense Arrays;
• hot paths where you cannot allocate and cannot preallocate a temporary.

In those situations the realistic alternative is a scalar Base loop, and SIMDMath is roughly 2–4x faster for Float32 and 1.2–2x for Float64.

using AppleAccelerate
using AppleAccelerate.SIMDMath: log

# weights applied on the fly: vForce would need a temporary for log.(X) first
function weighted_logsum(X, W)
    u = zero(eltype(X))
    @simd for i in eachindex(X, W)
        @inbounds u += W[i] * log(X[i])
    end
    u
end

X = collect(1.0:1000.0)
W = fill(0.5, 1000)
weighted_logsum(X, W) ≈ sum(W .* Base.log.(X))

true

Scoping it to one loop: @simdmath

using AppleAccelerate.SIMDMath: log replaces log for the whole enclosing module. Since these functions only have Float32/Float64 methods, every other log(2) or rem(i, n) in that module then throws a MethodError. SIMDMath.@simdmath confines the substitution to one expression instead, and imports nothing:

using AppleAccelerate
using AppleAccelerate.SIMDMath: @simdmath

function weighted_logsum(X, W)
    u = zero(eltype(X))
    @simdmath @simd for i in eachindex(X, W)
        @inbounds u += W[i] * log(X[i])^W[i]     # _simd_log_d2, _simd_pow_d2
    end
    u
end

X = rand(1000) .+ 1; W = rand(1000)
weighted_logsum(X, W) ≈ sum(W .* log.(X) .^ W)

true

It can go outside or inside @simd. It does not add @simd or @inbounds for you, and the loop still has to vectorise for any of this to matter.

A rewritten call uses the SIMD routine only when its arguments are all Float32 or all Float64; for anything else it falls back to the Base function. So index arithmetic like rem(i, 4), integer powers x^n, and complex or BigFloat values keep working inside the block. (nextafter and remainder have no Base counterpart to fall back to.)

The rewrite is purely syntactic: a local variable or argument that is named log and then called is rewritten too. The accuracy caveats below apply unchanged.

Accuracy

Confirming you got the vectorised form

A call only becomes a SIMD call if the surrounding loop actually vectorises, which in practice means @simd and usually @inbounds. Nothing warns you if it does not — the code stays correct and simply runs at scalar speed.

To check, look for the symbol in the generated code:

using InteractiveUtils
asm = sprint(code_native, weighted_logsum, (Vector{Float64}, Vector{Float64}))  # what `@code_native` prints
occursin("_simd_log_d2", asm)   # the scalar `log` became a 2-lane SIMD call

true

Float32 runs 4 lanes at a time (_simd_*_f4), Float64 runs 2 (_simd_*_d2).

Available functions

All are defined for Float32 and Float64 only.

nextafter, pow and remainder have no Base counterpart with matching semantics, so they are named after their C equivalents. Use pow(x, y) rather than x^y, or let @simdmath rewrite x^y for you.

That is every one- and two-argument routine <simd/math.h> offers that is worth calling. copysign, min/max (C's fmin/fmax) and fdim have no _simd_* entry point at all — the header implements them inline, and LLVM already vectorises the Base versions natively.

Deliberately absent:

• sqrt, floor, ceil, round, trunc, fma, abs — LLVM already lowers these to native instructions, so a library call would be slower.
• erf, erfc, tgamma — <simd/math.h> declares _simd_erf_f4 and friends, but they measure at 0.97–0.99x of a plain scalar libm loop: they are scalar loops behind a vector signature and buy nothing.
• lgamma — writes the global signgam, so it cannot be declared side-effect-free.

Measured speedups

out[i] = f(x[i]) over N = 100,000 on an M-series machine. SM/Base is SIMDMath against a scalar Base loop (higher is better); SM/vForce is against AppleAccelerate.f! (below 1.0 means vForce wins).

Reading this:

• SM/vForce is below 1.0 in every single row. vForce wins across the board, which is why the guidance above is to prefer it whenever it applies. nextafter is the extreme case at 0.18x — always prefer AppleAccelerate.nextafter!.
• The Float32 trigonometric functions are where this shines against Base, especially sinpi/cospi/tanpi, where Base is unusually slow.
• cosh and hypot barely beat Base in Float32 because Base implements them natively there — those loops vectorise with no library call at all, so there is little left to win.
• hypot and exp10 have no vForce equivalent, so SIMDMath is the only accelerated option for them.

How it works

LLVM's loop vectoriser can replace a scalar call inside a vectorising loop with a call to an equivalent SIMD routine, if it is told which routine to use. clang spells this -fveclib=Accelerate. The underlying mechanism is the LLVM vector function ABI (VFABI), which is driven entirely by IR metadata, so Julia can reach it through llvmcall — no compiler flag and no patched Julia required.

Each function is emitted as a scalar call to the real libm symbol (logf), carrying a "vector-function-abi-variant" attribute naming the SIMD replacement:

_ZGV_LLVM_N4v_logf(_simd_log_f4)

which reads as "_LLVM_ ISA, Not masked, 4 lanes, one vector operand, mapping scalar logf onto vector _simd_log_f4". The vector routine must also be declared in the module and anchored in @llvm.compiler.used, or it gets stripped before the vectoriser runs and the mapping is silently dropped.

The routines come from libsystem_m's <simd/math.h> rather than Accelerate's own vfp.h. Both are public Apple API, but vfp.h is Float32-only and covers fewer functions, whereas <simd/math.h> covers both widths — and since Float64 is Julia's default, using one source for both keeps accuracy consistent across types.

AppleAccelerate.SIMDMath

julia> using AppleAccelerate.SIMDMath: log

julia> function weighted_logsum(X, W)  # no temporary for log.(X), unlike the array API
           u = zero(eltype(X))
           @simd for i in eachindex(X, W)
               @inbounds u += W[i] * log(X[i])
           end
           u
       end;

julia> X = collect(1.0:1000.0); W = fill(0.5, 1000);

julia> weighted_logsum(X, W) ≈ sum(W .* Base.log.(X))
true

using InteractiveUtils
@code_native weighted_logsum(X, W)   # expect a call to _simd_log_d2

SIMDMath.@simdmath expr

julia> using AppleAccelerate.SIMDMath: @simdmath

julia> function weighted_logsum(X, W)
           u = zero(eltype(X))
           @simdmath @simd for i in eachindex(X, W)
               @inbounds u += W[i] * log(X[i])^W[i]
           end
           u
       end;

julia> X = [1.5, 2.5, 3.5, 4.5]; W = [0.5, 1.0, 1.5, 2.0];

julia> weighted_logsum(X, W) ≈ sum(W .* log.(X) .^ W)   # `log` and `^` are SIMDMath's here
true

julia> @simdmath rem(7, 4)                              # not Float32/Float64: falls back to Base
3