% The GPU Precision Trap: catastrophic cancellation in an fp32 prefix-sum WMA.
% Controlled, seeded numerical-analysis study. Every number in this manuscript
% comes from results/results.json (seed 0), produced by scripts/run_all.py.
% Primary arXiv category: cs.NA (Numerical Analysis).
% Cross-list: q-fin.CP (Computational Finance).
\documentclass[11pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amssymb}
\usepackage{booktabs}
\usepackage[numbers,round]{natbib}
\usepackage{microtype}
\usepackage[hidelinks]{hyperref}
\setlength{\emergencystretch}{2em}

\newcommand{\fp}[1]{\texttt{#1}}
\newcommand{\ulp}{\mathrm{ULP}}
\newcommand{\fltz}{\mathrm{fl}_{32}}

\title{The GPU Precision Trap:\\
Catastrophic Cancellation in an fp32 Prefix-Sum Moving Average,\\
and a Discrete-Invariant Test That Catches It}
\author{Eugen Soloviov}
\date{July 2026}

\begin{document}
\maketitle

\begin{abstract}
Apple's Metal GPU has no \fp{float64}: every array a backtest touches on Apple
Silicon is IEEE-754 binary32. The most tempting way to vectorize a linearly
weighted moving average (WMA) is the $O(n)$ prefix-sum trick, which precomputes
$S_2 = \operatorname{cumsum}(j\,x_j)$ and recovers each window by subtraction.
We show, on a seeded synthetic price series ($n=150{,}000$ bars, close
$\approx 30{,}000$), that this formulation is numerically catastrophic in
binary32 while looking perfectly healthy. The hidden intermediate $S_2$ climbs
to $3.259\times10^{14}$, which is $7.29$ orders of magnitude past binary32's
exact-integer ceiling $2^{24}=16{,}777{,}216$; one unit in the last place near
$S_2$ is $3.355\times10^{7}$, so recovering a window (true magnitude
$\approx 1.564\times10^{7}$) subtracts two $\sim\!10^{14}$ numbers whose
rounding noise is $2.14\times$ the answer --- textbook catastrophic
cancellation. The maximum relative error against a \fp{float64} reference
reaches $399$ at the raw-WMA level and $3416$ at the composed Hull-MA level, yet
the output is finite, smooth, and plausible. The cure is not more bits (there
are none) but a different sum: a direct windowed convolution keeps every partial
sum at price scale and drives the relative error to $7.1\times10^{-7}$, the
binary32 noise floor. We confirm the falsifiable mechanism through a length
sweep: the prefix-sum relative error grows superquadratically ($1059\times$ as
$n$ grows $15\times$, versus $n^2$'s $225\times$), while the convolution error
is essentially flat ($1.19\times$). Because a smooth-looking curve cannot be
audited by eye, we verify correctness through a \emph{discrete invariant} ---
the trade count of an HMA-cross strategy. Against the \fp{float64} reference,
the convolution disagrees on only $231$ of $389{,}891$ trade entries
($0.059\%$), a bounded handful of borderline crossovers, whereas the prefix-sum
pipeline rewrites $97.0\%$ of the trade list and produces $1{,}011{,}802$ trades
against the reference's $194{,}946$. We model binary32 GPU arithmetic with NumPy
\fp{float32} on the CPU: the failure is a property of IEEE-754 binary32, not of
any device, and we measure no GPU wall-clock time and make no speed claim. This
study accompanies a \href{https://marketmaker.cc}{marketmaker.cc} blog post.
\end{abstract}

\section{Introduction}
\label{sec:intro}

On a CPU, a NumPy backtest runs in \fp{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 \fp{float64} at
all. The Metal Shading Language has no scalar \fp{double} type
\citep{metal-msl}, and the array frameworks that target Apple Silicon ---
MLX~\citep{mlx2023}, 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 \citep{ieee754-2019}.

The trap is that binary32 is \emph{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)~\citep{hull2005active}, and compute the \emph{same} mathematical
quantity three ways on a seeded synthetic price series: (i) a \fp{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 \emph{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.

\paragraph{Scope and honesty constraints.} Two constraints bound the claims.
First, we do \emph{not} run on a GPU. binary32 GPU arithmetic is modeled with
NumPy~\citep{harris2020numpy} \fp{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 \emph{no} wall-clock time and
makes \emph{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 (\texttt{scripts/run\_all.py}), with the underlying methods in
\texttt{scripts/wma.py} and \texttt{scripts/strategy.py}.

\section{binary32, the exact-integer ceiling, and ULP arithmetic}
\label{sec:fp32}

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 \citep{ieee754-2019}. Two consequences drive everything
that follows.

\emph{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.

\emph{Unit in the last place (ULP).} For a value $v$ the gap between adjacent
binary32 numbers near $v$ is
\begin{equation}
\label{eq:ulp}
\ulp_{32}(v) = 2^{\lfloor \log_2 |v|\rfloor - 23},
\end{equation}
so the \emph{absolute} rounding error of storing $v$ in binary32 is on the order
of $\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 \emph{not} cancel;
they persist at the scale of the operands, while the true difference can be
arbitrarily smaller. If that difference falls below $\ulp_{32}$ of the operands,
it is buried in noise. This is \emph{catastrophic cancellation}
\citep{goldberg1991floating,higham2002accuracy}, and it is the entire failure.

The classical mitigation for accumulated summation error is compensated (Kahan)
summation \citep{kahan1965summation}, or accumulation in a wider type. Neither is
available here: Metal offers no \fp{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.

\section{The prefix-sum WMA and where it overflows}
\label{sec:trap}

\paragraph{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 \citep{hull2005active}:
\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.

\paragraph{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
\citep{blelloch1990prefix}. 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 \fp{float64})
arithmetic is exactly right. It is also the single worst thing to do in
binary32, and the reason hides in $S_2$.

\paragraph{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~\ref{tab:magnitudes} 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 $\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{\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 \emph{subtracting} two of them to recover something $10^{7}$ times
smaller than their magnitude.)

\begin{table}[t]
\centering
\caption{The hidden intermediate versus the binary32 ceiling ($n=150{,}000$,
close $\approx 30{,}000$, seed 0). The prefix-sum WMA builds $S_2$ to
$3.26\times10^{14}$, over seven orders past the last exactly representable
integer; one ULP there ($3.36\times10^{7}$) is $2.14\times$ the windowed
quantity (at $p=32$) it must recover by subtraction.}
\label{tab:magnitudes}
\begin{tabular}{lr}
\toprule
Quantity & Value \\
\midrule
binary32 exact-integer ceiling $2^{24}$        & $16{,}777{,}216$ \\
binary32 machine epsilon $2^{-23}$             & $1.192\times10^{-7}$ \\
mean price $\bar x$                             & $29{,}628.13$ \\
$S_2 = \operatorname{cumsum}(j\,x_j)$ (fp64)    & $3.2595\times10^{14}$ \\
$S_2$ in fp32, relative drift vs fp64          & $3.22\times10^{-5}$ \\
orders of magnitude of $S_2$ past $2^{24}$      & $7.29$ \\
$\ulp_{32}(S_2)$                                & $3.3554\times10^{7}$ \\
recovered window signal ($p=32$)               & $1.5644\times10^{7}$ \\
noise-to-signal $\ulp_{32}(S_2)/\text{signal}$ & $2.14$ \\
\bottomrule
\end{tabular}
\end{table}

\section{The fix: a direct windowed convolution}
\label{sec:fix}

The reflex, on seeing the error, is to reach for more precision --- accumulate
in \fp{float64}, or use Kahan summation \citep{kahan1965summation}. On Metal the
first is unavailable, and neither is needed, because the problem was never the
number of bits. The prefix-sum trick \emph{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 \fp{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~\ref{tab:trapfix} reports the maximum relative error against the
\fp{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
\fp{NaN}: it is finite and smooth, so it passes every sanity check short of a
direct comparison against a trusted reference.

\begin{table}[t]
\centering
\caption{The trap versus the fix: maximum, median, and minimum relative error
against the \fp{float64} reference, at the raw-WMA level (over $18$ periods) and
at the composed Hull level (over $80$ fast/slow combinations). The prefix-sum
formulation is wrong by factors of hundreds to thousands; the direct convolution
sits at the binary32 noise floor ($\sim\!10^{-7}$).}
\label{tab:trapfix}
\begin{tabular}{llrrr}
\toprule
Level & Pipeline & Max rel.\ err.\ & Median & Min \\
\midrule
WMA  & prefix-sum fp32 & $399.3$              & $57.92$              & $8.809$ \\
WMA  & conv fp32       & $7.141\times10^{-7}$ & $3.406\times10^{-7}$ & $1.691\times10^{-7}$ \\
\midrule
Hull & prefix-sum fp32 & $3416$               & $733.8$             & $429.5$ \\
Hull & conv fp32       & $7.781\times10^{-7}$ & $5.766\times10^{-7}$ & $4.370\times10^{-7}$ \\
\bottomrule
\end{tabular}
\end{table}

\section{The cliff: error versus length and price scale}
\label{sec:cliff}

The mechanism of Section~\ref{sec:trap} makes a falsifiable prediction. The
cancellation noise is $\ulp_{32}(S_2)\sim \varepsilon_{32}\, x\, n^2/2$, while
the recovered window stays at price scale. Therefore the prefix-sum
\emph{absolute} error should grow like $x\, n^2$ (linearly in price, quadratically
in length), and the \emph{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~\ref{tab:cliff}).

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 \emph{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 \emph{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.

\begin{table}[t]
\centering
\caption{The cliff at $p=32$ (seed 0). \emph{Length sweep} (top, close
$\approx 30{,}000$): the prefix-sum relative error grows superquadratically with
$n$ while the convolution error is flat. \emph{Price sweep} (bottom,
$n=150{,}000$): the prefix-sum \emph{absolute} error grows with price but the
\emph{relative} error is price-invariant; the convolution stays at the noise
floor throughout.}
\label{tab:cliff}
\begin{tabular}{rrrrr}
\toprule
$n$ / close & $S_2$ final & prefix rel.\ & prefix abs.\ & conv rel.\ \\
\midrule
\multicolumn{5}{l}{\emph{Length sweep} (close $\approx 30{,}000$)}\\
$10{,}000$  & $1.498\times10^{12}$ & $0.0854$ & $2{,}648$              & $2.330\times10^{-7}$ \\
$30{,}000$  & $1.411\times10^{13}$ & $0.7265$ & $23{,}340$             & $2.330\times10^{-7}$ \\
$75{,}000$  & $8.664\times10^{13}$ & $16.65$  & $5.000\times10^{5}$    & $2.647\times10^{-7}$ \\
$150{,}000$ & $3.259\times10^{14}$ & $90.46$  & $2.530\times10^{6}$    & $2.770\times10^{-7}$ \\
\midrule
\multicolumn{5}{l}{\emph{Price sweep} ($n=150{,}000$)}\\
close $30$     & $3.259\times10^{11}$ & $89.85$ & $2{,}496$           & $2.919\times10^{-7}$ \\
close $300$    & $3.259\times10^{12}$ & $68.97$ & $18{,}924$          & $2.913\times10^{-7}$ \\
close $3{,}000$  & $3.259\times10^{13}$ & $56.14$ & $1.555\times10^{5}$ & $2.874\times10^{-7}$ \\
close $30{,}000$ & $3.259\times10^{14}$ & $90.46$ & $2.530\times10^{6}$ & $2.770\times10^{-7}$ \\
\bottomrule
\end{tabular}
\end{table}

\section{Verification by a discrete invariant: trade-count parity}
\label{sec:oracle}

A moving average wrong by $399\times$ still \emph{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, \emph{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 \emph{correct} fp32 formulation
moves at most a handful of borderline crossovers by one bar; a numerically
\emph{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 \emph{shifted-trade fraction}, the symmetric difference of trade-entry
sets divided by their total.

Table~\ref{tab:oracle} 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 \emph{by how
much} the pipelines disagree, not \emph{whether} they do.

\begin{table}[t]
\centering
\caption{Discrete-invariant verification: HMA-cross trade parity against the
\fp{float64} reference over $80$ fast/slow combinations (seed 0). The fp64
reference produces $194{,}946$ trades. The convolution pipeline is
indistinguishable up to a bounded handful of borderline crossovers; the
prefix-sum pipeline rewrites the trade list wholesale.}
\label{tab:oracle}
\begin{tabular}{lrr}
\toprule
Metric & conv fp32 vs.\ fp64 & prefix-sum fp32 vs.\ fp64 \\
\midrule
Total trades (ref.\ $=194{,}946$)   & $194{,}945$ & $1{,}011{,}802$ \\
Max abs.\ trade-count drift (combo) & $2$         & $12{,}852$ \\
Max abs.\ PnL delta (p.p.)          & $0.246$     & $31.31$ \\
Shifted trade entries               & $231$       & $1{,}170{,}582$ \\
Total entries (both pipelines)      & $389{,}891$ & $1{,}206{,}748$ \\
Shifted-trade fraction              & $0.059\%$   & $97.0\%$ \\
\bottomrule
\end{tabular}
\end{table}

\section{Discussion}
\label{sec:discussion}

\paragraph{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.

\paragraph{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
\citep{higham2002accuracy}. 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 \citep{kahan1965summation}, 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.

\paragraph{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 \emph{by how much} two pipelines disagree, never \emph{whether}
they do at all.

\paragraph{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 \fp{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.

\section{Limitations}
\label{sec:limitations}

\begin{itemize}
\item \textbf{Modeled, not measured on a GPU.} We reproduce binary32 arithmetic
  with NumPy \fp{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.
\item \textbf{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.
\item \textbf{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 \emph{mechanism} --- catastrophic
  cancellation once $S_2$ crosses $2^{24}$ --- is invariant; only its magnitude
  moves with the data.
\item \textbf{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.
\end{itemize}

\section{Conclusion}
\label{sec:conclusion}

Apple's Metal GPU has no \fp{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.

\paragraph{Reproducibility.} All numbers derive from one seeded run.
\texttt{scripts/run\_all.py} regenerates \texttt{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 \texttt{scripts/wma.py}, and the strategy and trade-parity
oracle in \texttt{scripts/strategy.py}. The synthetic series is a geometric
random walk from \texttt{numpy.random.default\_rng(0)} (initial price $30{,}000$,
per-bar log-return standard deviation $5\times10^{-4}$, zero drift,
$n=150{,}000$). \texttt{tests/} contains deterministic invariant tests for the
three arithmetic pipelines and the parity report.

\bibliographystyle{plainnat}
\bibliography{refs}

\end{document}
