Introduction

On a CPU, a NumPy backtest runs in float64 by default: a 53-bit significand, integers exact to 2^{53}\approx 9\times10^{15}, relative precision \approx 1.1\times10^{-16}. At price scale — a cryptocurrency close near 30{,}000, a cumulative sum reaching 10^{14} — double precision has margin to spare, and one rarely thinks about it.

Apple's GPU removes that margin entirely, because it offers no float64 at all. The Metal Shading Language has no scalar double type [1], and the array frameworks that target Apple Silicon — MLX [4], and the Metal backends of the major tensor libraries — inherit the restriction. There is no flag to enable double precision and no slow-but-correct fallback: any computation that reaches the GPU happens in binary32, with a 24-bit significand (integers exact only to 2^{24}=16{,}777{,}216) and a relative precision of about 1.19\times10^{-7}, roughly seven significant decimal digits [8].

The trap is that binary32 is almost always fine. Prices, returns, PnL, and Sharpe ratios all live comfortably where seven digits is plenty, so a naive GPU port passes a smoke test and returns sane-looking output. The failure is localized to exactly one operation: the place where a hidden intermediate crosses the 2^{24} ceiling. And that operation is the one a performance-minded quant vectorizes first — the O(n) prefix-sum reduction of a moving average.

This paper is a controlled autopsy of that single operation. We take the linearly weighted moving average (WMA), the building block of the Hull moving average (HMA) [7], and compute the same mathematical quantity three ways on a seeded synthetic price series: (i) a float64 direct windowed sum (the trusted reference); (ii) the O(n) prefix-sum formulation entirely in binary32 (the trap); and (iii) a binary32 direct windowed convolution (the fix). We then quantify, from one seeded run, four things: where the prefix-sum intermediate actually lives relative to the binary32 ceiling; how large the resulting relative error is at the WMA and Hull levels; the falsifiable growth law of the error with series length and price scale; and — because a wrong moving average still looks like a moving average — a discrete-invariant verification, trade-count parity of an HMA-cross strategy, that turns a silent 399\times error into a screaming disagreement and an acceptable 10^{-7} error into a small, bounded, explainable one.

Scope and honesty constraints.

Two constraints bound the claims. First, we do not run on a GPU. binary32 GPU arithmetic is modeled with NumPy [5] float32 on the CPU; the cancellation is a property of IEEE-754 binary32, not of any particular device, and reproduces bit-for-bit in spirit without hardware. Second, this paper measures no wall-clock time and makes no speed or speedup claim: prior engineering work motivates moving such sweeps to the GPU for speed, but speed is out of scope here. The contribution is strictly the numerical-error mechanism and the discrete-invariant verification methodology. All results derive from one seeded run of one public harness (scripts/run_all.py), with the underlying methods in scripts/wma.py and scripts/strategy.py.

binary32, the exact-integer ceiling, and ULP arithmetic

IEEE-754 binary32 stores a number as a sign, an 8-bit exponent, and a 23-bit stored significand; with the implicit leading bit the significand is effectively 24 bits [8]. Two consequences drive everything that follows.

Exact-integer ceiling. Consecutive integers are representable exactly only up to 2^{24}=16{,}777{,}216. Beyond it the representable grid coarsens: between 2^{24} and 2^{25} only even integers exist, then only multiples of four, and so on.

Unit in the last place (ULP). For a value v the gap between adjacent binary32 numbers near v is \begin{equation} \label{eq:ulp} \mathrm{ULP}_{32}(v) = 2^{\lfloor \log_2 |v|\rfloor - 23}, \end{equation} so the absolute rounding error of storing v in binary32 is on the order of \mathrm{ULP}_{32}(v), which scales with |v|. The relative precision is fixed: NumPy reports the binary32 machine epsilon as 2^{-23}\approx1.19\times10^{-7} (the spacing above 1), and the unit roundoff is half of it. When two nearly equal large numbers are subtracted, the absolute errors do not cancel; they persist at the scale of the operands, while the true difference can be arbitrarily smaller. If that difference falls below \mathrm{ULP}_{32} of the operands, it is buried in noise. This is catastrophic cancellation [3, 6], and it is the entire failure.

The classical mitigation for accumulated summation error is compensated (Kahan) summation [9], or accumulation in a wider type. Neither is available here: Metal offers no float64 accumulator, and — as we show — the correct remedy is not more precision but a reformulation that never manufactures the large intermediate in the first place.

The prefix-sum WMA and where it overflows

The WMA and its Hull composition.

The linearly weighted moving average over a window of length p ending at bar t is \begin{equation} \label{eq:wma} \mathrm{WMA}_t = \frac{\sum_{k=1}^{p} k\, x_{t-p+k}}{p(p+1)/2}, \end{equation} a convolution of the price x with the pre-normalized linear kernel w_k = k/[p(p+1)/2], \sum_k w_k = 1. The Hull moving average nests three of these, and its HMA3 variant nests four [7]: \begin{equation} \label{eq:hma} \mathrm{HMA}(x,p) = \mathrm{WMA}\!\Big(2\,\mathrm{WMA}(x,\tfrac p2)-\mathrm{WMA}(x,p),\,\lfloor\sqrt p\,\rceil\Big). \end{equation} Because HMA/HMA3 are compositions of WMAs, each pipeline (fp64, prefix-sum fp32, conv fp32) is structurally identical: only the WMA backend changes, so the error of the backend propagates through the composition.

The tempting O(n) formulation.

Sweeping thousands of parameter combinations over 150{,}000 bars, the WMAs are the cost, so the instinct — correct on a CPU — is to make each WMA O(n) with prefix sums [2]. With global bar index j, precompute two cumulative sums once, \begin{equation} \label{eq:prefix} S_1[t] = \sum_{j\le t} x_j, \qquad S_2[t] = \sum_{j\le t} j\, x_j, \end{equation} and recover any window's weighted sum as differences and index shifts, \begin{equation} \label{eq:recover} \mathrm{wsum}_t = \big(S_2[t{+}1]-S_2[t{-}p{+}1]\big) - (t-p)\big(S_1[t{+}1]-S_1[t{-}p{+}1]\big), \end{equation} then \mathrm{WMA}_t = \mathrm{wsum}_t /[p(p+1)/2]. It vectorizes beautifully, maps onto the GPU's parallel-scan primitive, and in exact (or float64) arithmetic is exactly right. It is also the single worst thing to do in binary32, and the reason hides in S_2.

The hidden 10^{14} intermediate.

The term j\, x_j grows without staying at price scale: with j up to 150{,}000 and x\approx 30{,}000 the last term alone is \approx 4.5\times10^{9}, and S_2 is the running sum of n such terms, reaching \approx x\, n^2/2. Table 1 places the measured magnitudes next to the binary32 ceiling. The reference S_2 reaches 3.259\times10^{14}, which is 7.29 orders of magnitude past 2^{24}=1.678\times10^{7}. One ULP near S_2 is \mathrm{ULP}_{32}(S_2)=3.355\times10^{7} (exactly 2^{25}, since S_2 lies between 2^{48} and 2^{49}), so every stored value of S_2 is rounded to the nearest multiple of \approx 3.4\times10^{7}.

Now compare that noise to the signal being recovered. For a window of p=32 the true weighted sum \mathrm{wsum}_t has magnitude \approx \bar x\, p(p+1)/2 = 1.564\times10^{7} (mean price \bar x = 29{,}628). The ratio \begin{equation} \label{eq:noise-signal} \frac{\mathrm{ULP}_{32}(S_2)}{\text{signal}} = \frac{3.355\times10^{7}}{1.564\times10^{7}} = 2.99 \end{equation} exceeds one: the rounding of each subtracted operand is already larger than the entire windowed quantity being extracted. This is catastrophic cancellation in its purest form — the answer is smaller than the noise floor of the numbers it is subtracted from. (The end-to-end binary32 S_2 differs from its fp64 value by only 3.22\times10^{-5} in relative terms; the disaster is not in computing S_2 but in subtracting two of them to recover something 10^{7} times smaller than their magnitude.)

The fix: a direct windowed convolution

The reflex, on seeing the error, is to reach for more precision — accumulate in float64, or use Kahan summation [9]. On Metal the first is unavailable, and neither is needed, because the problem was never the number of bits. The prefix-sum trick manufactures 10^{14}-scale intermediates and then subtracts them back down; those magnitudes are an artifact of the algorithm, not of the answer. A formulation that never creates them is fine in binary32.

That formulation is the WMA definition itself, summed locally: \begin{equation} \label{eq:conv} \mathrm{WMA}_t = \sum_{k=0}^{p-1} w_k\, x_{t+k}, \qquad w_k = \frac{k+1}{p(p+1)/2}, \end{equation} a direct windowed convolution (the conv1d formulation). Each output is a sum of at most p terms, each of order x/p; every partial sum stays near price scale (\sim\!3\times10^{4}), never within six orders of magnitude of the binary32 ceiling. There is nothing to cancel because nothing ever inflated.

Table 2 reports the maximum relative error against the float64 reference across all 18 WMA periods and all 80 (\text{fast}, \text{slow}) Hull combinations of the sweep (periods 8–220). The prefix-sum WMA is wrong by a maximum relative factor of 399.3 (median 57.9), and the Hull composition amplifies it to a maximum of 3416 (median 734) as the error propagates through three or four nested WMAs and the internal 2\,\mathrm{WMA}-\mathrm{WMA} difference. The convolution, by contrast, holds a maximum relative error of 7.1\times10^{-7} at the WMA level and 7.8\times10^{-7} at the Hull level — both within an order of magnitude of the binary32 epsilon 1.19\times10^{-7}, i.e. at the noise floor. The same math, computed by a different sum, moves the relative error by nine orders of magnitude. Crucially, the prefix-sum output never overflows to infinity or NaN: it is finite and smooth, so it passes every sanity check short of a direct comparison against a trusted reference.

The cliff: error versus length and price scale

The mechanism of Section 3 makes a falsifiable prediction. The cancellation noise is \mathrm{ULP}_{32}(S_2)\sim \varepsilon_{32}\, x\, n^2/2, while the recovered window stays at price scale. Therefore the prefix-sum absolute error should grow like x\, n^2 (linearly in price, quadratically in length), and the relative error should grow like n^2 and be nearly price-invariant (signal and noise both scale with x). The convolution error, tied only to the local sum of p like-scale terms, should stay flat. We test this by sweeping series length at fixed price and price at fixed length, all at p=32 (Table 3).

The length sweep confirms the cliff. As n grows from 10{,}000 to 150{,}000 (15\times), the prefix-sum relative error grows from 0.0854 to 90.46 — a factor of 1060, which is even steeper than the n^2 prediction of 225 (the additional growth comes from the widening dominance of the accumulated cumsum noise over the price-scale signal as n increases). Over the same sweep the convolution relative error moves from 2.33\times10^{-7} to 2.77\times10^{-7}, a factor of 1.19: flat. The price sweep confirms the price-invariance of the relative error: the prefix-sum absolute error grows roughly 1000\times across the price decades (2{,}496 at close 30 to 2.53\times10^{6} at close 30{,}000), while the relative error stays in the band 56–90, and the convolution relative error stays pinned near 2.9\times10^{-7} at every scale. The mechanism is not asserted; it is measured and it obeys its own growth law.

Verification by a discrete invariant: trade-count parity

A moving average wrong by 399\times still looks like a moving average — it is finite, smooth, and roughly in range on the bars where cancellation is mild — so it cannot be audited by eye. Correctness needs an invariant that a downstream, discrete part of the pipeline exposes. A backtest supplies a natural one: its trades.

We run a deliberately simple, fully deterministic HMA-cross strategy: go long when \mathrm{HMA}(\text{close},\text{fast}) is above \mathrm{HMA3}(\text{close},\text{slow}) and flat otherwise, with one-bar-delayed execution at the close. The strategy's purpose is not alpha but verification: crossovers quantize the smooth indicator into discrete trade events, so trade-count parity between an fp32 pipeline and the fp64 reference is an oracle that a smooth-looking equity curve is not. A correct fp32 formulation moves at most a handful of borderline crossovers by one bar; a numerically broken one rewrites the trade list wholesale. We compare, per (\text{fast},\text{slow}) combination and in aggregate over all 80: the maximum absolute trade-count drift; the maximum PnL delta in percentage points; and the shifted-trade fraction, the symmetric difference of trade-entry sets divided by their total.

Table 4 is decisive. The convolution pipeline produces 194{,}945 trades against the reference's 194{,}946 — a total drift of one trade — with a maximum per-combo drift of 2 trades, a maximum PnL delta of 0.246 percentage points, and only 231 of 389{,}891 trade entries shifted (0.059\%). Every one of those is a borderline crossover where h-h3 was itself within the \sim\!10^{-6} relative wobble of zero; small, bounded, and explainable. The prefix-sum pipeline is a different animal entirely: it produces 1{,}011{,}802 trades — more than five times the reference — with a maximum per-combo drift of 12{,}852 trades, a maximum PnL delta of 31.3 percentage points, and 1{,}170{,}582 of 1{,}206{,}748 entries shifted, a 97.0\% rewrite of the trade list. The discrete invariant converts a silent 399\times indicator error into an unmissable failure, and simultaneously certifies the convolution pipeline as correct: it is the oracle the equity curve refused to be. We note that “fraction of combinations that differ” is a decoy metric — with thousands of crossovers per combo, essentially every combo shows some difference under both pipelines — which is why we report by how much the pipelines disagree, not whether they do.

Discussion

Know where your intermediates live, not just your inputs and outputs.

The inputs (prices \sim\!10^{4}) and outputs (a WMA \sim\!10^{4}) were both comfortably inside binary32's exact range. The disaster lived entirely in a hidden intermediate, S_2\sim10^{14}, that neither the API nor the array types made visible. Before trusting any binary32 reduction, ask what the largest partial sum reaches; if it crosses \sim\!10^{7}, the formulation — not the precision — is the problem.

Prefer formulations that keep magnitudes bounded.

Direct convolution over prefix sums; local windows over global scans; centering or differencing before summing rather than after. “Big-then-cancel” is the anti-pattern [6]. The numerically safe algorithm is often the one that looks asymptotically worse on paper (O(np) here versus O(n)) precisely because it never manufactures a quantity it must later cancel away. This is why compensated summation [9], which would help a genuine long-sum accuracy problem, is beside the point: there is no accuracy problem in computing S_2 (its relative drift is 3\times10^{-5}); the problem is the subtraction of two such S_2 values, and no accumulator width short of one that holds the difference exactly repairs a cancellation whose operands are correct.

Validate against an fp64 oracle through a discrete invariant.

Do not compare curves; compare something quantized and downstream — trade count, number of crossovers, position-change events. A discrete invariant turns a silent 399\times error into an assertion failure and an acceptable 10^{-7} error into a small, bounded, explainable delta. The convolution's 231 shifted entries out of 389{,}891 are visibly the fingerprint of correct binary32 arithmetic; the prefix-sum's 97\% rewrite is visibly not. The health metric must be by how much two pipelines disagree, never whether they do at all.

The general regime shift.

The seductive pitch of GPU backtesting is “one big matrix”: stack every parameter combination into a tensor and let the hardware eat it. Moving the same array expression to Metal silently changes the numerical regime underneath the code — from float64 with a 2^{53} integer ceiling to binary32 with a 2^{24} one — and the identical line, \operatorname{cumsum}(j\,x), goes from exact to garbage with nothing in the syntax to warn you. binary32 does not fail loudly; it fails politely, with finite, plausible numbers.

Limitations

• Modeled, not measured on a GPU. We reproduce binary32 arithmetic with NumPy float32 on the CPU. This is faithful to IEEE-754 binary32 semantics but does not capture device-specific quirks (fused multiply-add contraction, non-IEEE fast-math flags, reduction-order nondeterminism) that a real Metal/MLX kernel might add. Those could only make a naive prefix-sum worse, not better; the qualitative conclusion is unchanged.

• No performance claim. We measure no wall-clock time and quote no speedup. The motivation for GPU sweeps is speed, established by prior engineering work, but this study is purely about the numerical mechanism and its verification.

• One synthetic DGP. Numbers are specific to a single seeded geometric-random-walk price series (close \approx 30{,}000, per-bar volatility 5\times10^{-4}, n=150{,}000). The magnitudes are data-dependent: the relative error scales with n^2 and the absolute error with price, so a shorter series or a lower-priced instrument reduces the error, and a longer or higher-priced one worsens it. The mechanism — catastrophic cancellation once S_2 crosses 2^{24} — is invariant; only its magnitude moves with the data.

• One indicator family and one strategy. We study the linear-WMA / Hull-MA family and a single HMA-cross strategy as the discrete oracle. The prefix-sum cancellation and the direct-convolution fix generalize to any windowed linear filter, but the specific error and trade-parity numbers are tied to these choices.

Conclusion

Apple's Metal GPU has no float64, so a vectorized backtest ported to it runs in IEEE-754 binary32, where integers are exact only to 2^{24}=16{,}777{,}216. The most tempting way to vectorize a weighted moving average — the O(n) prefix-sum reduction — builds a hidden intermediate S_2=\operatorname{cumsum}(j\,x_j) that reaches 3.26\times10^{14}, seven orders past that ceiling, and recovers each window by subtracting two such \sim\!10^{14} numbers whose one-ULP rounding noise (3.36\times10^{7}) is 2.14\times the answer. The result is catastrophic cancellation: a maximum relative error of 399 at the WMA level and 3416 at the Hull level, delivered as finite, smooth, plausible output. The cure is not more bits, of which Metal has none, but a different sum — a direct windowed convolution that keeps every partial sum at price scale and drives the relative error to 7\times10^{-7}, the binary32 noise floor. The mechanism is falsifiable and obeys its growth law: the prefix-sum relative error explodes superquadratically with series length (1060\times as n grows 15\times) while the convolution stays flat (1.19\times). And because a wrong moving average cannot be audited by eye, correctness is established through a discrete invariant: HMA-cross trade-count parity shows the convolution disagreeing with the fp64 reference on 231 of 389{,}891 trade entries (0.059\%), all borderline crossovers, while the prefix-sum pipeline rewrites 97.0\% of the trade list. If a backtest gets dramatically faster and its numbers still look fine, that is not confirmation: on Metal, “looks fine” is exactly what a relative error of hundreds looks like.

Reproducibility.

All numbers derive from one seeded run. scripts/run_all.py regenerates results/results.json from seed 0 (Python 3.12.3, NumPy 2.5.1) deterministically; the WMA/HMA backends and the error metrics are in scripts/wma.py, and the strategy and trade-parity oracle in scripts/strategy.py. The synthetic series is a geometric random walk from numpy.random.default_rng(0) (initial price 30{,}000, per-bar log-return standard deviation 5\times10^{-4}, zero drift, n=150{,}000). tests/ contains deterministic invariant tests for the three arithmetic pipelines and the parity report.

References