Fast Stencil Computations on a Single Arbitrarily Moving Interval
arXiv:2609.14879v1 [cs.DS] 14 Sep 2026
Aaron Gregory∗ Stony Brook University
Abstract A stencil computation repeatedly updates every cell of a grid from its neighbours’ values at the previous timestep. Simulating T steps on N cells directly costs Θ(N T ), and a line of work beginning with Ahmad et al. [4] reduces this by composing many timesteps into one linear operator and applying it with a Fast Fourier Transform. That technique needs to know which cells will still obey the same operator when the composed step ends. In a free-boundary problem they do not: the region governed by a given rule is determined by the solution and moves as it evolves. We study the case of one spatial dimension, a three-point stencil with time-varying coefficients, and a computed region that is a single interval whose two endpoints move by arbitrary amounts at every step, revealed online. Let B be the horizon plus the total variation of the boundary trajectory. We give a schedule whose total work is O((B + N ) log T log(N + B)) and whose span is O(T log T log(N + B)), and we prove that the values it computes are exact. The best existing bound for a region that moves requires its boundary to travel at most one cell per timestep. We drop that requirement and lose nothing by it: a boundary obeying the requirement has B ≤ 3T , so our bound stays near-linear on every trajectory the earlier result covers. On the trajectories it does not cover, B grows only by the distance the boundary actually travels — one jump of width N costs T + 2N . The reason total variation suffices is that everything the two endpoints touch over a time window of any length lies in two intervals, one per endpoint. We also √ show this cannot be relaxed: with p regions the bound degrades by a factor p, and at p = T there is an instance on which the work is Θ(T 3/2 ) while B + N = Θ(T ). All results are machine-checked in Lean 4, apart from the classical convolution bound of Theorem 3.5, which is imported as an interface.
1
Introduction
A stencil computation evaluates a recurrence on a spatial grid over a time horizon: the value at a cell at time t is a fixed function of the values at that cell and its neighbours at time t − 1. Write N for the number of cells and T for the number of timesteps. Evaluating the recurrence as written — every cell at every step — costs Θ(N T ), and for many scientific and financial kernels that is the dominant cost. When the update is linear and does not vary across space, k consecutive steps compose into a single linear operator whose support has radius k, and applying that operator to the whole grid is one convolution, which a Fast Fourier Transform evaluates in O(N log N ) work rather than Θ(N k). We call one such composed application a superstep. Trading k timesteps for one superstep is the idea behind the FFT-based stencil algorithms of [4, 6, 5]. ∗
1
The difficulty. A superstep can only be taken over a set of cells that obey the same operator for all k of the steps it spans. In a free-boundary problem this set is not known in advance: the region governed by a given rule is determined by the solution itself and moves as the computation proceeds. An algorithm that wants long supersteps must therefore decide, before evaluating them, which cells will still be governed by the same operator when the superstep ends — and must do so without inspecting the whole grid at every step, since that would already cost Θ(N T ). The setting. We work in one spatial dimension. The computed region at time t is a single half-open interval I(t) = [at , bt ) ⊂ Z, nonempty. Its two endpoints are otherwise unconstrained: each may move by an arbitrary amount in either direction at each step, and the trajectory is revealed online, one step at a time. The stencil support at each step is contained in {−1, 0, 1}, with coefficients that may vary with time. Cells whose dependencies reach outside the previous region are supplied by a boundary oracle (Section 3). The algorithm must produce every value on I(t) for every t ≤ T . The parameter.
Write ∆at = |at − at−1 | and ∆bt = |bt − bt−1 |, and set B = T +
T X
∆at + ∆bt ,
(1)
t=1
the horizon plus the total variation of the boundary trajectory. The T term is forced: with coefficients that change every step, the schedule must be read, so Ω(T ) work is unavoidable even for a stationary region.
1.1
Prior work and what it assumes
The FFT-based line optimizes the cost of a homogeneous superstep. [4] applies it to periodic linear stencils and [6] generates implementations; [5] handles supplied boundary conditions by a top-down recursive decomposition of the spacetime domain. Each of these fixes the spatial partition in advance, which is what a moving region violates. Formulations in which the update is nonlinear, so that the set obeying a given rule is determined by the solution, are treated in [2]. [7] allows the operator to vary in time and to differ between regions, with coefficients homogeneous inside each region, and allows the regions themselves to vary in time; it introduces the binary time-product tree we use in Section 4 to obtain composed operators. We take the binary-forking cost model and its work-efficient FFT from [3]. Of these, only [7] admits a region that moves at all, and it is the result we compare against. Its Theorem 4.8 states that when T ≤ N/B0 , AperiodicSolve computes aT in O(N log N log T + BΣ ) work and O(T log N log log N ) span,
(2)
where Bt counts the P boundary cells of every region together with the cells governed by the boundary condition, BΣ = t≤T Bt , and B0 is that count at t = 0. It carries a regularity condition on how a region may move: the boundary at time t + 1 must lie inside the region of influence of the boundary at time t, so an endpoint travels at most one stencil radius per step. Specialized to the present setting — one dimension, one interval, a three-point stencil — the regularity condition says each endpoint moves by at most one cell per step, BΣ = Θ(T ), and (2) reads O(N log N log T ) work for T = O(N ). We remove the regularity condition, and the restriction on T with it. An endpoint may move any distance at any step, and the horizon is unbounded. What replaces the condition 2
Direct simulation Bentley et al. [7] This paper
Work
Span
Θ(N T ) O(N log N log T ) O((B+N ) log T log(N +B))
Θ(T ) O(T log N log log N ) O(T log T log(N +B))
Boundary may move freely one cell per step freely
Table 1: Bounds for one dimension, one interval and a three-point stencil; row two is Theorem 4.8 of [7] specialized to that setting, and holds only for T = O(N ). Where row two applies, row three matches it to within a logarithmic factor, because a boundary moving one cell per step has B ≤ 3T ; where it does not, row three still applies and charges the distance the boundary travels (Proposition 3.8). is the parameter: a trajectory that obeys it has B ≤ 3T , so (2) and our bound describe the same regime there, while a trajectory that breaks it even once leaves (2) with nothing to say and is charged by us at B (Proposition 3.8).
1.2
Results
Our results are the following. (i) A closed form for the cursor. The set of cells at time f that can be computed from time f − w without consulting anything outside the region is an interval, and both of its endpoints are an explicit maximum and minimum over the intervening slices (Definition 4.1). For arbitrary regions this set is an intersection of erosions, and a radius bound for it must be obtained by induction; here it is read off a formula. (ii) Exactness (Theorem 5.2), which uses neither linearity of the update, nor any bound on the boundary’s speed, nor any relation between the stencil and the boundary’s motion. (iii) A write count of 12(L + 1)B, with no hypothesis on the trajectory (Theorem 6.6), where L = ⌊log2 T ⌋. (iv) Work and span (Theorems 6.8 and 6.9), as in Table 1. (v) A necessity √ result (Proposition 7.1): √ with p regions the write count degrades by a factor p, and at p = T the bound fails by Θ( T ). (vi) An identity (Theorem 3.7) showing that the total variation is within a factor two of the total outward movement plus the initial width, which is what justifies charging contraction as well as expansion. (vii) The cost of the generality (Proposition 3.8): under the regularity condition B ≤ 3T , and a single unconstrained step is enough to leave that condition behind. Scope. Three limitations, stated once. The stencil is three-point and the dimension is one; Section 8 explains why the argument does not survive d ≥ 2. The span is linear in T , because the region at time t is not determined until time t − 1 has been computed, so the timesteps serialize; when the trajectory is instead supplied in advance and the update is linear, Corollary 6.10 gives polylogarithmic span. And what we prove about correctness is Theorem 5.2 together with the covering identity of Proposition 4.3: every cell is supplied by the oracle, filled by exactly one shell, or deferred, and each shell’s values are determined by the slice it is solved from. Algorithms 1 and 2 are presented but not formalized: we prove neither that the recursive and per-time views coincide, nor an end-to-end statement that running the recursion produces every requested value.
3
left strip
right strip
window
I(t), at a time t with ν(t) = 3
S0
a
a
b
S1
S2 C23 (t) S2
S1
S0
b
(a) whatever the endpoints do, everything they touch lies in two strips
(b) the shells are differences of nested cursors, hence two arms each
Figure 1: (a) The computed region I(t) = [at , bt ) over a window of nine steps, including one large jump of the left endpoint. Every cell the endpoints sweep lies in one of two strips, whose widths are the endpoints’ total travel over the window; this holds for a window of any length, which is Lemma 6.1. (b) At a trigger time t, the nested cursors C1 (t) ⊃ C2 (t) ⊃ C4 (t) ⊃ C8 (t) cut the region into shells, each a pair of arms, plus a deep cursor handed to the next level. Roadmap. Section 2 places the result among neighbouring lines of work. Section 3 fixes the model, the oracle, and the cost of a superstep. Section 4 gives the cursor, the schedule, and a worked example; a reader who wants only the algorithm can stop there. Section 5 proves exactness. Section 6 is the accounting, and runs in one chain: a shell is charged to a backward window, the activity in a window is two intervals, and the windows at a given level tile. Section 7 shows the single-interval hypothesis cannot be dropped, and Section 8 discusses span and dimension.
2
Related work
Section 1.1 covered the FFT-based line this paper extends and the bound it compares against. Three neighbouring lines are worth separating from it. Direct algorithms. The standard stencil algorithms evaluate the recurrence as written and perform Θ(N T ) work: nested loops, tiled loops — time skewing [14] and diamond tiling [9] among them — and cache-oblivious recursive decomposition of the spacetime domain [11], realized in the Pochoir compiler [13]. These differ in the order they visit the N T spacetime cells, which is what determines their locality and parallelism, but not in how many they visit. Their generality is the reason: they place no condition on the dimension, the support, the linearity of the update, or the shape of the region, and so they remain the fallback whenever the hypotheses of Section 3 fail. Everything asymptotically faster buys the speedup with a structural assumption, and the assumption is the interesting part. Tracking a moving interface. A separate literature is devoted to representing a boundary that moves. Level-set methods [12] carry the interface as the zero set of an auxiliary field, and narrowband variants [1] restrict the update to a neighbourhood of that set so the per-step cost scales with the interface rather than the grid. The concern there is the geometry: how to represent a front that merges, splits or develops curvature, and how to advance it stably. That is not the question here.
4
We take the boundary as given — an oracle reports it, one step at a time — and ask how cheaply the interior can be advanced once it is known. The two are complementary, and neither subsumes the other: a narrow band still advances the interior one timestep at a time, which is exactly the Θ(T ) factor a superstep removes. Cost model. Work and span are measured in the binary-forking model [8, 3], in which spawning n threads costs Θ(log n) span. We use it only through Theorem 3.5; nothing in Section 6 depends on the choice.
3
Model and cost
3.1
Regions and trajectories
Definition 3.1 (Region and trajectory). A region is a pair of integers a < b, written I = [a, b), with cells {x ∈ Z : a ≤ x < b} and width wid(I) = b − a. A trajectory assigns a region I(t) = [at , bt ) to each t ∈ {0, . . . , T }. We write N = wid(I(0)). Definition 3.2 (Sweep, changed cells, frontier). The sweep of an endpoint from u to v is sw(u, v) = [min(u, v), max(u, v)), which has |u − v| cells. The changed cells at step t are Dt = sw(at−1 , at ) ∪ sw(bt−1 , bt ), and the frontier of a region is Γ(I) = {a − 1, a, b − 1, b}. Sweeps are half-open so that consecutive sweeps of the same endpoint meet exactly: [u, v) ∪ [v, w) = [u, w). With inclusive endpoints they would leave a one-cell gap at each v, and Lemma 6.1 would be false. I(t) = [at , bt ) N = wid(I(0)) Dt Γ(I) B ν(t), L
the region at time t initial width cells that changed at t frontier of a region boundary cost, (1) 2-adic valuation, ⌊log2 T ⌋
Cw (t) (w) (w) at , bt Si (t) A[s, w] HR [s, w]
cursor: solvable from t − w its endpoints shell filled at level i activity over a window its hull, dilated by R
Objects indexed by a time are written as functions of it, with the level or lookback as a subscript; (0) the cursor’s endpoints decorate the region’s own endpoints, and at = at .
3.2
The update, the oracle, and the output
Definition 3.3 (Stencil). A stencil over a value space V is a family of maps φt : (Z → V) → Z → V, one per timestep, each local : if g and h agree on {y : |y − x| ≤ 1} then φt (g)(x) = φt (h)(x). We (k) write φs for the k-fold iterate from time s. Locality is all that correctness uses. Linearity of φt is needed only for the cost model, since it is what allows k steps to be composed into one convolution. Definition 3.4 (Boundary oracle and required output). A cell of I(t) whose radius-one neighbourhood is not contained in I(t − 1) cannot be obtained from the previous slice; its value is supplied by a boundary oracle at unit cost per cell. The algorithm must produce the value at every cell of I(T ).
5
Only the final slice is required. Values at intermediate times are computed where they are needed and nowhere else: a superstep that spans many timesteps produces its output without ever materializing the slices it crosses, and that is what makes a bound below Θ(N T ) possible at all. Requiring every intermediate value would force Θ(N T ) output on its own. The oracle’s workload is part of the cost and is charged in Lemma 4.4: it is at most 2B calls over the horizon.
3.3
Cost model and the superstep primitive
We use the binary-forking model [3] and count arithmetic operations in V at unit cost; the bounds below are operation counts, not bit counts. “Exact” throughout means exact with respect to the recurrence and the schedule: every value produced is the value the recurrence defines, given that each ring operation is performed exactly. Theorem 3.5 (Convolution by FFT [10, 3]). In the binary-forking model, the convolution of two sequences of total length n over a commutative ring admitting a principal 2⌈log2 n⌉ -th root of unity can be computed with O(n log n) work and O(log n) span. Consequently, when the update is linear, a superstep of k timesteps can be evaluated on a block of v contiguous cells with O((v + k + 1) log(2 + v + k)) work and O(log(2 + v + k)) span. Composing k linear updates gives one operator of support radius at most k, so applying it to v contiguous cells is a single convolution of a sequence of length O(v + k) against a kernel of length O(k); Theorem 3.5 is the only property of the update we use beyond locality. The kernel term is carried separately because it does not follow from the block term: a shell can be small while the level that fills it is deep.
3.4
Why total variation is the right charge
One might object that only expansion can force work, since a cell that leaves the region need never be revisited, and that charging contraction inflates B. It does, but only by a factor of two. P Definition 3.6 (Outward movement). out(T ) = Tt=1 (at−1 − at )+ + (bt − bt−1 )+ is the total outward movement of the two endpoints. Theorem 3.7 (Contraction is paid for by expansion). For every trajectory and every T , T X
∆at + ∆bt
= 2 out(T ) + wid(I(0)) − wid(I(T )),
t=1
and hence T + out(T ) ≤ B ≤ 2 T + out(T ) + N . Proof. Apply |x| = 2x+ − x to each endpoint’s increments and telescope; the two telescoping sums contribute wid(I(T )) − wid(I(0)). For the upper bound use wid(I(T )) ≥ 1, which holds since aT < bT ; for the lower bound use x+ ≤ |x|. Proposition 3.8 (What the regularity condition costs). Let σ ≥ 1. (a) If every step moves each endpoint by at most σ — the regularity condition of [7] for a stencil of radius σ — then B ≤ (1 + 2σ)T . For the three-point stencil, σ = 1 and B ≤ 3T . (b) For every N ≥ 1, every T ≥ 1 and every s < T , the trajectory that holds I(t) = [0, N ) for t ≤ s and I(t) = [N, 2N ) thereafter has B = T + 2N , its single largest endpoint displacement being N . 6
Proof. (a) Each step contributes ∆at + ∆bt ≤ 2σ to the sum in (1), and there are T of them. (b) Only the step from s to s + 1 moves either endpoint, and it moves each by exactly N . Part (a) is why the generality is free where the earlier result applies: on any trajectory satisfying the regularity condition, B = Θ(T ), so our bound is O((N +T ) log T log(N +T )) there — within a logarithmic factor of (2), with no restriction on T . Part (b) is the other side: the trajectory sits still, jumps once by its own width, and sits still again, which for N > σ violates the regularity condition at exactly one step and so falls outside (2) entirely. Our bound still applies, and charges T + 2N . Remark 3.9 (B is not a universal lower bound). Theorem 3.7 says B is not wasteful relative to expansion; it does not say it is necessary. Let the region widen by M and narrow again at every step. Then every slice has at most N + M cells, yet B = (M + 1) T. So no bound in terms of the slice size and the horizon controls B, and B is not a lower bound on the work any algorithm must do. It is the right parameter when every cell entering the region carries independent data, and not otherwise.
4
The algorithm
4.1
The cursor in closed form
Definition 4.1 (Cursor). For an end time f and a lookback w, Cw (f ) = (w)
af
(w) (w) af , bf , (w)
= max af −j + j ,
bf
0≤j≤w
= min bf −j − j . 0≤j≤w
We only ever use w ≤ f . Each j contributes the constraint imposed by slice f − j: to have its j-step dependency cone inside that slice, a cell must sit at least j inside it. So Cw (f ) is exactly the set of cells at time f computable from slice f − w using only values interior to the region throughout. Lemma 5.1 gives one direction; for the other, a cell outside the cursor falls short of the binding constraint at some slice f − j, and then the cell j away from it on that side lies outside I(f − j). For example, with w = 3 and left endpoints af −3 , . . . , af = 0, 5, 4, 4, the four constraints are (3) 0 + 3, 5 + 2, 4 + 1, 4 + 0, so af = 7: the slice three steps back is not the binding one. Lemma 4.2 (Nesting). If w ≤ w′ then Cw′ (f ) ⊆ Cw (f ); and Cw (f ) ⊆ I(f ). (w)
(w)
Proof. af is a maximum and bf term of each.
4.2
a minimum over a larger index set; the second claim is the j = 0
The trigger schedule
Let ν(t) be the 2-adic valuation of t, so ν(12) = 2, and let L = ⌊log2 T ⌋. At time t the algorithm fires levels 0, . . . , ν(t) − 1; level i fills the shell Si (t) = C2i (t) \ C2i+1 (t) 7
Algorithm 1 Solve(s, f, X): produce the time-f values on X, given time s Require: [s, f ) an aligned dyadic interval; the time-s values on the radius-(f −s) neighbourhood of X within I(s) are available 1: if f = s + 1 then 2: for x ∈ X do 3: x ← φs (slice s)(x) if x ∈ C1 (f ), else x ← oracle(f, x) 4: end for 5: return 6: end if 7: if X ⊆ Cf −s (f ) then 8: evaluate one superstep from s to f onto X ▷ Theorem 3.5; no oracle value is needed, by Lemma 5.1 9: return 10: end if 11: m ← (s + f )/2 12: Y ← X ⊕ [−(f −m), f −m] ∩ I(m) ▷ what time m must supply 13: Solve(s, m, Y ); Solve(m, f, X) Algorithm 2 The same schedule, viewed one timestep at a time 1: query the oracle for the endpoints at , bt , and for the value at every cell of I(t) \ C1 (t) 2: for i = 0, . . . , ν(t) − 1 do 3: evaluate one superstep from time t − 2i onto Si (t) = C2i (t) \ C2i+1 (t) 4: end for 5: leave C2ν(t) (t) to a longer superstep, which lands at a later time divisible by a higher power of two and skips time t for those cells by one superstep from time t − 2i . By Lemma 4.2 a shell is the difference of nested intervals, hence a union of two intervals — a left arm and a right arm (Figure 1b). Each arm is evaluated by its own superstep, so a level-i call is two calls to the primitive of Theorem 3.5; this doubles the call count, which the bounds below absorb. The schedule fires level i exactly when 2i+1 divides t, that is, once every 2i+1 steps. So a cell is re-derived at level i only O T /2i+1 times over the horizon, while each level-i superstep spans 2i timesteps. Summing over the L + 1 levels is what converts timesteps into logarithmic factors. The top-level call is Solve(0, T, I(T )). The guard on line 8 is the only place a decision is made, and Lemma 5.1 is what justifies it: if X ⊆ Cf −s (f ) then nothing influencing X leaves the region during [s, f ), so the ordinary stencil applies throughout and no boundary value is consulted. Otherwise the interval is halved and the cells are recovered at finer granularity. Flattening the recursion gives the per-time view. At an intermediate time t the cells materialized are those that no longer aligned superstep landing at t can reach, and since an aligned superstep of length 2i landing at t requires 2i | t, the longest available is 2ν(t) . That is the shell decomposition: Proposition 4.3 (Shell decomposition). For every t and every n, the shells S0 (t), . . . , Sn−1 (t) are pairwise disjoint and [ C1 (t) = Si (t) ∪ C2n (t). i<n
Proof. Induction on n, using C2k (t) = Sk (t) ∪ C2k+1 (t), which is Lemma 4.2. For disjointness, a cell of Sj (t) lies in C2j (t) ⊆ C2i+1 (t) for i < j, while cells of Si (t) lie outside C2i+1 (t). 8
So the region splits into three families at each step: the cells outside C1 (t), which the oracle supplies; the shells, which the chain fills, each cell in exactly one; and the deep cursor, solvable from far enough back to be handed to the enclosing level. Nothing is missed and nothing is filled twice. Lemma 4.4 (The oracle’s workload). I(t) \ C1 (t) ≤ 2 + ∆at + ∆bt for every t ≥ 1, so the oracle is called at most 2B times over the horizon. Proof. C1 (t) = [max(at , at−1 +1), min(bt , bt−1 −1)) by Definition 4.1, and it lies inside I(t) by + + Lemma 4.2, so the difference has wid(I(t)) P − |C1 (t)| ≤ (at−1 +1 − at ) + (bt − bt−1 +1) cells. Summing over t and using (1) gives 2T + t (∆at + ∆bt ) ≤ 2B.
4.3
Composed operators
A superstep over an aligned dyadic time interval needs the composed operator and the composed support for that interval. These are maintained in a binary time-product tree [7]: a leaf for [t, t + 1) stores that step’s operator and support, and an internal node for [u, w) with children [u, v), [v, w) stores the composition of the two operators and the Minkowski sum of the two supports. Every call in Algorithm 2 uses an aligned dyadic interval, so each requested interval is already a node of the tree and no recomposition is needed at query time. Because the support here is three-point, composed supports are intervals and their Minkowski sums are computed by adding endpoints, in constant time per node. With supports of unbounded width one must instead convolve zero-one support masks and test positivity, to avoid mistaking coefficient cancellation for the absence of a dependency.
4.4
A worked example
Take T = 8 and I(0) = [4, 8), so N = 4. Let the region be stationary for three steps, then jump: a4 = 0, after which it contracts by one on each of two steps and is stationary thereafter; let bt be fixed at 8 throughout. t
at
bt
∆at + ∆bt
ν(t)
levels fired
0 1 2 3 4 5 6 7 8
4 4 4 4 0 1 2 2 2
8 8 8 8 8 8 8 8 8
— 0 0 0 4 1 1 0 0
— 0 1 0 2 0 1 0 3
— none 0 none 0, 1 none 0 none 0, 1, 2
P Here B = 8 + 6 = 14, against t wid(I(t)) = 49 cells for direct simulation. The jump at t = 4 moves the left endpoint four cells at once, so this trajectory is outside the regularity condition and outside (2); it contributes 4 to B and nothing else. At t = 4 the oracle supplies the four cells the jump exposed plus the two at the ends; at t = 8 the chain fires three levels, whose shells partition C1 (8) and leave C8 (8) to the enclosing level.
9
5
Correctness
Lemma 5.1 (Cone containment). If x ∈ Cw (f ), j ≤ w, and |y − x| ≤ j, then y ∈ I(f − j). (w)
Proof. af
is a maximum including the term af −j + j, so x ≥ af −j + j and y ≥ x − j ≥ af −j . (w)
Dually x < bf
≤ bf −j − j gives y ≤ x + j < bf −j . (k)
(k)
Theorem 5.2 (Exactness). Let g, h : Z → V agree on I(s). Then φs (g)(x) = φs (h)(x) for every x ∈ Ck (s + k). (k)
Proof. By locality and induction on k, φs (g)(x) depends on g only through its restriction to {y : |y − x| ≤ k}. By Lemma 5.1 with j = k, that ball lies in I(s), where g and h agree. A cell solvable from time s is therefore determined by the region’s values at s alone: no boundary datum can change it, and the ordinary stencil applies throughout. The proof uses no induction along dependency chains, no bound on the boundary’s speed, and no linearity.
6
The accounting
We prove Theorem 6.6: the schedule performs O(B log T ) writes, with no hypothesis on the trajectory. The argument is a chain of four steps, and it is worth having in view before the notation arrives. Each shell Si (t) is written at trigger time t, and we charge it to the backward window of states [t − 2i+1 , t]. Lemma 6.3 shows the shell lies in the two-interval hull of the window’s endpoint excursions, dilated by 2i+1 . Lemma 6.1 is what makes that hull the right object: everything the endpoints touch over a window — of any length — lies in those same two intervals, each no wider than the window’s total variation. Fattening two intervals by a radius R adds 4R cells; fattening 2p of them adds 4pR, and that factor p is the whole of Section 7. Finally Lemma 6.4 shows the windows of a given level are disjoint, so each level costs one copy of the total variation, and there are L + 1 levels.
6.1
Activity over a window is two intervals
Fix a window of states s, s + 1, . . . , s + w and let A[s, w] be all activity over it, namely Γ(I(s)) together with Du+1 ∪ Γ(I(u + 1)) for s ≤ u < s + w. Write a, a for the least and greatest of as , . . . , as+w , and b, b likewise. For R ≥ 0 put HR [s, w] = a − 1 − R, a + R ∪ b − 1 − R, b + R . Lemma 6.1 (Chaining). A[s, w] ⊆ H0 [s, w], and every cell within distance R of A[s, w] lies in HR [s, w]. Proof. Induction on w. For every u in the window the frontier cells au − 1, au lie in the first interval and bu − 1, bu in the second, by definition of the extremes. The changed cells Du+1 split into a left sweep between au and au+1 and a right sweep between bu and bu+1 ; each has both endpoints within the corresponding extremes, and being half-open it meets the previous step’s sweep rather than starting a new component. A point within R of a member of an interval lies in that interval’s R-dilation.
10
The window is indexed by its states rather than by its activity times because activity at time u refers to I(u − 1), so a window of w activity times spans w + 1 states. P Lemma 6.2 (The hull is small). HR [s, w] ≤ 4 + 2 s+w u=s+1 (∆au + ∆bu ) + 4R. Proof. Two intervals, each of size (extreme spread) + 2R + 2. Each endpoint’s spread over a window is at most the window’s total variation, by induction on the window length: adjoining a state moves the running maximum and minimum apart by at most the new increment. Two intervals, each charged the whole variation, give the factor 2. The two summands are the two terms of the complexity: the variation is charged to B, and the 4R to the schedule. The coefficient of R is 4 — twice the number of intervals — and, crucially, it does not grow with w.
6.2
The shells are local
Lemma 6.3 (Shell locality). If w ≤ w′ ≤ f then Cw (f ) \ Cw′ (f ) ⊆ Hw′ [f − w′ , w′ ]. In particular Si (t) ⊆ H2·2i [t − 2i+1 , 2i+1 ]. Proof. A cell of the difference fails the longer-lookback cursor on one side. If on the left, then (w′ ) (w) x < af = af −j + j for some j ≤ w′ , so x ≤ a + w′ ; and x ≥ af ≥ af ≥ a. The right side is dual. Without a closed form for the cursor, a bound of this shape has to be obtained by induction along dependency chains, and the constant it yields is larger. Here the radius is the window length itself, read off Definition 4.1.
6.3
Charging
Lemma 6.4 (Level windows tile). Fix ℓ ≥ 1. As t ranges over the multiples of 2ℓ in [1, T ], the backward windows [t − 2ℓ , t] have pairwise disjoint increment sets, so X X (∆au+1 + ∆bu+1 ) ≤ B − T. t
u∈[t−2ℓ , t)
Proof. Distinct multiples of 2ℓ differ by at least 2ℓ , so the half-open increment ranges [t − 2ℓ , t) are disjoint and contained in [0, T ). Charging by trigger time rather than by written slice is what keeps the window backwardlooking. A chronological algorithm on an online trajectory cannot justify a write by activity that has not happened yet, so a window extending forward from the written slice would be unavailable; charging by the pair (trigger time, level) makes the direction explicit, and the windows then tile instead of merely overlapping boundedly. PT ν(t) ≤ (L + 1) T . Lemma 6.5 (Total band radius). t=1 2 P P Proof. 2ν(t) ≤ j≤L [ 2j | t ] 2j . Exchanging the order of summation gives j≤L 2j ⌊T /2j ⌋ ≤ (L + 1)T . Theorem 6.6 (Maintenance write count). The cells materialized at intermediate times satisfy T X X Si (t) ≤ 12 (L + 1) B, with no hypothesis on the trajectory. t=1 i<ν(t)
11
Proof. By Lemma 6.3 each shell lies in a hull, whose size Lemma 6.2 bounds by X 4 + 2 (∆au + ∆bu ) + 8 · 2i . u∈(t−2i+1 , t]
P Sum over i < ν(t) and t ≤ T . The constant terms give 4 t ν(t) ≤ 4(L + 1)T . The variation terms give at most 2(L + 1)(B − T ): by Lemma 6.4 each level contributes one copy, and there are L + 1 levels. The radius terms give at most 8(L + 1)T by Lemma 6.5. The two terms proportional to T contribute (4 + 8)(L + 1)T and the variation term contributes 2(L + 1)(B − T ); majorizing both coefficients by 12 and using T + (B − T ) = B gives 12(L + 1)B. Theorem 6.6 counts the cells materialized at intermediate times. Two families sit outside it and are bounded separately: the oracle’s cells, at most 2B by Lemma 4.4, and the output slice I(T ), at most N + B by Lemma 6.7. What it costs to produce these values is the superstep primitive, charged next.
6.4
Work and span
Lemma 6.7 (Block size). Every block the schedule evaluates has at most N + B cells. P Proof. A shell lies in C2i (t) ⊆ I(t) by Lemma 4.2, and wid(I(t)) ≤ wid(I(0)) + u≤t (∆au + ∆bu ) ≤ N + B by telescoping. Theorem 6.8 (Work). The schedule’s total work is O B + N (L + 1) log(N +B) . Proof. By Theorem 3.5 a level-i call on a block of v cells costs O (v + 2i + 1) log(2 + v + 2i ) . Lemma 6.7 bounds v by N + B, and 2i < 2ν(t) ≤ t ≤ T ≤ B, so every logarithm is O(log(N + B)). Summing over the intermediate calls, the v terms give the write count of Theorem 6.6; the kernel P ν(t) P P i ≤ 2 ≤ (L + 1)T by Lemma 6.5; and the +1 terms give the number 2 terms give t P t i<ν(t) of calls, t ν(t) ≤ (L + 1)T . Both are at most (L + 1)B since T ≤ B. The output slice adds one superstep on at most N + B cells, and the oracle contributes O(B) by Lemma 4.4. Theorem 6.9 (Span). The schedule’s total span is O(log(N +B) · (L + 1) · T ). Proof. Each call contributes span O(log(N +B)) by P Theorem 3.5, using the same bounds on v and 2i as above, and time t makes ν(t) calls; sum using t ν(t) ≤ (L + 1)T . The span is linear in T . That is not an artefact of the analysis: the region at time t is not determined until time t − 1 has been computed, so the timesteps serialize. We do not know a schedule that does better, and we do not prove a matching lower bound. What the bound does not contain is B — an arbitrarily large jump costs work, not depth. When the trajectory is known in advance the serialization disappears: Corollary 6.10 (Depth of a balanced schedule). A schedule organized as a balanced binary tree over 2h supersteps, each of depth δ, has depth (h + 1) δ. Proof. Every root-to-leaf path costs the same, so the depth of a node is the maximum of its children’s depths plus δ; induction on the tree gives (height + 1)δ, and the balanced tree over 2h leaves has height h.
12
When the trajectory is supplied in advance and the update is linear, operator composition is associative over time, so the supersteps may be arranged as such a tree rather than as a chronological loop. With h = ⌈log2 T ⌉ and δ = O(log(N +B)) this gives span O(log T log(N +B)) at the same work. With L = O(log T ) the two theorems read work O((B + N ) log T log(N + B)) ,
span O(T log T log(N + B)) ,
as in Table 1.
7
A single interval is necessary
Lemma 6.1 gives two intervals for one region and 2p for p of them, so the radius term of Lemma 6.2 becomes 4pR and the same argument would give O(p (L + 1)B) for Theorem 6.6. That degradation is attained. Proposition 7.1 (The radius is paid once per component). Let R ≥ 0. (a) If p centers are pairwise at least 2R + 1 apart, their radius-R dilation has exactly p (2R + 1) cells. (b) If every center lies in [u, v], the dilation has at most v − u + 2R + 1 cells, whatever p is. (c) With p = 2n such components of radius R = 2n , fired at each of 2n trigger times, the dilation total is at least 2 · 8n ; for T = 4n this is 2 T 3/2 . Proof. (a) The balls around separated centers are pairwise disjoint and each has 2R + 1 cells. (b) The dilation is contained in [u − R, v + R]. (c) Multiply: 2n trigger times, 2n components, 2 · 2n + 1 cells each, and 2n · 2n · (2 · 2n + 1) = 2 · 8n + 4n . Parts (a) and (b) are the whole of it: a single interval pays the radius once, p separated components pay it p times, and Lemma 6.2’s coefficient 4R becomes 4pR. Part (c) turns that factor into a T 3/2 total, against a boundary cost that such a family keeps at Θ(T ) — the √ singleinterval bound of Theorem 6.6 would give√O(T log T ) for the same data, so the gap is Θ( T ). No polylogarithmic factor absorbs Θ( T ). Note that such a family can keep its mean activity at O(1), so bounding average activity does not rescue the bound. What fails is that nothing ties the number of components to the 2-adic valuation, and a single interval forecloses that by having boundedly many components in every window rather than boundedly many cells.
8
Discussion
Higher dimensions. Lemma 6.1 replaces a union of swept bands by two intervals, which is available only in one dimension. The corresponding statement for d ≥ 2 would convert a (d − 1)dimensional surface into a d-dimensional volume, and the component count that makes Lemma 6.2 work has no analogue there. We do not know whether the obstruction is avoidable. Span. Corollary 6.10 settles the offline case. In the online case we believe the right bound has the shape Õ(m), where m is the number of timesteps at which the boundary actually moves, since the dependency chain through the boundary decisions has depth m; obtaining it appears to require speculating over blocks and verifying, which is sound only under a locality hypothesis on the criterion that moves the boundary. We prove neither that upper bound nor a matching lower bound. 13
When the bound is strong. The work is O((B + N ) log T log(N +B)), so it is near-linear exactly when the boundary’s total travel is comparable to the horizon. A boundary that moves O(1) cells per step on average has B = O(T ) and total work O((N + T ) log T log(N +T )) — whatever its largest single step, since one arbitrarily long jump is absorbed into the sum. The bound is weak when the total travel is superlinear in T , and Remark 3.9 shows that in that regime it can also be far from optimal.
References [1] D. Adalsteinsson and J. A. Sethian. A fast level set method for propagating interfaces. Journal of Computational Physics, 118(2):269–277, 1995. [2] Z. Ahmad, R. Browne, R. Chowdhury, R. Das, Y. Huang, and Y. Zhu. Fast American option pricing using nonlinear stencils. In Proceedings of the 29th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming, pages 316–332, 2024. [3] Z. Ahmad, R. Chowdhury, R. Das, P. Ganapathi, A. Gregory, and M. M. Javanmard. Low-span parallel algorithms for the Binary-Forking model. In Proceedings of the 33rd ACM Symposium on Parallelism in Algorithms and Architectures, pages 22–34, 2021. [4] Z. Ahmad, R. Chowdhury, R. Das, P. Ganapathi, A. Gregory, and Y. Zhu. Fast stencil computations using Fast Fourier Transforms. In Proceedings of the 33rd ACM Symposium on Parallelism in Algorithms and Architectures, pages 8–21, 2021. [5] Z. Ahmad, R. Chowdhury, R. Das, P. Ganapathi, A. Gregory, and Y. Zhu. A fast algorithm for aperiodic linear stencil computation using Fast Fourier Transforms. ACM Transactions on Parallel Computing, 10(4):1–34, 2023. [6] Z. Ahmad, M. M. Javanmard, G. Croisdale, A. Gregory, P. Ganapathi, L.-N. Pouchet, and R. Chowdhury. Fourst: A code generator for FFT-based fast stencil computations. In 2022 IEEE International Symposium on Performance Analysis of Systems and Software (ISPASS), pages 99–108. IEEE, 2022. [7] R. Bentley, R. Chowdhury, A. Gregory, and M. Santomauro. Applying Fast Fourier Transforms to accelerate spatially and temporally inhomogeneous stencil computations. In Proceedings of the 37th ACM Symposium on Parallelism in Algorithms and Architectures (SPAA), pages 17–33, 2025. [8] G. E. Blelloch, J. T. Fineman, Y. Gu, and Y. Sun. Optimal parallel algorithms in the binary-forking model. arXiv preprint arXiv:1903.04650, 2019. [9] U. Bondhugula, V. Bandishti, and I. Pananilath. Diamond tiling: Tiling techniques to maximize parallelism for stencil computations. IEEE Transactions on Parallel and Distributed Systems, 28(5):1285– 1298, 2017. [10] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90):297–301, 1965. [11] M. Frigo and V. Strumpen. Cache oblivious stencil computations. In Proceedings of the 19th International Conference on Supercomputing, pages 361–366, 2005. [12] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988. [13] Y. Tang, R. A. Chowdhury, B. C. Kuszmaul, C.-K. Luk, and C. E. Leiserson. The Pochoir stencil compiler. In Proceedings of the 23rd ACM Symposium on Parallelism in Algorithms and Architectures, pages 117–128, 2011. [14] D. Wonnacott. Achieving scalable locality with time skewing. International Journal of Parallel Programming, 30(3):181–221, 2002.
14