8 March 2026 · 6 min
Suppose that you have a record of your sales per day. You might want to get a running record where, for each day, you are told how many sales you have made since the start of the year.
| day | sales per day | running sales |
|---|---|---|
| 1 | 10$ | 10 $ |
| 2 | 15$ | 25 $ |
| 3 | 5$ | 30 $ |
Such an operation is called a prefix sum or a scan.
Implementing it in C is not difficult. It is a simple loop.
for (size_t i = 1; i < length; i++) {
data[i] += data[i - 1];
}
How fast can this function be? We can derive a speed limit rather simply: to compute the current value, you must have computed the previous one, and so forth.
data[0] -> data[1] -> data[2] -> ...
At best, you require one CPU cycle per entry in your table. Thus, on a 4 GHz processor, you might process 4 billion integer values per second. It is an upper bound but you might be able to reach close to it in practice on many modern systems. Of course, there are other instructions involved such as loads, stores and branching, but our processors can execute many instructions per cycle and they can predict branches effectively. So you should be able to process billions of integers per second on most processors today.
Not bad! But can we do better?
We can use SIMD instructions. SIMD instructions are special instructions that process several values at once. All 64-bit ARM processors support NEON instructions. NEON instructions can process four integers at once, if they are packed in one SIMD register.
But how do you do the prefix sum on a 4-value register? You can do it with two shifts and two additions. In theory, it scales as log(N) where N is the number elements in a vector register.
input = [A B C D]
shift1 = [0 A B C]
sum1 = [A A+B B+C C+D]
shift2 = [0 0 A B+A]
result = [A A+B A+B+C A+B+C+D]
You can then extract the last value (A+B+C+D) and broadcast it to all positions so that you can add it to the next value.
Is this faster than the scalar approach? We have 4 instructions in sequence, plus at least one instruction if you want to use the total sum in the next block of four values.
Thus the SIMD approach might be worse. It is disappointing.
A solution might be the scale up over many more integer values.
Consider ARM NEON which has interleaved load and store instructions. If you can load 16 values at once, and get all of the first values together, all of the second values together, and so forth.
original data : ABCD EFGH IJKL MNOP
loaded data : AEIM BFJN CGKO DHLP
Then I can do a prefix sum over the four blocks in parallel. It takes three instructions. At the end of the three instructions, we have one register which contains the local sums:
A+B+C+D E+F+G+H I+J+K+L M+N+O+P
And then we can apply our prefix sum recipe on this register (4 instructions). You might end up with something like 8 sequential instructions per block of 16 values.
It is theoretically twice as fast as the scalar approach.
In C with instrinsics, you might code it as follows.
void neon_prefixsum_fast(uint32_t *data, size_t length) {
uint32x4_t zero = {0, 0, 0, 0};
uint32x4_t prev = {0, 0, 0, 0};
for (size_t i = 0; i < length / 16; i++) {
uint32x4x4_t vals = vld4q_u32(data + 16 * i);
// Prefix sum inside each transposed ("vertical") lane
vals.val[1] = vaddq_u32(vals.val[1], vals.val[0]);
vals.val[2] = vaddq_u32(vals.val[2], vals.val[1]);
vals.val[3] = vaddq_u32(vals.val[3], vals.val[2]);
// Now vals.val[3] contains the four local prefix sums:
// vals.val[3] = [s0=A+B+C+D, s1=E+F+G+H,
// s2=I+J+K+L, s3=M+N+O+P]
// Compute prefix sum across the four local sums
uint32x4_t off = vextq_u32(zero, vals.val[3], 3);
uint32x4_t ps = vaddq_u32(vals.val[3], off);
off = vextq_u32(zero, ps, 2);
ps = vaddq_u32(ps, off);
// Now ps contains cumulative sums across the four groups
// Add the incoming carry from the previous 16-element block
ps = vaddq_u32(ps, prev);
// Prepare carry for next block: broadcast the last lane of ps
prev = vdupq_laneq_u32(ps, 3);
// The add vector to apply to the original lanes is the
// prefix up to previous group
uint32x4_t add = vextq_u32(prev, ps, 3);
// Apply carry/offset to each of the four transposed lanes
vals.val[0] = vaddq_u32(vals.val[0], add);
vals.val[1] = vaddq_u32(vals.val[1], add);
vals.val[2] = vaddq_u32(vals.val[2], add);
vals.val[3] = vaddq_u32(vals.val[3], add);
// Store back the four lanes (interleaved)
vst4q_u32(data + 16 * i, vals);
}
scalar_prefixsum_leftover(data, length, 16);
}
Let us try it out on an Apple M4 processor (4.5 GHz).
| method | billions of values/s |
|---|---|
| scalar | 3.9 |
| naive SIMD | 3.6 |
| fast SIMD | 8.9 |
So the SIMD approach is about 2.3 times faster than the scalar approach. Not bad.
My source code is available on GitHub.
Appendix. Instrinsics
| Intrinsic | What it does |
|---|---|
vld4q_u32 |
Loads 16 consecutive 32-bit unsigned integers from memory and deinterleaves them into 4 separate uint32x4_t vectors (lane 0 = elements 0,4,8,12,…; lane 1 = 1,5,9,13,… etc.). |
vaddq_u32 |
Adds corresponding 32-bit unsigned integer lanes from two vectors (a[i] + b[i] for each of 4 lanes). |
vextq_u32 |
Extracts (concatenates a and b, then takes 4 lanes starting from lane n of the 8-lane concatenation). Used to implement shifts/rotates by inserting zeros (when a is zero vector). |
vdupq_laneq_u32 |
Broadcasts (duplicates) the value from the specified lane (0–3) of the input vector to all 4 lanes of the result. |
vdupq_n_u32 (implied usage) |
Sets all 4 lanes of the result to the same scalar value (commonly used for zero or broadcast). |
Daniel Lemire, "Prefix sums at tens of gigabytes per second with ARM NEON," in Daniel Lemire's blog, March 8, 2026, https://lemire.me/blog/2026/03/08/prefix-sums-at-tens-of-gigabytes-per-second-with-arm-neon/.
[BibTeX]
I’m compelled to show this program in APL, whose serious implementations also use SIMD instructions, but implicitly:
+\
I know which one I prefer.
I don’t think any APLs have SIMD plus-scans, although I can see why you might think so. My own BQN (which spells it +`) has supported it in CBQN for not quite 3 years now. Dyalog doesn’t, and it’s the only commercial APL that’s any good at SIMD, as others like APL2 and APL*PLUS have not seen development in a long time (free ones like GNU APL and Kap aren’t too performance-focused and don’t do manual SIMD). Nor does J (where it’s +/\). I expect at least Shakti K does, though I don’t have access to it or K4. No support in ngn/k or the fork growler/k; I took part in some discussions on how to do it with SWAR since these implementations don’t use intrinsics, but no one went forward with it. I do see +\ support in Goal, but SIMD there is limited to 16-byte SSE registers and the integers are 8 bytes so that’s pretty underwhelming for scans.
This operation is a little tougher than you might think because APL only exposes smaller integer types as a subset of floats for optimization, so it needs to detect when the sum exceeds the range of the type and widen appropriately. I’ve only known a good way to do the check since last year when Nick Nickolov (the ngn of ngn/k) suggested using the result values: the trick is to compare whether the result minus the argument overflows at any position, so that you implicitly test against the previous result value but no shifting is needed. That made it into the strided scans (e.g. along the first column) I did a little after but I haven’t gone back to our 1-dimensional scans, so those are still using a conservative ahead-of-time check on absolute values. I remember I did some cleaning up of scans when I worked at Dyalog but that was more about removing indirection and using the basic scalar C loop whenever possible. Even min- and max-scans, which don’t have any overflow difficulties, use scalar loops.
I was curious if I could shorten the dependency chain in the lanewise prefix sum from 3 to 2, but it seems to have no effect (on an M1). The compiler appears to restore the original 3-step scan in the assembler output.
// 2 concurrents
vals.val[1] = vaddq_u32(vals.val[1], vals.val[0]);
vals.val[3] = vaddq_u32(vals.val[3], vals.val[2]);
// 2 concurrents
vals.val[2] = vaddq_u32(vals.val[2], vals.val[1]);
vals.val[3] = vaddq_u32(vals.val[3], vals.val[1]);