Primitive Recursion without Composition Dynamical Characterizations, from Neural Networks to Polynomial ODEs
Olivier Bournez # Institut Polytechnique de Paris, Ecole Polytechnique, LIX, Palaiseau, France
arXiv:2604.24356v1 [cs.CC] 27 Apr 2026
Abstract What computational mechanisms do recurrent neural networks, polynomial ordinary differential equations, and discrete polynomial maps each bring to the table, and what do they lack? All three are models of computation over the continuum: they operate on real-valued states and evolve by real-valued dynamics, even when the functions we ask them to compute are ultimately discrete. We investigate how these models compare, their strengths, their limitations, and the precise resources on which each one relies, through the lens of primitive recursive functions. We prove that the classical notion of primitive recursion admits equivalent characterizations in all three dynamical frameworks: bounded iteration of a fixed recurrent ReLU network, robust computation by a fixed polynomial ordinary differential equation, and iteration of a fixed polynomial map in discrete time with an externally supplied step-size parameter. In each case, the time bound is itself primitive recursive, composition is not postulated as a closure rule but emerges from the dynamics, and the input is given as a raw integer vector with no auxiliary encoding. At the proof level, every primitive recursive function is first compiled into bounded iteration of a single threshold-affine normal form map, which is then interpreted as a recurrent ReLU computation on the one hand, and as a robust polynomial ODE on the other. The equivalences expose a structural asymmetry between discrete and continuous polynomial computation. We prove that no fixed polynomial map can round uniformly toward the nearest integer, and that none can realize exact phase selection: two operations that polynomial ODEs perform robustly through their continuous-time flow. Each formalism compensates for a limitation that the others do not share: the ReLU gate provides exact branching, continuous time provides autonomous rounding and control, and the step-size parameter recovers both at the cost of discretization precision. Our equivalence theorem characterizes what each resource contributes, and opens the way to dynamical characterizations of subrecursive hierarchies and complexity classes by restricting the time bounds, polynomial degrees, or discretization resources within the same framework. More broadly, the constructions reveal that these real-valued models do not compute by composing subroutines in the classical sense: they compute by shaping the trajectory of a dynamical system, through clocks, phase selectors, stabilization mechanisms, and error correction built into the dynamics itself. This is a mode of computation that differs structurally from symbolic programming, and our equivalence theorem provides a precise framework in which the difference can be studied. 2012 ACM Subject Classification Theory of computation → Models of computation; Theory of computation → Computability; Theory of computation → Complexity classes; Mathematics of computing → Ordinary differential equations; Computer systems organisation → Analog computers; Machine Learning → Neural networks Keywords and phrases Discrete ordinary differential equations, Finite Differences, Implicit complexity, Recursion scheme, Ordinary differential equations, Models of computation, Analog Computations, Formal neural networks Digital Object Identifier 10.4230/LIPIcs...
1
Introduction
Over the integers, the theory of computation rests on a remarkably stable foundation: Turing machines, register machines, the λ -calculus, and many other formalisms are all known to define the same notion of effective computation. Over the reals, the situation is less settled. Several models have been proposed with different motivations, including computable analysis [40] and algebraic models [7], but they are provably inequivalent [40, 17]. Within this landscape, two families of models have © Olivier Bournez; licensed under Creative Commons License CC-BY 4.0 Leibniz International Proceedings in Informatics Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
Olivier Bournez
XX:1
recently emerged as genuinely computational rather than merely descriptive, and both have developed substantial theory. The first comes from neural computation. Already in the early 1990s, recurrent networks with threshold-like or analog activations were shown to possess significant computational power and nontrivial complexity-theoretic structure [37, 4, 30, 36]. More recent architectures, most notably transformer models, rely on substantially different mechanisms, and only make these foundational questions more pressing. What matters here is not just that such models are expressive, but that they are programmable, in a way quite unlike classical symbolic programming. The second family comes from the analog-computing (General Purpose Analog Computer [35]) tradition and, more broadly, from polynomial Ordinary Differential Equations (ODEs). A long line of work has established polynomial ODEs as a robust and machine-independent framework for computability and complexity over the reals [32, 16, 8, 15, 13, 11, 6]. Again, the point is not merely that such systems can simulate classical models, but that computations can be synthesized through mechanisms natural to continuous dynamics: oscillation, stabilization, transfer of intermediate quantities, and correction of errors along the flow. Despite their apparent differences, these two families share a common feature: they realize computation by shaping the evolution of a dynamical system. This is a very different constructive logic from the classical one based on explicit symbolic composition. In recurrent networks and in polynomial ODEs, one often obtains a computation not by writing down its decomposition into subroutines, but by designing a state space, a feedback mechanism, a thresholding behavior, or a flow with the right dynamical properties. A first indication that this viewpoint is genuinely different is that it changes what is easy and what is hard. Rounding to the nearest integer, for instance, is almost free in continuous time (a polynomial ODE built on y′′ = −y contracts to the lattice) but is provably costly or impossible for a pure polynomial map in discrete time. Section 2 makes this asymmetry precise; it is the dynamical fact that organizes the rest of the paper. This leads to a natural foundational question. Instead of taking machines or syntactic closure principles as the primary point of departure, can one find a simple classical class of functions that already sits at the intersection of these dynamical modes of computation? Through this paper, we argue that the right place to begin is to get back to the historical roots of computable theory: the old but remarkably robust class of primitive recursive functions. We assume the reader has some familiarity with computability theory, particularly with primitive recursive functions: [33] provides a thorough reference. Primitive recursion as a dynamical system The class PR of primitive recursive functions is usually defined as the smallest class of functions over the natural numbers containing the zero, successor, and projection functions, and closed under composition and primitive recursion. Many subrecursive and complexity-theoretic classes are obtained by restricting these schemes, including the Grzegorczyk hierarchy [27, 26], polynomial-time functions [22], and polynomial-space functions [39]; see [21] for a broad overview. In all these definitions, closure under composition is taken for granted: this is somehow the one rule that no formalism considers optional. Primitive recursion already admits a dynamical reading. The scheme f (0, x) = g(x),
f (n + 1, x) = h(n, f (n, x), x)
can be viewed as the bounded iteration of the transition (n, a, x) 7→ (n + 1, h(n, a, x), x), where a plays the role of an accumulator: in this view, the computation is not assembled from composed subroutines, but unfolds as a single trajectory. This is reminiscent of how polynomial ODEs compute in the analog tradition, where a single continuous flow replaces any composition rule, and the required
XX:2
Primitive Recursion without Composition
control emerges from the dynamics itself through basins of attraction, stabilization, and autonomous switching mechanisms. The same pattern appears in recurrent neural networks, though this is not always made explicit. A recurrent network with activation ρ repeatedly applies a fixed update rule to a state vector, where one step is a finite composition of affine maps and coordinatewise applications of ρ. The computation is the trajectory zn+1 = R(zn ) of a fixed dynamical system. When ρ is ReLU(t) = max(t, 0), the system is piecewise linear. When ρ has bounded range, as for the sigmoid activations used in practice, the dynamics is confined to a bounded region of state space. This leads to a concrete version of our foundational question: can primitive recursive functions be characterized as the functions computable by bounded iteration of such a recurrent network, with no explicit closure under composition as a basic external assumption, but composition emerging from the dynamics. Main results. The answer is yes, and the characterization extends to all the dynamical frameworks discussed above. Our main theorem states: ▶ Theorem 1 (Main theorem, informal1 ). For a function f : Nd → Ne , the following are equivalent: 1. f is primitive recursive in the classical sense; 2. f is computable by bounded iteration of a fixed recurrent ReLU network; 3. f is computable by bounded iteration of a fixed recurrent ρ-activation network with bounded precision, over a compact domain; 4. f is computable robustly by a fixed polynomial ODE; 5. f is computable by iteration of a fixed polynomial map in discrete time, provided a sufficiently small step size is supplied as an external parameter. In each case, “bounded” means that the relevant quantity, e.g. the number of iterations, the observation time, the domain diameter, or the precision, is itself primitive recursive as a function of the input. The activation ρ is only required to be primitive recursively computable in a natural sense, a condition that is satisfied by all sigmoid functions considered in practice. In Items (2) and (4)–(5), a crucial point is that the computation starts from the raw input state (x, 0, . . . , 0): the integer input x is embedded directly in Rm , with no auxiliary encoding. Item (3) necessarily departs from this convention, since a bounded domain cannot accommodate arbitrarily large integers as raw values; there, the input is supplied as a rational number under a fixed, standard encoding, but, in the same spirit, no problem-specific preprocessing is involved. Namely, we consider ν(n) = 2−n to encode integer n, i.e we encode n as a dyadic, but this choice is arbitrary. ▶ Remark 2. Our raw-input assumption is essential, to reveal the strength of the statement and, not merely a stylistic choice. Without it, the results would be considerably easier to prove but would say far less: the computational work could be hidden in the encoding step rather than exposed in the dynamics. As stated, our theorems guarantee that the entire computation is carried out by the dynamics of the system itself, with no external help. The equivalence (1) ⇔ (2) says that composition is redundant: bounded iteration of a single ReLU network already generates the full class. The equivalence (2) ⇔ (3) shows that any well-behaved activation ρ yields the same computational power, provided domain and precision are controlled. The passage (2) ⇒ (4) eliminates the threshold gates by moving to continuous time, where a polynomial ODE autonomously generates the required control through mechanisms like the one described above. The passage (4) ⇒ (5) discretizes the ODE, but the control that was free in continuous time now reappears as the step-size parameter: a genuine computational resource. The converse (5) ⇒ (2)
1
Formal version is Theorem 19
Olivier Bournez
states that once the map and the time bound are fixed, the induced evolution on rational states is primitive recursive. These equivalences make visible genuinely different costs of computation: the ReLU model treats threshold switching as free, the ρ-model pays with domain and precision, the polynomial ODE pays with continuous time, and the discrete polynomial model pays with step size. What these characterizations reveal The five formulations are equivalent, but they make visible genuinely different costs of computation. In the ReLU model (2), the implicit threshold gate provides exact switching for free, but this is a strong assumption on the activation. The ρ-activation model (3) shows that any reasonable nonlinearity suffices, at the price of controlling domain and precision. In the polynomial ODE (4), no discrete activation is needed at all: the dynamics itself generates rounding, switching, and error correction, but requires continuous time. In the discrete polynomial simulator (5), continuous time is also removed, and the step size becomes the price paid for its absence. Related work Our work sits at the intersection of three traditions. Classical subrecursive function theory characterizes complexity classes by restricting recursion schemes — Cobham for P [22], Bellantoni–Cook [5], Leivant–Marion for FPSPACE [29], see [21] for a survey — but always takes composition as a built-in closure. The neural-computation line [38, 37, 4, 30] showed early on that recurrent analog systems compute through state update, feedback, and thresholding rather than symbolic composition. The GPAC and polynomial-ODE tradition, from the Grzegorczyk-level analog characterizations of [19, 20, 16, 25] to the polynomial-time and polynomial-space characterizations of [14, 34, 23, 12] and the discrete-ODE lines of [9, 10, 1, 2, 3], shows that polynomial dynamics provide a machine-independent framework over the reals. Our paper combines these: machine-free like implicit complexity, polynomial dynamics like the ODE literature, discrete-time and raw-input like neither, and the internalization of composition as a consequence of the dynamics rather than a postulate.
2
The asymmetry between discrete and continuous polynomial dynamics
Computation by a polynomial dynamical system Before stating the results of this section, we make precise the common shape of computation shared by items (4) and (5) of the main theorem, and by the dynamical reading of items (2) and (3) as well. The data is a polynomial map p : Rm → Rm , viewed as a fixed vector field on the state space Rm . This single object can be evolved in two ways. In discrete time, it defines the iteration xn+1 = p(xn ). In continuous time, it defines the autonomous ODE ẋ = p(x). In both cases a computation on input x ∈ Nd is a trajectory started from the raw initial state (x, 0, . . . , 0), and the output is read off a designated coordinate at some time n (respectively t) that depends on the input. What is striking, and what this section is about, is that these two evolutions of the same polynomial data are not computationally interchangeable. Continuous time makes available mechanisms, such as contraction toward the integers and autonomous phase switching, that a discrete iteration of the same kind of object provably cannot realize uniformly. Throughout this section, f [n] := f ◦ · · · ◦ f (n times) denotes the n-fold iterate of a map f , with f [0] := id. Two primitives of dynamical computation To make the asymmetry quantitative, we examine two operations that are natural building blocks for any computation on integer inputs carried out by a dynamical system. The first is rounding: contracting a real state that is close to an integer back onto that integer, so that integer values can be maintained along a long trajectory. The second is phase selection: switching autonomously between phases of the computation after a number of steps
XX:3
XX:4
Primitive Recursion without Composition
supplied as input. Both are ubiquitous in the GPAC and polynomial-ODE literature, where they serve as the foundation on which clocks, error correction, and bounded recursion are built [24]. We now state two impossibility results showing that both of these operations, which are essentially free in continuous time for polynomial vector fields, are provably unattainable by pure polynomial maps in discrete time. The proofs are short and direct, but the consequence is conceptual: the discrete and continuous settings rest on genuinely different assumptions about what is elementary, and this explains why additional computational resources (the ρ-activation gate in the neural models, the step-size parameter in the discrete simulator) become mandatory when passing from one world to the other. ▶ Remark 3. We focus here on the contrast between discrete-time and continuous-time polynomials, but the conclusions bear directly on the neural-network setting. The choice of activation function is itself a commitment about which operations are taken as primitive. Our impossibility results make this commitment explicit: the activation ρ is not a convenience but a necessary compensation for what polynomials alone cannot do in discrete time. Since affine maps are particular polynomial maps, any impossibility for polynomial maps a fortiori rules out affine maps, and thus pure feedforward linear networks, as well. Uniform rounding One basic operation when passing between integers and reals is rounding. As observed in [24], in continuous time, this has an elegant and essentially costless solution: The map σ (x) = x −
1 sin(2πx) 2π
fixes every integer and contracts every real number toward the nearest integer: for any ε < 1/2, there exists λ < 1 such that |σ (n + δ ) − n| ≤ λ |δ | whenever |δ | ≤ ε. Since sin is natively generated in the polynomial ODE world by y′′ = −y, this contraction is available at no additional cost. This mechanism, introduced as the very first statement (Lemma 1) of [24], is the foundational building block of the entire GPAC computation theory: all subsequent control structures, such as clocks, holding phases, error correction, are built on top of it in [24]. In discrete time, the analogous task is provably impossible from a fixed polynomial map: ▶ Proposition 4 (No uniform polynomial rounder in discrete time). Let 0 < ε < δ < 1/2. There is no polynomial p : R → R and no integer N ≥ 1 such that for every m ∈ Z and every u with |u − m| ≤ δ , |p[N] (u) − m| ≤ ε. Proof. Suppose such p and N exist. The polynomial q(x) = p[N] (x) − x satisfies |q(m)| ≤ ε for every integer m. A polynomial bounded on infinitely many integers is constant: q ≡ c. But then |δ + c| ≤ ε and |−δ + c| ≤ ε, which is impossible when ε < δ . ◀ To achieve approximate rounding in discrete time, one can approximate σ on a bounded interval [−B, B] by a polynomial whose degree depends on B and on the target accuracy. An operation that is free in continuous time becomes structurally expensive in discrete time: the polynomial degree is the currency in which the discrete construction pays. Exact phase selection A second fundamental operation is switching between phases of a computation. Given a counter C supplied as input, one needs an indicator that outputs 1 for the first C steps and 0 thereafter. In continuous time, a polynomial ODE comparator achieves this autonomously, driving a register exponentially fast toward 0 or 1 depending on the sign of its input. This can be done using a polynomial map [24].
Olivier Bournez
▶ Proposition 5 (No exact polynomial phase selector). There is no polynomial map P : Rm → Rm and coordinate projection π : Rm → R, (z1 , . . . , zm ) 7→ zi , such that for every integer C ∈ N, ( 1 if n < C, π(P[n] (C, 0, . . . , 0)) = for all n ∈ N. 0 if n ≥ C, Proof. For fixed n, the function Qn (C) = π(Pn (C, 0, . . . , 0)) is a polynomial in C. Since Qn (C) = 1 for all integers C > n, the polynomial Qn − 1 has infinitely many roots, so Qn ≡ 1. But Qn (n) = 0 by hypothesis—a contradiction. ◀
3
Formal framework and equivalence for hard models
We now define the five dynamical models, state the equivalence theorems, and sketch the proof architecture. Notation Throughout the paper, an affine map from R p to Rq is a function A : R p → Rq of the form A(u) = Mu + b with M ∈ Qq×p and b ∈ Qq . Unless stated otherwise, all affine maps considered in this paper have rational coefficients. Given an ordered tuple of distinct indices I = (i1 , . . . , ie ) ∈ {1, . . . , m}e , we write πI : Rm → Re for the associated coordinate projection πI (u1 , . . . , um ) = (ui1 , . . . , uie ). Where the index tuple is clear from context, we abbreviate πI by πe .
3.1
The discrete exact model: recurrent ReLU computation
▶ Definition 6 (Feedforward ReLU block). A feedforward ReLU block on Rm is a map R : Rm → Rm of the form R = AL ◦ σ ◦ AL−1 ◦ σ ◦ · · · ◦ A1 ◦ σ ◦ A0 , where L ≥ 0 is a fixed integer (the depth); the Aℓ are affine maps with rational coefficients Aℓ : Rmℓ → Rmℓ+1 , with m0 = mL+1 = m and arbitrary intermediate widths m1 , . . . , mL ∈ N; and σ denotes the coordinatewise application of ReLU(t) := max(t, 0). Equivalently, R is the function computed by a feedforward neural network with L hidden layers, rational weights and biases, and ReLU activations. The map R is piecewise affine with rational pieces. ▶ Definition 7 (Recurrent ReLU computation). A function f : Nd → Ne admits a recurrent ReLU computation if there exist an integer m ≥ d, a feedforward ReLU block R : Rm → Rm , a primitive recursive observation time T : Nd → N, and an ordered tuple of designated output-coordinate indices I ∈ {1, . . . , m}e such that, for every input x ∈ Nd , the orbit z0 (x) := (x, 0, . . . , 0) ∈ Rm , zn+1 (x) := R(zn (x)) satisfies f (x) = πI zT (x) (x) . In other words: a single fixed piecewise-affine map with rational coefficients is iterated from the raw input (x, 0, . . . , 0), and the exact value of f (x) is read off by projecting onto the output coordinates after T (x) steps.
3.2
The neural network model: recurrent ρ-computation
For bounded-range activations (sigmoid, tanh, . . . ), the compact state space cannot accommodate arbitrarily large integers as raw values. We fix the encoding ν : N → (0, 1] defined by ν(n) := 2−n , and write ν(x) := (2−x1 , . . . , 2−xd ) for a tuple x ∈ Nd . Note that ν(0) = 1 and ν(n) ≤ 12 for n ≥ 1. ▶ Definition 8 (Feedforward ρ-block). Let ρ : [0, 1] → [0, 1]. A feedforward ρ-block of width m is a map R : [0, 1]m → [0, 1]m of the form R = AL ◦ ρ ◦ AL−1 ◦ ρ ◦ · · · ◦ A1 ◦ ρ ◦ A0 , where L ≥ 0 is a fixed integer depth; the Aℓ are affine maps with rational coefficients Aℓ : Rmℓ → Rmℓ+1 , with m0 = mL+1 = m and arbitrary intermediate widths, satisfying Aℓ [0, 1]mℓ ⊆ [0, 1]mℓ+1 for every ℓ; and ρ denotes the coordinatewise application of ρ to a vector of [0, 1]mℓ .
XX:5
XX:6
Primitive Recursion without Composition
The inclusion condition on the Aℓ simply ensures that the composition is well-defined on [0, 1]m and takes values in [0, 1]m . ▶ Definition 9 (Admissible activation). A function ρ : [0, 1] → [0, 1] is admissible if: (i) ρ is computable in the sense of computable analysis [40, 28], with a primitive recursive evaluation oracle and a primitive recursive modulus of uniform continuity; (ii) there exist a feedforward ρ-block Z : [0, 1] → [0, 1] and some η ∈ (0, 18 ) such that |Z(u)| ≤ η for u ∈ [0, 58 ] and |Z(u) − 1| ≤ η for u ∈ [ 78 , 1]. Condition (i) ensures primitive recursive simulability. Condition (ii) provides a robust zerodetection mechanism: the network can distinguish ν(0) = 1 from all other codes. Both conditions hold for all standard sigmoid activations. ▶ Definition 10 (Recurrent ρ-computation). Let ρ be admissible. A function f : Nd → Ne admits a recurrent ρ-computation if there exist an integer m ≥ d, a feedforward ρ-block R : [0, 1]m → [0, 1]m , a primitive recursive observation time T : Nd → N, a primitive recursive output precision S : Nd → N with S(x) ≥ 2, and an ordered tuple I = (i1 , . . . , ie ) ∈ {1, . . . , m}e of designated output-coordinate indices, such that the orbit z0 (x) := (ν(x), 0, . . . , 0) ∈ [0, 1]m , zn+1 (x) :=R(zn (x)) satisfies, for every output index j ∈ {1, . . . , e} and every x ∈ Nd : πI (zT (x) (x)) j − ν f j (x) ≤ 2−S(x) ν f j (x) . The bound S ≥ 2 guarantees unique decodability: for every n ∈ N, the nearest ν-code to ν(n) is at distance 12 ν(n), so a relative error of 2−S ≤ 14 suffices. The precision function S(x) plays the role of a readout resource: the network and initial state are fixed once x is given, and only the required output precision varies with the input.
3.3
The continuous model: polynomial ODE computation
Here and below, a polynomial map (resp. polynomial vector field) from R p to Rq is a function whose q components are polynomials in the p input variables with rational coefficients. Thus every affine map is in particular a polynomial map. ▶ Definition 11 (Robust polynomial ODE representation). A function f : Nd → Ne admits a robust polynomial ODE representation if there exist an integer M ≥ d, a polynomial vector field F : RM → RM with rational coefficients, a primitive recursive observation time τ : Nd → N, a primitive recursive safety bound Be : Nd → N, and an ordered tuple I = (i1 , . . . , ie ) ∈ {1, . . . , M}e of designated output-coordinate indices, such that for every x ∈ Nd the solution Y = Yx of the Cauchy problem Y ′ (t) = F Y (t) , Y (0) = (x, 0, . . . , 0) ∈ RM , e satisfies (i) Y is defined on [0, τ(x) + 1]; (ii) ∥Y (t)∥∞ ≤ B(x) on this interval; (iii) πI (Y (t)) − 1 f (x) ∞ < 2 for all t ∈ [τ(x), τ(x) + 1]. The vector field F is fixed; only the observation time, the safety bound, and the choice of output indices depend on the problem. The output is recovered by coordinatewise rounding of πI (Y (τ(x))) to the nearest integer: condition (iii) guarantees that every output coordinate stays within distance 12 of the integer f j (x) throughout the entire window [τ(x), τ(x) + 1], so the reading is insensitive to the exact observation time. This is the sense in which the representation is robust.
3.4
The intermediate model: step-size-controlled polynomial simulation
▶ Definition 12 (Step-size-controlled polynomial representation). A function f : Nd → Ne admits a step-size-controlled polynomial representation if there exist an integer M ≥ d + 1, a polynomial
Olivier Bournez
XX:7
map P : RM → RM with rational coefficients, a primitive recursive precision threshold S : Nd → N, a primitive recursive observation count N : Nd × N → N, and an ordered tuple I = (i1 , . . . , ie ) ∈ {1, . . . , M}e of designated output-coordinate indices, such that for every x ∈ Nd and every precision parameter s ∈ N with s ≥ S(x), the orbit defined by Z0 := (x, 2−s , 0, . . . , 0) ∈ RM , Zn+1 := P(Zn ) satisfies πI (ZN(x,s) ) − f (x) ∞ < 12 . We call models where the orbit remains integer-valued hard (the ReLU model and the thresholdaffine normal form below) and models where it generically leaves the integers smooth (the ρ-model, the polynomial ODE, and the step-size-controlled simulator). The map P is purely polynomial: no threshold gate or transcendental activation is involved. The step size h = 2−s is supplied as part of the initial state, stored in a dedicated coordinate. We isolate the proof-level normal form through which all equivalences pass. Fix Θ : R → R with Θ(u) = 0 for u ≤ 0 and Θ(u) = 1 for u ≥ 1; its values on (0, 1) are irrelevant here, since all gate arguments will be integers. ▶ Definition 13 (Threshold-affine normal form). A function f : Nd → Ne admits a threshold-affine normal form if there exist an integer m, affine maps A0 , . . . , As : Rm → Rm and q1 , . . . , qs : Rm → R with integer coefficients, a primitive recursive time bound T : Nd → N, and output coordinates πe : Rm → Re such that, writing s
P(y) = A0 (y) + ∑ Θ(q j (y)) A j (y),
(1)
j=1
and defining y0 (x) := (x, 0, . . . , 0), yn+1 (x) := P(yn (x)), one has f (x) = πe (yT (x) (x)) for every x ∈ Nd . Since the A j and q j have integer coefficients and y0 (x) ∈ Zm , every iterate yn (x) lies in Zm and every gate argument q j (yn (x)) is an integer; in particular, Θ is never evaluated in its ambiguous region (0, 1). ▶ Theorem 14 (Main equivalence: hard models). For a function f : Nd → Ne , the following are equivalent: 1. f is primitive recursive; 2. f admits a recurrent ReLU computation; 3. f admits a threshold-affine normal form. The proof is based on the fact that primitive recursive functions are known to coincide with those computable by a LOOP program [31], and that a single LOOP instruction compiles directly into the form given by Equation (1). The interest of the theorem lies not in the proof but in what it singles out: the threshold-affine normal form of Item 3. This is the central object of the paper, and the hinge through which the smooth models will also be analyzed. More details on proof architecture. (1) ⇒ (3) (Section 5). A LOOP program computing f is compiled one instruction at a time into (1): arithmetic is carried by the affine maps A j , and zero-testing by the threshold Θ applied to an integer argument q j (y). (3) ⇒ (2) (Section 6). Each threshold gate Θ(q j (y)) is replaced by ReLU(q j (y)) − ReLU(q j (y) − 1), an expression that coincides with Θ on every integer argument. (2) ⇒ (1) (Section 7). A fixed ReLU block with rational parameters induces, on rational initial states, a trajectory whose coordinates are primitive recursive in the iteration count. A subtlety. The compilation (3) ⇒ (2) produces a ReLU block that agrees with the underlying threshold-affine map only along the integer-valued trajectories relevant to the computation, not on all real inputs. This is enough here, since the equivalence concerns functions on Nd .
XX:8
Primitive Recursion without Composition
4
Equivalence for smooth models
4.1
A dynamical systems perspective
The equivalences for hard models in Theorem 14 establish more than the bare statement. The constructions endow the threshold-affine map P of (1), together with its orbits, with four structural properties that are most naturally phrased in the language of dynamical systems. (i) The integer lattice is forward-invariant under P. The affine maps A j have integer coefficients, and Θ takes values in {0, 1} on integer arguments; hence P(Zm ) ⊆ Zm . Every trajectory issued from the raw integer initial condition y0 (x) = (x, 0, . . . , 0) therefore stays on the lattice. (ii) The reference trajectory is bounded, uniformly in the input. There exists a primitive recursive function B : Nd → N such that, for every x ∈ Nd and every n ≤ T (x), the state yn (x) lies in the cube [−B(x), B(x)]m . Combined with (i), this confines the effective dynamics to a finite subset of the lattice Zm , whose size is primitive recursively controlled by the input. (iii) A distinguished orbit carries the computational content. For each input x, the finite sequence y0 (x), . . . , yT (x) (x) is a specific trajectory of P along which the output is read off at time T (x). The computation is encoded by the family of these reference orbits, indexed by x ∈ Nd . (iv) The reference orbit is robust to perturbations of the threshold. Every gate argument q j (yn (x)) on the reference orbit is an integer; it lies in the region where Θ is fully determined (u ≤ 0 or u ≥ 1) and never enters the ambiguous interval (0, 1). The orbit is therefore insensitive to any modification of Θ that preserves its values on Z. Property (iv) is the hinge between the hard and smooth worlds: any alternative dynamics agreeing with P on Zm reproduces the computation faithfully. Formally, we have: ▶ Theorem 15 (Uniform-robust characterization). A function f : Nd → Ne is primitive recursive if and only if it admits a threshold-affine normal form in which Θ(u) = 0 for u ≤ 14 and Θ(u) = 1 for u ≥ 34 . e Proof. (⇐) Suppose f admits a threshold-affine normal form in which the threshold function Θ 1 3 1 e e satisfies Θ(u) = 0 for u ≤ 4 and Θ(u) = 1 for u ≥ 4 . Since {u : u ≤ 0} ⊆ {u : u ≤ 4 } and {u : u ≥ e in particular satisfies Θ(u) e e e is a 1} ⊆ {u : u ≥ 34 }, Θ = 0 for u ≤ 0 and Θ(u) = 1 for u ≥ 1. Hence Θ valid instance of the threshold function fixed in Definition 13, and f admits a threshold-affine normal form in the sense of that definition. By Theorem 14, f ∈ PR. (⇒) Suppose f ∈ PR. By Theorem 14, there exist an integer m, affine maps A0 , . . . , As : Rm → Rm , affine maps q1 , . . . , qs : Rm → R, a primitive recursive time bound T : Nd → N, an index tuple I ∈ {1, . . . , m}e of output coordinates, and a threshold function Θ : R → R satisfying Θ(u) = 0 for u ≤ 0 and Θ(u) = 1 for u ≥ 1, such that, writing s P(y) = A0 (y) + ∑ Θ q j (y) A j (y), j=1
and y0 (x) := (x, 0, . . . , 0), yn+1 (x) := P yn (x) , one has f (x) = πI yT (x) (x) for every x ∈ Nd . Moreover, the compilation of Theorem 14 produces the A j and q j with integer coefficients (the affine parts encode LOOP assignments, control-register updates, and Boolean combinations of one-hot program-counter bits; the q j encode zero-tests, all acting on integer-valued registers). Step 1: the reference orbit remains in Zm We show by induction on n ∈ {0, 1, . . . , T (x)} that yn (x) ∈ Zm . For n = 0, y0 (x) = (x, 0, . . . , 0) ∈ Zm since x ∈ Nd . For the inductive step, assume yn (x) ∈ Zm . The integer coefficients of A0 and each A j give A0 (yn (x)) ∈ Zm and A j (yn (x)) ∈ Zm ; the integer coefficients of q j give q j (yn (x)) ∈ Z. Since any integer k satisfies either k ≤ 0 or k ≥ 1, and Θ
Olivier Bournez
XX:9
takes value 0 on the first set and 1 on the second, Θ(q j (yn (x))) ∈ {0, 1}. Therefore s
yn+1 (x) = A0 (yn (x)) + ∑ Θ(q j (yn (x))) A j (yn (x)) ∈ Zm . {z } | {z } j=1 | ∈{0,1}
∈Zm
In particular, every gate argument q j (yn (x)) for n ∈ {0, . . . , T (x) − 1} and j ∈ {1, . . . , s} is an integer. e : R → R by Step 2: substituting a margin threshold Define Θ e Θ(u) := ReLU(2u − 12 ) − ReLU(2u − 23 ). e Direct verification shows: for u ≤ 14 , both ReLU’s vanish, so Θ(u) = 0; for u ≥ 34 , both arguments e e meets the margin condition. On integers are ≥ 0, yielding Θ(u) = (2u − 12 ) − (2u − 32 ) = 1. Hence Θ k ∈ Z: e if k ≤ 0 then k ≤ 41 , so Θ(k) = 0 = Θ(k); 3 e if k ≥ 1 then k ≥ 4 , so Θ(k) = 1 = Θ(k). e and Θ agree on Z. Consequently Θ Step 3: the orbit is unchanged Let s e q j (y) A j (y), e P(y) := A0 (y) + ∑ Θ j=1
e yn (x)). We keeping the same A j and q j , and define the new orbit ye0 (x) := (x, 0, . . . , 0), yen+1 (x) := P(e show by induction on n ∈ {0, . . . , T (x)} that yen (x) = yn (x). The base case is immediate. For the inductive step, assume yen (x) = yn (x); by Step 1, this common e j (yn (x))) = Θ(q j (yn (x))) for every j. value is in Zm , and each q j (yn (x)) is an integer. By Step 2, Θ(q Therefore s
s
e j (yn (x))) A j (yn (x)) = A0 (yn (x))+ ∑ Θ(q j (yn (x))) A j (yn (x)) = P(yn (x)) = yn+1 (x). e yn (x)) = A0 (yn (x))+ ∑ Θ(q P(e j=1
j=1
Taking n = T (x) yields yeT (x) (x) = yT (x) (x), hence πI (e yT (x) (x)) = πI (yT (x) (x)) = f (x). The tuple e is therefore a threshold-affine normal (m, A0 , . . . , As , q1 , . . . , qs , T, I), together with the threshold Θ, form for f of the required kind. ◀ A proof is given later. The idea: the margin condition is strictly stronger than the threshold condition of Definition 13, so the "if" direction is immediate from Theorem 14. For the "only if" direction, the compilation of that theorem produces integer gate arguments (property (iv)), so any e threshold agreeing with Θ on Z, in particular the margin version Θ(u) = ReLU(2u − 21 ) − ReLU(2u − 3 2 ), leaves the reference orbit unchanged. Properties (i)–(iv) motivate the smooth models of the next subsection: it is enough for the smooth dynamics to send each integer state into the basin of correct decoding of its successor, on the bounded trapping region.
4.2
Integer-faithful shadowing
We formalize requied concepts using a notion of shadowing. ▶ Definition 16 (Integer-faithful approximation). Let D ⊆ Rm , let F, Fe : D → Rm , and assume e −F(y)∥∞ < 1 for every y ∈ D ∩ Zm . F(D ∩Zm ) ⊆ Zm . The map Fe is integer-faithful to F on D if ∥F(y) 2
XX:10
Primitive Recursion without Composition
The map Fe need not preserve the lattice or coincide with F on non-integer points. It only needs to land in the correct decoding basin on the integer skeleton of D. ▶ Lemma 17 (Decoded shadowing). Under the hypotheses above, if additionally F(D ∩ Zm ) ⊆ e yn )) started from ye0 = y0 ∈ D ∩ Zm satisfies yen = yn D ∩ Zm , then the decoded orbit yen+1 = round(F(e for all n such that y0 , . . . , yn ∈ D. e n ) lies within 1 of F(yn ) ∈ Zm , so rounding recovers Proof. Induction: if yen = yn ∈ D ∩ Zm , then F(y 2 yn+1 . ◀ For the smooth models, that do not round exactly between steps, the integer-faithful property alone (Definition 16) is no longer sufficient. What we need in addition is a quantitative control on how a one-step error propagates along the orbit. The following discrete analogue of Grönwall’s inequality provides it. ▶ Lemma 18 (Discrete Grönwall shadowing). Let D ⊆ Rm and let P : D → D admit an orbit y0 , y1 , . . . , yT (so yn+1 = P(yn ) for 0 ≤ n ≤ T − 1). Let L ≥ 0, ε ≥ 0, r > 0 satisfy the tube condition S −1 −1 k ε ∑Tk=0 L ≤ r, and let Φ be defined at least on Tn=0 B∞ (yn , r) and satisfy, for every n ∈ {0, . . . , T −1} and every z ∈ B∞ (yn , r), Φ(z) − P(yn ) ∞ ≤ L ∥z − yn ∥∞ + ε. Then the orbit defined by z0 := y0 and zn+1 := Φ(zn ) is well-defined for 0 ≤ n ≤ T and satisfies k ∥zn − yn ∥∞ ≤ ε ∑n−1 k=0 L (0 ≤ n ≤ T ), where the empty sum at n = 0 is taken to be zero. Proof sketch. Set en := ∥zn − yn ∥∞ . Then e0 = 0, and the one-step bound applied to z = zn gives the k linear-plus-constant recurrence en+1 ≤ L en + ε. Induction yields en ≤ ε ∑n−1 k=0 L ; the tube condition guarantees en ≤ r throughout, so zn never leaves the domain where the inequality holds. ◀ Proof. We repeat the statement, and involved conditions: Let D ⊆ Rm and let P : D → D admit an orbit y0 , y1 , . . . , yT (so yn+1 = P(yn ) for 0 ≤ n ≤ T − 1). Let L ≥ 0, ε ≥ 0, r > 0 satisfy the tube condition T −1
ε ∑ Lk ≤ r,
(2)
k=0
and let Φ be defined at least on z ∈ B∞ (yn , r),
ST −1
n=0 B∞ (yn , r) and satisfy, for every n ∈ {0, . . . , T − 1} and every
Φ(z) − P(yn ) ∞ ≤ L ∥z − yn ∥∞ + ε.
(3)
Then the orbit defined by z0 := y0 and zn+1 := Φ(zn ) is well-defined for 0 ≤ n ≤ T and satisfies n−1
∥zn − yn ∥∞ ≤ ε ∑ Lk (0 ≤ n ≤ T ),
(4)
k=0
where the empty sum at n = 0 is taken to be zero. We now go to the proof: Set en := ∥zn − yn ∥∞ for 0 ≤ n ≤ T . We prove by induction on n ∈ {0, . . . , T } the joint statement n−1
(a)
en ≤ ε ∑ Lk ,
(b)
en ≤ r.
(5)
k=0
Base case n = 0. Since z0 = y0 , we have e0 = 0. The empty sum at n = 0 is zero, so (a) holds with equality. Since r > 0, (b) holds as well.
Olivier Bournez
XX:11
Inductive step. Suppose (5) holds at some index n ∈ {0, . . . , T − 1}, and let us prove it at index n + 1. By the induction hypothesis (b), zn ∈ B∞ (yn , r), so Φ(zn ) is defined and (3) applies at z = zn : ∥Φ(zn ) − P(yn )∥∞ ≤ L ∥zn − yn ∥∞ + ε = L en + ε. Since zn+1 = Φ(zn ) and yn+1 = P(yn ), this gives the one-step error bound en+1 ≤ L en + ε.
(6)
Combining (6) with the induction hypothesis (a), and using L ≥ 0, n−1 n n en+1 ≤ L ε ∑ Lk + ε = ε ∑ Lk + ε = ε ∑ Lk , k=0
k=1
k=0
which is (a) at index n + 1. j For (b) at index n + 1: the partial sums ∑k=0 Lk are non-negative and non-decreasing in j (as L ≥ 0), so n
T −1
en+1 ≤ ε ∑ Lk ≤ ε ∑ Lk ≤ r k=0
k=0
by the tube condition (2). This closes the induction. The conclusion (4) is precisely (a) for all n ∈ {0, . . . , T }.
4.3
◀
Main theorem: all models
▶ Theorem 19 (Main equivalence: hard and smooth models). For a function f : Nd → Ne , the following are equivalent: 1. f is primitive recursive; 2. f admits a recurrent ReLU computation; 3. f admits a recurrent ρ-computation; 4. f admits a uniform threshold-affine normal form; 5. f admits a polynomial ODE representation; 6. f admits a step-size-controlled polynomial representation. Proof architecture By Theorems 14 and 15, it suffices to pass from the uniform-robust hard normal form to the smooth models and to close the loop back to primitive recursion. The guiding principle is that the hard dynamics comes with a reference orbit and a uniform margin, so each smooth model only has to stay in the correct decoding basin along that orbit. (4) ⇒ (5) (Section 8). Each threshold becomes a polynomial ODE comparator, and each affine update is realized by continuous-time transport in a sample-and-hold cycle; correctness follows from Lemma 18 on the trapping region. (4) ⇒ (3) (Section 9). The hard computation is transported to the ν-coded domain, and each threshold is replaced by a smooth gate built from ρ; the margin condition ensures that the smooth gate takes the same branch decisions along the reference orbit. The remaining three implications follow standard arguments and are based on the following: (5) ⇒ (6) (Section 10) is Euler discretization on the primitive-recursively bounded time-space region; (6) ⇒ (1) (Section 11.1) uses that bounded iteration of a fixed polynomial map with primitive recursive parameters is itself primitive recursive; and (3) ⇒ (1) (Section 11.2) combines effective ρ-approximation with primitive recursive simulation of the recurrent dynamics.
XX:12
Primitive Recursion without Composition
5
From primitive recursion to threshold-affine normal form
This section proves (1) ⇒ (3) of Theorem 14: every primitive recursive function admits a thresholdaffine normal form. The proof uses the standard characterization of primitive recursive functions by LOOP programs [31]. A LOOP program uses registers R1 , . . . , Rs with atomic instructions Ri := 0, Ri := R j , Ri := R j + 1, composed sequentially and by bounded loops for Ri do Q. A function is primitive recursive if and only if it is computed by a fixed LOOP program, and the running time of a fixed LOOP program is itself primitive recursive [31, 21]. From Θ , and its properties, we can derive an exact zero-test: Zero(u) := 1 − Θ(u) − Θ(−u), which satisfies Zero(u) = 1 iff u = 0 for u ∈ Z. Compiling one step into a threshold-affine map. Fix a LOOP program Π computing f : Nd → Ne . Since Π is fixed, the number of registers, instructions, and the maximal nesting depth are all constants. We encode a complete configuration by a vector y ∈ Nm containing the data registers, a one-hot program counter, a finite loop-counter stack, and a few auxiliary registers. The initial configuration for input x is (x, 0, . . . , 0). Each instruction step updates the configuration as follows. Assignments (Ri := 0, Ri := R j , Ri := R j +1) are affine. Loop entry and exit require testing whether a counter is zero or positive, which is expressed exactly by Zero and Θ on integer-valued registers. The program counter, stored in one-hot Boolean form, is combined using affine operations and Θ (for instance, a ∧ b = Θ(a + b − 1)). Since all quantities are fixed and finite, the global one-step transition is P(y) = A0 (y) + ∑sj=1 Θ(q j (y)) A j (y): a threshold-affine normal form. Notice that the nonlinearity of the construction comes entirely from the threshold tests on integer control registers; the arithmetic updates themselves are affine.
6
From threshold-affine normal form to recurrent ReLU computation
This section proves the implication (3) ⇒ (2) of Theorem 14. A fixed ReLU realization of the threshold gate. Recall that ReLU(t) := max(t, 0). Define the associated gate gadget Gate(t) := ReLU(t) − ReLU(t − 1). The key observation is that Gate behaves exactly like the proof-level threshold primitive on integer-valued inputs. ▶ Lemma 20. For every integer u, one has Gate(u) = 0 if u ≤ 0 and Gate(u) = 1 if u ≥ 1. In particular, if Θ is any threshold function satisfying Θ(u) = 0 for u ≤ 0 and Θ(u) = 1 for u ≥ 1, then Gate(u) = Θ(u) for every u ∈ Z. Proof. If u ≤ 0, then both ReLU(u) and ReLU(u − 1) are zero, so Gate(u) = 0. If u ≥ 1, then ReLU(u) = u and ReLU(u − 1) = u − 1, hence Gate(u) = 1. ◀ Thus, whenever the arguments of the threshold tests are integer-valued, the threshold primitive can be replaced exactly by the fixed two-ReLU gadget Gate. Compiling one threshold-affine normal form into a ReLU block. Let P : Rm → Rm be a map of the form P(y) = A0 (y) + ∑sj=1 Θ(q j (y)) A j (y), where the A j are affine maps and the q j are affine maps. Assume that, along the trajectory from every raw input (x, 0, . . . , 0) relevant to the computation, all values q j (yn ) are integers. We first compute the affine quantities q j (y), then apply the ReLU gadget Gate to each of them, and finally use the resulting gate values in the affine recombination A0 (y) + ∑ j Gate(q j (y)) A j (y).
Olivier Bournez
Since each Gate(q j (y)) is exactly equal to Θ(q j (y)) on the intended trajectory, the resulting ReLU block reproduces exactly the same one-step update on that trajectory. Formally, we obtain the following. ▶ Lemma 21. Let P : Rm → Rm be a threshold-affine normal form map of the form P(y) = A0 (y) + ∑sj=1 Θ(q j (y)) A j (y). Assume that, along the orbit from every raw input (x, 0, . . . , 0) with x ∈ Nd , all gate arguments q j (yn ) are integers. Then there exist an integer M ≥ m, a feedforward ReLU block R : RM → RM , an embedding ι : Rm → RM , and a projection π : RM → Rm such that, for every such raw input x and every n ∈ N, one has π(Rn (ι(x, 0, . . . , 0))) = Pn (x, 0, . . . , 0). Proof. Replace each occurrence of Θ(q j (y)) in the expression of P by the ReLU gadget Gate(q j (y)) = ReLU(q j (y))−ReLU(q j (y)−1). By Lemma 20, this replacement is exact on the intended trajectories, because all gate arguments are integers there. The resulting one-step update is obtained from affine maps and coordinatewise applications of ReLU, hence is computed by a fixed feedforward ReLU block, possibly after enlarging the state by finitely many auxiliary coordinates that store intermediate affine quantities and gate values. Since the construction is local and uniform, the same block R works for all inputs. By construction, its iterates agree exactly with those of P on the embedded raw-input trajectories. ◀ Application to the threshold-affine normal form. We now apply the previous lemma to the prooflevel normal form of Section 5. The crucial hypothesis of Lemma 21 is automatically satisfied there: the gate arguments arise from integer-valued counters, one-hot instruction pointers, and loop-stack registers, so they remain integers throughout the execution. ▶ Theorem 22. Every threshold-affine normal form unfolding admits an exact recurrent ReLU computation. Equivalently, Item (3) implies Item (2) in Theorem 14. Proof. Let f : Nd → Ne admit an threshold-affine normal form. By Definition, there exist an integer m, a threshold-affine normal form map P : Rm → Rm , a primitive recursive time bound T : Nd → N, and designated output coordinates πe such that f (x) = πe (PT (x) (x, 0, . . . , 0)) for every input x. By the construction of Section 5, the gate arguments along the relevant raw-input trajectories are integer-valued. Hence Lemma 21 applies. We obtain a feedforward ReLU block R, an embedding ι, and a projection π such that the iterates of R agree exactly with those of P on the embedded raw-input trajectories. Composing the projections if necessary, we conclude that f (x) = πe′ (RT (x) (ι(x, 0, . . . , 0))) for some designated output coordinates πe′ . Thus f admits an exact recurrent ReLU computation in the sense of Definition 7. ◀ ▶ Remark 23. The theorem should not be read as saying that the discontinuous threshold-affine normal form transition map and the ReLU block coincide as functions on all of Rm . What is preserved is the exact orbit on the raw integer inputs relevant to the computation. This is exactly what Item (2) of the main theorem requires.
7
From recurrent ReLU computation back to primitive recursion.
We next prove the converse implication (2) ⇒ (1) of Theorem 14. This is elementary, because a fixed feedforward ReLU block with rational parameters induces a primitive recursive evolution on rational states. ▶ Theorem 24. Every exact recurrent ReLU computation computes a primitive recursive function. Equivalently, Item (2) implies Item (1) in Theorem 14.
XX:13
XX:14
Primitive Recursion without Composition
Proof. Assume that f : Nd → Ne admits a recurrent ReLU computation. Thus there exist a feedforward ReLU block R : Rm → Rm , a primitive recursive observation time T : Nd → N, and designated output coordinates πe such that, for every input x ∈ Nd , if z0 (x) := (x, 0, . . . , 0) and zn+1 (x) := R(zn (x)), then f (x) = πe (zT (x) (x)). A feedforward ReLU block is built from affine maps with rational parameters and coordinatewise applications of ReLU(t) = max(t, 0). On rational inputs, affine maps with rational parameters are primitive recursive, and ReLU is primitive recursive as well, since it is just the maximum of its input and 0. It follows that the one-step update (x, n) 7→ zn+1 (x) is primitive recursive on rational encodings. Hence, by primitive recursion on n, the full evolution (x, n) 7→ zn (x) is primitive recursive. Since T is primitive recursive and πe is a projection, x 7→ πe (zT (x) (x)) is primitive recursive. By assumption, this is exactly f (x). ◀
8
From uniform threshold-affine normal form to polynomial ODE
We now prove (4) ⇒ (5) of Theorem 19. The strategy Let P(y) = A0 (y) + ∑sj=1 Θ(q j (y)) A j (y) be a uniform threshold-affine normal form with reference orbit y0 (x), . . . , yT (x) (x) ⊆ Zm ∩ [−B(x), B(x)]m , integer gate arguments (property (iv) of Section 4.1), and output projection πI . Our goal is a polynomial vector field FP whose flow samples, at a fixed period κ0 > 0, a sequence of states zn that tracks yn (x) to within a constant fraction of the unit spacing of Zm . This is exactly the sort of tracking that Section 4.2 calls integer-faithful shadowing, applied at a resolution that allows decoding by rounding to the nearest integer. Concretely, we will build FP so that the sampled map zn 7→ zn+1 agrees with P up to a contraction-plus-noise error of the form ∥zn+1 − P(zn )∥∞ ≤ λ ∥zn − yn (x)∥∞ + η, with λ < 1 and η/(1−λ ) < 14 . Since P is integer-faithful to itself, the discrete Grönwall shadowing lemma (Lemma 18) then keeps the sampled orbit within the 14 -decoding basin of the reference orbit throughout time. The substance of the section is that two standard polynomial-ODE gadgets suffice to produce such a sampled map: a smooth clock, which organizes the cycle, and a polynomial comparator, which realizes Θ up to the η-noise tolerated by the shadowing lemma. ▶ Remark 25. The resulting construction can be read as a particular improved instance of the Branicky [18] alternating-targeting trick used, for example, in [24, 11]. We provide an original view on it. Rather than analyzing the continuous flow step by step, we treat it as a candidate integer-faithful approximation of P in the sense of Section 4.2 and let the discrete Grönwall lemma handle the propagation of error. The one computation that remains specific to the continuous construction is the per-cycle estimate (7); everything else is a direct application of the shadowing framework. The smooth clock Let (c1 , c2 ) satisfy the rational-coefficient ODE c′1 = −c2 , c′2 = c1 with initial condition (c1 (0), c2 (0)) = (1, 0); its solution (c1 (t), c2 (t)) = (cost, sint) is periodic of period κ0 := 2π. A polynomial driver h(t), computed as a fixed polynomial expression in (c1 , c2 ), can be arranged to be non-negative, to satisfy h(t) ≥ 1 on each compute half-period [nκ0 , nκ0 + κ0 /2], and to vanish on each hold half-period [nκ0 + κ0 /2, (n+1)κ0 ], for every integer n ≥ 0 (see [11, §3] for a standard realization). Everything below is phrased as an ODE whose right-hand side multiplies h(t): active during compute phases, frozen during hold phases. The polynomial comparator For every δ ∈ (0, 14 ) and every contraction rate Λ > 0, there is a polynomial vector field on an extra coordinate r such that for every input u ∈ R with u ∈ / ( 14 , 34 ) ′ and every compute half-period [t0 ,t0 + κ0 /2], the driven equation r = h(t) · Λ · (Θ∞ (u) − r) satisfies |r(t0 + κ0 /2) − Θ∞ (u)| ≤ e−Λκ0 /2 |r(t0 ) − Θ∞ (u)|, which is below δ once Λ is chosen large enough. Here Θ∞ (u) ∈ {0, 1} is the target value determined by whether u ≤ 14 or u ≥ 34 ; it is computed by a
Olivier Bournez
XX:15
polynomial in u approximating Θ on the two separated regions (details in [11, §3]). In our application u will be one of the gate arguments q j (π(Z)), tracked continuously by an affine expression in the state. Property (iv) of Section 4.1 guarantees that every comparator call on the reference orbit satisfies the margin hypothesis q j (yn (x)) ∈ / ( 14 , 43 ), hence falls in a determined regime. Composing the step Given a cycle-start state zn ∈ Rm close to yn (x), we want a cycle-end state zn+1 close to P(yn (x)). Per clock period κ0 we allocate s comparator coordinates r1 , . . . , rs driven by the comparator ODE on u j := q j (π(Z)) (so r j → Θ(q j (yn (x))) exponentially during the compute phase); m target-tracking coordinates w1 , . .. , wm driven by the sample-and-hold equation w′k = h(t) · Λ · (A0 (π(Z)) + ∑ j r j A j (π(Z)))k − wk ; and a projection π that reads w(t) at the end of each hold phase. Assembling these ODEs with the clock gives a polynomial vector field FP with rational coefficients on RM , M := 2 + s + m, with initial condition (1, 0, 0, . . . , 0, x, 0, . . . , 0). An explicit oneperiod calculation yields, with δ := e−Λκ0 /2 , λ := δ · max j ∥A j ∥∞ , and η := δ (1 + s · max j ∥A j ∥∞ ), the estimate ∥zn+1 − P(yn (x))∥∞ ≤ λ ∥zn − yn (x)∥∞ + η, zn := w(nκ0 ). (7) Choosing Λ large enough makes λ < 1 and η/(1−λ ) < 41 . ▶ Theorem 26. Every uniform threshold-affine normal form admits a robust polynomial ODE representation. Proof sketch. Set en := ∥π(Z(nκ0 )) − yn (x)∥∞ . By (7), e0 = 0 and en+1 ≤ λ en + η; Lemma 18 with L := λ , ε := η, r := 14 gives en < 14 for every n ≤ T (x), so the sampled state is integer-faithful e := RP (B(x)+1) (a polynomial Grönwall to the reference orbit. Setting τ(x) := ⌈κ0 T (x)⌉ and B(x) majorant for FP on one cycle), the interval [τ(x), τ(x)+1] lies inside the final hold half-period, on which w is frozen within 14 of yT (x) (x). Projecting by πI yields ∥πI (Z(t)) − f (x)∥∞ < 14 throughout [τ(x), τ(x)+1], which is Definition 11. ◀ Proof. Let FP : RM → RM be the polynomial vector field assembled in Section 8, with state coordinates (c1 , c2 , r1 , . . . , rs , w1 , . . . , wm ) and M = 2 + s + m. Let π : RM → Rm be the projection (c, r, w) 7→ w, let πI : Rm → Re be the output-coordinate projection of the normal form P, and write Π := πI ◦ π for the composite output projection. Let λ ∈ [0, 1) and η > 0 be as in (7), with η/(1−λ ) < 14 secured by the contraction rate Λ. The initial state is Z(0) := (1, 0, 0, . . . , 0, x, 0, . . . , 0): the clock at phase zero, the comparators at zero, and the w-block at y0 (x) = (x, 0, . . . , 0). Set the sampling times tn := nκ0 , and the sampled error en := ∥π(Z(tn )) − yn (x)∥∞ whenever Z(tn ) is defined. Step 1: existence of the flow and sampled error We prove by induction on n ∈ {0, 1, . . . , T (x)} the joint statement n−1
(a) Z is defined on [0,tn +κ0 ];
(b) en ≤ η ∑ λ k ;
(c) en ≤ 14 .
(8)
k=0
Base case n = 0. FP admits a local solution from Z(0) that extends to [0, κ0 ] (local existence for polynomial ODEs, with no escape to infinity on one cycle), so (a) holds. Since π(Z(0)) = y0 (x), we have e0 = 0; the empty sum at n = 0 equals zero, so (b) holds; and (c) holds trivially. Inductive step. Assume (8) at some n < T (x). By (c), ∥π(Z(tn )) − yn (x)∥∞ ≤ 41 , and by property (iv), q j (yn (x)) ∈ Z for every j, so a fortiori q j (π(Z(tn ))) ∈ / ( 14 , 43 ), the comparator margin hypothesis is satisfied, and the polynomial flow extends to [tn ,tn+1 +κ0 ], giving (a) at n+1. The cycle-to-cycle estimate (7) gives en+1 ≤ λ en + η; combined with (b) and λ ≥ 0, en+1 ≤ η ∑nk=0 λ k , giving (b) at n+1. For (c), en+1 ≤ η/(1−λ ) < 14 by the choice of λ and η. This closes the induction.
XX:16
Primitive Recursion without Composition
Step 2: the output window Set τ(x) := ⌈κ0 T (x)⌉, primitive recursive in x since κ0 = 2π is a primitive recursive real and T is primitive recursive. Since κ0 T (x) ≤ τ(x) ≤ κ0 T (x) + 1 and κ0 /2 = π > 2, the interval [τ(x), τ(x)+1] is contained in the final hold half-period [κ0 T (x), κ0 T (x)+κ0 /2]. On this halfperiod, h(t) = 0 freezes every w-coordinate, so π(Z(t)) = π(Z(tT (x) )) throughout. By Step 1 at n = T (x), ∥π(Z(tT (x) )) − yT (x) (x)∥∞ ≤ 14 , hence ∥π(Z(t)) − yT (x) (x)∥∞ ≤ 14 on [τ(x), τ(x) + 1]. Projecting by πI and using πI (yT (x) (x)) = f (x) gives ∥Π(Z(t)) − f (x)∥∞ ≤ 41 < 12 , which is condition (iii) of Definition 11. Step 3: the safety bound Since FP is a fixed polynomial vector field, a standard Grönwall argument gives a fixed polynomial RP : N → N such that every solution Z of Z ′ = FP (Z) on a time interval of length at most κ0 + 1 satisfies ∥Z(t)∥∞ ≤ RP (∥Z(t0 )∥∞ + 1). Applying this cycle by cycle and using ∥yn (x)∥∞ ≤ B(x) together with (c) of Step 1, every sampled state satisfies ∥π(Z(tn ))∥∞ ≤ B(x) + 1; (T (x)+1) iterating the bound T (x) + 1 times yields ∥Z(t)∥∞ ≤ RP (B(x) + 1) on [0, τ(x) + 1]. Setting (T (x)+1) e B(x) := RP (B(x) + 1), which is primitive recursive in x (as an iterated composition of a fixed polynomial with primitive recursive arguments), gives condition (ii). e and πI as defined, so f Conditions (i), (ii), (iii) of Definition 11 are thus all satisfied with FP , τ, B, admits a robust polynomial ODE representation. ◀
9
From uniform threshold-affine normal form to recurrent ρ-computation
We now prove (4) ⇒ (3) of Theorem 19. Let P(y) = A0 (y) + ∑sj=1 Θ(q j (y)) A j (y) be a uniform threshold-affine normal form on Rm , with A j and q j of integer coefficients (as produced by Theorem 14), reference orbit y0 (x), . . . , yT (x) (x) ⊆ Zm ∩ [−B(x), B(x)]m (properties (i)–(ii) of Section 4.1), and output coordinates πI . We use a strategy with similarities with previous section, except that bassins are not anymore corresponding to integers, but to encoded integers. Encoding and basins Encode integers by ν(n) = 2−|n| , with the sign carried on an auxiliary coordinate (suppressed in the notation below; the construction is symmetric). For y ∈ Zm write ⌊y⌋ν := (ν(y1 ), . . . , ν(ym )) ∈ [0, 1]m , and define the ε-basin around ⌊y⌋ν as the product of relative m intervals Bε (y) := ∏m k=1 [(1−ε)ν(yk ), (1+ε)ν(yk )] ⊆ [0, 1] . The relative tolerance is the natural 1 one for dyadic codes, and ε < 8 ensures distinct states have disjoint basins. The coded one-step update The three ingredients of P transport to the coded domain as follows. Integer-coefficient affine maps transport to feedforward ρ-blocks with dyadic-rational coefficients (an integer multiplier c acts on ν(n) by ν(n) 7→ 2−c ν(n)). The threshold Θ(q j (y)) ∈ {0, 1} is replaced by Z(q̌ j (z)), where Z is the admissible detector of Definition 9 and q̌ j is the coded form of q j ; property (iv) of Section 4.1 guarantees that ν(q j (yn )) ∈ [0, 85 ] ∪ [ 78 , 1] throughout the reference orbit, the region where Z is unambiguous. Writing Ǎ j for the coded form of A j , the global coded one-step update is s R(z) = Ǎ0 (z) + ∑ Z q̌ j (z) · Ǎ j (z), (9) j=1
which is a feedforward ρ-block as required by Definition 8. ▶ Lemma 27 (One robust smooth coded step). There exists ε0 ∈ (0, 81 ), depending only on the normal form and on the detector margin η, such that for every ε ≤ ε0 and every y ∈ Zm with all gate arguments q j (y) integer, R(Bε (y)) ⊆ Bε (P(y)). Proof. We prove the inclusion by controlling, step by step, the relative error between the coded output R(z) and the target ⌊P(y)⌋ν . Throughout the proof, constants CA ,Cq ,CR below depend only
Olivier Bournez
XX:17
on the fixed normal form (the number s of gates, the dimension m, and the integer coefficients of A0 , . . . , As and q1 , . . . , qs ) and not on y, z, or ε. Convention on the sign coordinate Recall that ν is extended to Z by ν(−n) := 2−n , with an auxiliary sign bit σk ∈ {0, 1} per state coordinate indicating whether yk < 0. The basin Bε (y) fixes every sign bit to its correct value, i.e. the signs are absolute (not relative) coordinates. Throughout the proof, we fix ε0 < 14 , which is small enough to ensure that sign bits are recovered unambiguously from any z ∈ Bε (y), and we treat the sign bookkeeping silently in the coded affine and gate transports. Step 1: relative error is preserved by coded affine maps Fix an integer affine map A : Z p → Z, A(u) = ∑i ci ui + b with ci , b ∈ Z. Its coded form Ǎ : [0, 1] p → [0, 1] is the map Ǎ(ν(u1 ), . . . , ν(u p )) = ν A(u1 , . . . , u p ) = 2−b ∏i 2−ci ui = 2−b ∏i ν(ui )ci when all ui , A(u) ≥ 0 (the negative case is handled by the sign-bit bookkeeping, symmetrically). We claim that if (ũ1 , . . . , ũ p ) ∈ R p satisfies |ũi − ν(ui )| ≤ ε ν(ui ) for every i, then Ǎ(ũ) − ν(A(u)) ≤ CA ε ν(A(u)),
(10)
for a constant CA depending only on A. Indeed, writing ũi = (1 + δi ) ν(ui ) with |δi | ≤ ε, c Ǎ(ũ) = 2−b ∏ (1 + δi )ν(ui ) i = ν(A(u)) · ∏(1 + δi )ci . i
i
For ε ≤ 14 and |ci | ≤ C (a bound fixed by the normal form), the product satisfies | ∏i (1 + δi )ci − 1| ≤ CA ε for CA := 2pC (using |1 − (1 + δ )c | ≤ 2|c| |δ | when |cδ | ≤ 12 ), which gives (10). Applying this componentwise to each A j : Zm → Zm yields: for every z ∈ Bε (y), Ǎ j (z) ∈ BCA ε A j (y) , j = 0, . . . , s. (11) Step 2: detector output is close to the true threshold The same argument applied to the coded gate argument q̌ j , which is a feedforward ρ-block encoding the same integer affine map, gives a constant Cq such that q̌ j (z) − ν(q j (y)) ≤ Cq ε ν(q j (y))
for z ∈ Bε (y).
(12)
By hypothesis, q j (y) is an integer, so ν(q j (y)) ∈ {1} (when q j (y) = 0) or ν(q j (y)) ∈ [0, 21 ] (when q j (y) ̸= 0). Case 1: q j (y) = 0, so Θ(q j (y)) = 1. Then ν(q j (y)) = 1 and (12) gives q̌ j (z) ≥ 1−Cq ε. Choosing ε0 small enough that 1 −Cq ε0 ≥ 78 (any ε0 ≤ 8C1 q works), the detector hypothesis gives Z(q̌ j (z)) ∈ [1 − η, 1], i.e. Z(q̌ j (z)) − Θ(q j (y)) ≤ η.
(13)
Case 2: q j (y) ̸= 0, so Θ(q j (y)) = 0. Then ν(q j (y)) ≤ 12 , and (12) gives q̌ j (z) ≤ 12 (1 + Cq ε). Choosing ε0 further small enough that 12 (1 +Cq ε0 ) ≤ 85 (any ε0 ≤ 4C1 q works), the detector hypothesis gives Z(q̌ j (z)) ∈ [0, η], so (13) holds in this case as well. Thus in both cases, Z(q̌ j (z)) − Θ(q j (y)) ≤ η,
j = 1, . . . , s.
(14)
Step 3: assembling the global one-step update Fix z ∈ Bε (y) and set, coordinatewise, the error of the coded output R(z) against the target ⌊P(y)⌋ν : i s h ∆ := R(z) − ⌊P(y)⌋ν = Ǎ0 (z) − ⌊A0 (y)⌋ν + ∑ Z(q̌ j (z)) Ǎ j (z) − Θ(q j (y)) ⌊A j (y)⌋ν . j=1
XX:18
Primitive Recursion without Composition
The first bracket is bounded by CA ε ∥⌊A0 (y)⌋ν ∥∞ by (11). For the j-th summand, write a j := Θ(q j (y)) ∈ {0, 1} and ã j := Z(q̌ j (z)); then ã j Ǎ j (z) − a j ⌊A j (y)⌋ν = (ã j − a j ) Ǎ j (z) + a j Ǎ j (z) − ⌊A j (y)⌋ν , which is bounded in ∥ · ∥∞ by η · ∥Ǎ j (z)∥∞ + CA ε ∥⌊A j (y)⌋ν ∥∞ using (14) and (11). Since Ǎ j (z) ∈ BCA ε (A j (y)) ⊆ [0, 1]m , we have ∥Ǎ j (z)∥∞ ≤ 1, and a fortiori ∥⌊A j (y)⌋ν ∥∞ ≤ 1. Summing over j gives ∥∆∥∞ ≤ (s+1)CA ε + s η =: CR ε + s η,
(15)
with CR := (s+1)CA . Step 4: the error fits inside Bε (P(y)) The target basin Bε (P(y)) is the product of relative intervals of radius ε ν(P(y)k ) around ν(P(y)k ). Since ν(P(y)k ) ≥ 2−B(x) on the reference orbit (property (ii) of Section 4.1), it suffices to check that every coordinate of ∆ has absolute value at most ε ν(P(y)k ). The reference orbit stays in [−B(x), B(x)]m , so ν(P(y)k ) ≥ 2−B(x) and the target radius is at least ε · 2−B(x) in every coordinate. By (15), we need CR ε + s η ≤ ε · 2−B(x)
coordinatewise.
(16)
Since s and CR are fixed by the normal form, this is a condition on ε alone once η is fixed by Definition 9. In particular, we may take ε0 := min 8C1 q , 4C1 q , 2C1R , Cs ηR · 2B(x) (the first two bounds come from Step 2, the third from the coefficient in (16), and the last ensures the additive s η term fits). Since the detector hypothesis η < 81 is absolute, and B(x) is bounded on the reference orbit, ε0 is bounded away from zero. For every ε ≤ ε0 , (16) holds with slack, hence ∆k ∈ −ε ν(P(y)k ), ε ν(P(y)k ) for every coordinate k, which is exactly the statement R(z) ∈ Bε (P(y)). ◀ Proof sketch. For z ∈ Bε (y), each coded affine map Ǎ j (z) lies in a basin of radius O(ε) around ⌊A j (y)⌋ν , and each coded gate argument q̌ j (z) differs from ν(q j (y)) by a relative error O(ε). By property (iv), ν(q j (y)) lies in an unambiguous region of Z, so Z(q̌ j (z)) differs from Θ(q j (y)) by at most η. Assembling these bounds, the coordinatewise error of R(z) against ⌊P(y)⌋ν is O(ε) + O(η), with constants depending only on the normal form. Since the target basin has radius proportional to ε, fitting requires ε and η to be matched: it suffices to take ε above a fixed multiple of η (to absorb the constant η-contribution) and below a fixed upper bound (to absorb the ε-contribution). Both conditions are compatible because η < 18 . ◀ ▶ Theorem 28. Every uniform threshold-affine normal form yields a recurrent ρ-computation. Proof sketch. The reference orbit stays in the cube [−B(x), B(x)]m (property (ii)) and has integer gate arguments throughout (property (iv)), so Lemma 27 applies at every step with a fixed basin radius ε. Choosing ε := 2−s0 (x) with s0 (x) primitive recursive in B(x) makes the induction zn ∈ Bε (yn (x)) go through from the raw encoded input. At n = T (x) the output coordinates approximate ν( f j (x)) within relative error 2−s0 (x) , which is Definition 10 with Tρ := T and Sρ := s0 . ◀ Proof. Let P be a uniform threshold-affine normal form for f with reference orbit y0 (x), . . . , yT (x) (x), trapping radius B : Nd → N primitive recursive (property (ii) of Section 4.1), and output projection πI . Let R : [0, 1]m → [0, 1]m be the coded one-step update of (9), and let η < 18 be the margin of the admissible detector Z (Definition 9).
Olivier Bournez
Step 1: uniform basin radius on the reference orbit Set β (x) := B(x) + max j ∥A j ∥∞ · m, where ∥A j ∥∞ is the ℓ∞ -operator norm of the integer matrix of A j . Then for every n ≤ T (x), ∥yn (x)∥∞ ≤ B(x) ≤ β (x); moreover, since yn+1 (x) = P(yn (x)) = A0 (yn (x)) + ∑ j Θ(q j (yn (x)))A j (yn (x)), every coordinate of yn+1 (x) is bounded by β (x) as well, so the reference orbit is confined to [−β (x), β (x)]m . Apply Lemma 27 with range parameter β (x): there exists ε0 (x) ∈ (0, 81 ), depending on the normal form, on η, and on β (x), such that for every ε ≤ ε0 (x) and every y ∈ Zm ∩ [−β (x), β (x)]m with integer gate arguments, R(Bε (y)) ⊆ Bε (P(y)). Tracing the formula for ε0 from the proof of Lemma 27, we have ε0 (x) ≥ 2−c(P)β (x) for a constant c(P) depending only on the normal form. Set s0 (x) := max c(P)β (x), 2 , ε(x) := 2−s0 (x) . Then ε(x) ≤ ε0 (x), and s0 (x) is primitive recursive in x (as a fixed function of B(x), itself primitive recursive). The condition s0 (x) ≥ 2 is the one required by Definition 10. Step 2: induction on the coded orbit Let z0 := (ν(x), 0, . . . , 0) ∈ [0, 1]m be the raw encoded input, and define zn+1 := R(zn ) for n ≥ 0. We prove by induction on n ∈ {0, 1, . . . , T (x)} that zn ∈ Bε(x) yn (x) . (17) Base case n = 0. We have z0 = (ν(x), 0, . . . , 0) = ⌊y0 (x)⌋ν exactly, where we use ν(0) = 1 on the zero-padding coordinates. Hence the coordinate-wise difference z0,k − ν(y0,k (x)) is zero for every k, which is trivially within relative error ε(x). So z0 ∈ Bε(x) (y0 (x)). Inductive step. Assume zn ∈ Bε(x) (yn (x)) for some n < T (x). By property (iv), every gate argument q j (yn (x)) is an integer. By property (ii) and Step 1, yn (x) ∈ Zm ∩ [−β (x), β (x)]m . By Step 1, ε(x) ≤ ε0 (x), so Lemma 27 applies: zn+1 = R(zn ) ∈ R(Bε(x) (yn (x))) ⊆ Bε(x) (P(yn (x))) = Bε(x) (yn+1 (x)), which is (17) at index n + 1. Step 3: verifying Definition 10 At n = T (x), (17) gives zT (x) ∈ Bε(x) (yT (x) (x)). Projecting onto the output coordinates I = (i1 , . . . , ie ), every output coordinate j ∈ {1, . . . , e} satisfies πI (zT (x) ) j − ν f j (x) ≤ ε(x) ν f j (x) = 2−s0 (x) ν f j (x) , using that πI (yT (x) (x)) = f (x) (the normal form computes f ) and that the basin Bε is a product of relative intervals. This is exactly the bound required by Definition 10 with observation time Tρ (x) := T (x) and output precision Sρ (x) := s0 (x), both primitive recursive in x. Hence f admits a recurrent ρ-computation. ◀
10
From robust polynomial ODE computation to step-size-controlled polynomial simulation
This section proves the implication 5) ⇒ 6) of Theorem 19. The idea is standard: once a function is computed robustly by a fixed polynomial ODE, one may discretize that ODE by the Euler scheme. The resulting discrete update is polynomial, provided the step size is stored as an additional register. The only subtlety is to make the dependence on the step size explicit and primitive recursive. The step-size-controlled model. Recall that Item 6) of Theorem 19 is the following notion. ▶ Definition 29 (Step-size-controlled polynomial representation). A function f : Nd → Ne admits a step-size-controlled polynomial representation if there exist a polynomial map P : RM → RM with rational coefficients, a primitive recursive precision threshold S : Nd → N, a primitive recursive observation count N : Nd × N → N, and designated output coordinates πe such that, for every input x ∈ Nd and every s ≥ S(x), if z0 (x, s) := (x, 2−s , 0, . . . , 0) and zn+1 (x, s) := P(zn (x, s)), then ∥πe (zN(x,s) (x, s)) − f (x)∥∞ < 21 .
XX:19
XX:20
Primitive Recursion without Composition
The Euler simulator. Let F : RM → RM be a fixed polynomial vector field. Its explicit Euler step with step size h is Y 7→ Y + hF(Y ). If h is stored as an extra register, then this becomes a polynomial self-map. ▶ Definition 30 (Euler simulator). Given a polynomial vector field F : RM → RM , define the associated Euler simulator EF : RM+1 → RM+1 by EF (Y, h) := (Y + hF(Y ), h). The map EF is polynomial with rational coefficients whenever F is. Uniform bounds for a fixed polynomial vector field. Let F : RM → RM be fixed. Since F is polynomial, there exist explicit integer-valued polynomials MF♯ , LF♯ : N → N such that, for every R ∈ N and every Y, Z ∈ [−R, R]M , one has ∥F(Y )∥∞ ≤ MF♯ (R) and ∥F(Y ) − F(Z)∥∞ ≤ LF♯ (R) ∥Y − Z∥∞ . We shall use these as coarse majorants for the size and local Lipschitz constant of the vector field on bounded regions. ▶ Lemma 31 (Global Euler error bound on a bounded region). Let F : RM → RM be a fixed polynomial vector field, and let Yx : [0, T ] → RM be a solution of Y ′ (t) = F(Y (t)) such that ∥Yx (t)∥∞ ≤ R for all t ∈ [0, T ]. Let Y0 := Yx (0) and Yk+1 := Yk + hF(Yk ) be the explicit Euler orbit with step size h > 0. Then, as long as all Euler iterates remain in [−R, R]M , the global error ek := ∥Yk −Yx (kh)∥∞ ♯ satisfies ek ≤ C(R, T ) h for all kh ≤ T , where C(R, T ) := LF♯ (R)MF♯ (R) T eLF (R)T . Proof. For one Euler step, the standard local truncation estimate gives ∥Yx ((k + 1)h) − Yx (kh) − hF(Yx (kh))∥∞ ≤ LF♯ (R)MF♯ (R)h2 , because on the relevant region the vector field has size at most MF♯ (R) and Lipschitz constant at most LF♯ (R). Writing ek := ∥Yk −Yx (kh)∥∞ , one gets ek+1 ≤ (1 + hLF♯ (R))ek + LF♯ (R)MF♯ (R)h2 . Since e0 = 0, a discrete Grönwall estimate yields ek ≤ C(R, T )h for all kh ≤ T , with the stated constant. ◀ Discretizing the robust ODE computation. We now apply the previous lemma to the robust polynomial ODE representation of Item 5). ▶ Theorem 32. Every robust polynomial ODE representation yields a step-size-controlled polynomial representation. Equivalently, Item 5) implies Item 6) in Theorem 19. Proof. Assume that f : Nd → Ne admits a robust polynomial ODE representation. Thus there exist a polynomial vector field F : RM → RM , a primitive recursive observation time τ : Nd → N, a primitive recursive safety bound Be : Nd → N, and designated output coordinates πe such that the solution of e e M for all Y ′ (t) = F(Y (t)) with initial condition Y (0) = (x, 0, . . . , 0) remains in the box [−B(x), B(x)] 1 t ∈ [0, τ(x) + 1], and moreover ∥πe (Yx (t)) − f (x)∥∞ < 4 for all t ∈ [τ(x), τ(x) + 1]. As usual, the margin 14 is obtained simply by tightening the precision in the ODE construction; only a margin strictly below 12 is needed. e + 1, Let EF be the Euler simulator of Definition 30. For a given input x, set R(x) := B(x) ♯ ♯ L (τ(x)+1) x Lx := LF (R(x)), Mx := MF (R(x)), and Cx := Lx Mx (τ(x) + 1)e . Choose a primitive recursive threshold S : Nd → N such that 2−S(x)Cx < 14 for every input x. This is possible because all quantities involved are primitive recursive. For every s ≥ S(x), let h := 2−s and define N(x, s) := ⌈τ(x) 2s ⌉. Since 2s and τ(x) are primitive recursive, so is N. We now iterate the Euler simulator from the initial state z0 (x, s) := ((x, 0, . . . , 0), 2−s ). By Lemma 31, for every s ≥ S(x) the discrete state of the Euler orbit after N(x, s) steps is within distance 14 of the exact ODE state at time N(x, s)2−s . Since |N(x, s)2−s − τ(x)| ≤ 2−s ≤ 1, this time belongs to the final validity interval [τ(x), τ(x) + 1]. On that interval, the ODE output is already within distance 14 of f (x). Therefore the total error is strictly less than 12 . This proves that the Euler simulator provides a step-size-controlled polynomial representation of f. ◀
Olivier Bournez
▶ Remark 33. The important point is that the step size is not a hidden numerical detail, but an explicit computational resource. The continuous-time polynomial ODE uses a fixed flow. After discretization, the same mechanism survives only provided enough numerical resolution is supplied. This is precisely what Item 6) makes visible.
11
Converse implications and completion of the proof
It remains to prove the easy converse implications. These are much more direct than the forward constructions, since in both cases the dynamics is fixed in advance and the evolution takes place over rational states.
11.1
From step-size-controlled polynomial simulation back to primitive recursion.
We prove the converse implication 6) ⇒ 1) of Theorem 19. ▶ Theorem 34. Every step-size-controlled polynomial representation computes a primitive recursive function. Equivalently, Item 6) implies Item 1) in Theorem 19. Proof. Assume that f : Nd → Ne admits a step-size-controlled polynomial representation. Thus there exist a polynomial map P : RM → RM with rational coefficients, a primitive recursive precision threshold S : Nd → N, a primitive recursive observation count N : Nd × N → N, and designated output coordinates πe such that, for every x ∈ Nd and every s ≥ S(x), if z0 (x, s) := (x, 2−s , 0, . . . , 0) and zn+1 (x, s) := P(zn (x, s)), then ∥πe (zN(x,s) (x, s)) − f (x)∥∞ < 12 . Fix x. Since the coefficients of P are rational and the initial state z0 (x, s) is rational, every iterate zn (x, s) is rational. Under the standard encoding of rationals, polynomial evaluation with rational coefficients is primitive recursive. Hence the map (x, s, n) 7→ zn (x, s) is primitive recursive by primitive recursion on n. Since S and N are primitive recursive, so is the map x 7→ zN(x,S(x)) (x, S(x)). Composing with the output projection πe and coordinatewise rounding, we obtain a primitive recursive function of x. By the defining property of the representation, the rounded output is exactly f (x). Therefore f is primitive recursive. ◀ From recurrent ReLU computation back to primitive recursion. We next prove the converse implication 2) ⇒ 1). This is again elementary, because a fixed feedforward ReLU block with rational parameters induces a primitive recursive evolution on rational states. ▶ Theorem 35. Every exact recurrent ReLU computation computes a primitive recursive function. Equivalently, Item 2) implies Item 1) in Theorem 19. Proof. Assume that f : Nd → Ne admits an exact recurrent ReLU computation. Thus there exist a feedforward ReLU block R : Rm → Rm , a primitive recursive observation time T : Nd → N, and designated output coordinates πe such that, for every input x ∈ Nd , if z0 (x) := (x, 0, . . . , 0) and zn+1 (x) := R(zn (x)), then f (x) = πe (zT (x) (x)). A feedforward ReLU block is built from affine maps with rational parameters and coordinatewise applications of ReLU(t) = max(t, 0). On rational inputs, affine maps with rational parameters are primitive recursive, and ReLU is primitive recursive as well, since it is just the maximum of its input and 0. It follows that the one-step update (x, n) 7→ zn+1 (x) is primitive recursive on rational encodings. Hence, by primitive recursion on n, the full evolution (x, n) 7→ zn (x) is primitive recursive. Since T is primitive recursive and πe is a projection, x 7→ πe (zT (x) (x)) is primitive recursive. By assumption, this is exactly f (x). ◀
XX:21
XX:22
Primitive Recursion without Composition
11.2
From recurrent ρ-computation back to primitive recursion
We prove the converse implication (3) ⇒ (1) of Theorem 19. The argument is straightforward once one separates the two ingredients built into Definition 10. First, the recurrent dynamics is generated by iterating a fixed feedforward ρ-block on the compact cube [0, 1]m . Second, the admissibility assumptions on ρ imply that this block is primitive recursively approximable on rational inputs, with a primitive recursive modulus of continuity. Since the observation time and the required output precision are themselves primitive recursive, the whole orbit can be simulated primitive recursively up to the precision needed for the final ν-decoding. ▶ Lemma 36 (Effective evaluation of a fixed ρ-block). Let ρ be admissible, and let R : [0, 1]m → [0, 1]m be a fixed feedforward ρ-block. Then there exist primitive recursive functions EvalR : (Q ∩ [0, 1])m × N → Qm , µR : N → N, such that: 1. for every q ∈ (Q ∩ [0, 1])m and every s ∈ N, ∥EvalR (q, s) − R(q)∥∞ ≤ 2−s ; 2. for every u, v ∈ [0, 1]m and every s ∈ N, ∥u − v∥∞ ≤ 2−µR (s) =⇒ ∥R(u) − R(v)∥∞ ≤ 2−s . Proof. One proves both statements by induction on the depth of the fixed network R. Affine layers are exactly computable on rational inputs and have an obvious primitive recursive modulus. Coordinatewise applications of ρ are handled using the primitive recursive evaluator Evalρ and modulus µρ from Definition 9. Composing the finitely many layers yields the required evaluator EvalR and modulus µR . ◀ We also isolate the decoding step. ▶ Lemma 37 (Primitive recursive decoding of readable ν-codes). There exists a primitive recursive partial decoder Decodeν : Q∩[0, 1] → N with the following property: whenever q ∈ Q∩[0, 1] satisfies |q − ν(n)| < 14 ν(n) for some n ∈ N, then Decodeν (q) = n. Proof. By the remark following Definition 10, the integer n is uniquely determined by the displayed inequality. One may recover it by primitive recursive bounded search, for example as the least integer k such that 2k q ≥ 34 . ◀ ▶ Theorem 38. Every recurrent ρ-computation computes a primitive recursive function. Proof. Assume that f : Nd → Ne admits a recurrent ρ-computation. Fix the corresponding block R : [0, 1]m → [0, 1]m , observation time T : Nd → N, output-precision function S : Nd → N, and output coordinates πe from Definition 10. Let z0 (x) := (ν(x), 0, . . . , 0),
zn+1 (x) := R(zn (x)).
We first show that the state zT (x) (x) can be approximated primitive recursively to any prescribed precision. Let p : Nd → N be any primitive recursive function. Starting from the exact rational initial state z0 (x), define rational approximants e zn (x) recursively by e z0 (x) := z0 (x),
e zn+1 (x) := EvalR e zn (x), qn (x) ,
where the internal precisions qn (x) are chosen backwards using the modulus µR so that the final error at time T (x) is at most 2−p(x) . Since R is fixed and µR is primitive recursive, this backward precision schedule is itself primitive recursive. Hence there is a primitive recursive procedure which, given x, computes a rational vector e z(x) such that ∥e z(x) − zT (x) (x)∥∞ ≤ 2−p(x) .
Olivier Bournez
XX:23
We now choose p(x) large enough so that the corresponding approximation of the output coordinates still lies in the readable basin of ν( f (x)). By Definition 10, πe,i (zT (x) (x)) − ν( fi (x)) ≤ 2−S(x) ν( fi (x))
for every i.
Since the output is readable, we may choose p(x) primitive recursively large enough that the additional simulation error is smaller than the spare margin in this decoding condition. Then for every output coordinate i, 1 |πe,i (e z(x)) − ν( fi (x))| < ν( fi (x)). 4 Lemma 37 therefore yields fi (x) = Decodeν πe,i (e z(x)) . Thus each output coordinate of f (x) is obtained from x by a primitive recursive simulation followed by a primitive recursive decoder. Hence f itself is primitive recursive. ◀
11.3
Completion of the proof of the main theorem Theorem 19
We confirm here, if needed, that we can indeed assemble the implications established in the previous sections. Proof of Theorem 19. The implication (1) ⇒ (4) follows from Theorem 14 together with Theorem 15, which shows that the threshold-affine normal form can be taken with the uniform margin required by Item (4). The implication (4) ⇒ (2) is Theorem 22 (via Theorem 14). The implication (4) ⇒ (3) is Theorem 28. The implication (4) ⇒ (5) is Theorem 26. The implication (5) ⇒ (6) is Theorem 32. The converse implications back to Item (1) are: (2) ⇒ (1) is Theorem 24; (3) ⇒ (1) is Theorem 38; and (6) ⇒ (1) is Theorem 34. Together, these implications close the cycle, so Items (1)–(6) of Theorem 19 are all equivalent. ◀ ▶ Corollary 39 (Theorem 1, in the introduction). For a function f : Nd → Ne , the following are equivalent: 1. f is primitive recursive; 2. f admits a recurrent ReLU computation; 3. f admits a recurrent ρ-computation; 4. f admits a robust polynomial ODE representation; 5. f admits a step-size-controlled polynomial representation. Proof. This is Theorem 19 after eliminating the proof-level threshold-affine normal form of Item (4). ◀
12
Conclusion
We have characterized the primitive recursive functions as the class of functions computable by bounded iteration of a feedforward ReLU network, by bounded-range recurrent ρ-networks under ν-encoding, by robust polynomial ODEs, and by step-size controlled iteration of a polynomial map. The equivalences are mediated by a threshold-affine normal form that isolates the minimal branching structure behind a primitive recursive computation. The five formalisms realize the same class through genuinely different mechanisms, and the paper also clarifies a real asymmetry between them: exact recursive unfolding is natural in discrete time, while uniform rounding and autonomous phase organization are natural in continuous time.
XX:24
Primitive Recursion without Composition
The primitive recursive functions sit at the top of a rich hierarchy of classes obtained by restricting the recursion scheme: the Grzegorczyk hierarchy [27, 26], FP [22, 5], FPSPACE [29, 39], and more generally the subrecursive hierarchies surveyed in [21]. Each of our five characterizations offers a lens through which these subclasses may be redescribed, by restricting the corresponding dynamical resource (iteration count, domain diameter, precision, observation time, or degree of the polynomial map). Carrying out this programme systematically is a concrete direction for future work. Of particular interest in this respect is the step-size-controlled polynomial model. Although polynomial, it is genuinely distinct from the discrete-ODE frameworks of [9, 10, 1, 2, 3], which take the sign function and discrete derivatives as primitive operations: our model is strictly polynomial, with the step size, not necessarily unitary here, playing the role of an externally supplied precision parameter. This opens a parallel route to polynomial-only characterizations of the standard complexity classes. References 1
2
3
4
5 6
7
8
9
Melissa Antonelli, Arnaud Durand, and Juha Kontinen. A new characterization of FAC0 via discrete ordinary differential equations. In Rastislav Královic and Antonín Kucera, editors, 49th International Symposium on Mathematical Foundations of Computer Science, MFCS 2024, Bratislava, Slovakia, August 26-30, 2024, volume 306 of LIPIcs, pages 10:1–10:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2024. URL: https://doi.org/10.4230/LIPIcs.MFCS.2024.10, doi:10.4230/ LIPICS.MFCS.2024.10. Melissa Antonelli, Arnaud Durand, and Juha Kontinen. Characterizing small circuit classes from FAC0 to FAC1 via discrete ordinary differential equations. In Pawel Gawrychowski, Filip Mazowiecki, and Michal Skrzypczak, editors, 50th International Symposium on Mathematical Foundations of Computer Science, MFCS 2025, Warsaw, Poland, August 25-29, 2025, volume 345 of LIPIcs, pages 10:1–10:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2025. URL: https://doi.org/10.4230/LIPIcs.MFCS. 2025.10, doi:10.4230/LIPICS.MFCS.2025.10. Melissa Antonelli, Arnaud Durand, and Juha Kontinen. Towards new characterizations of small circuit classes via discrete ordinary differential equations. Theor. Comput. Sci., 1062:115655, 2026. URL: https://doi.org/10.1016/j.tcs.2025.115655, doi:10.1016/J.TCS.2025.115655. José L. Balcázar, Ricard Gavaldà, Hava T. Siegelmann, and Eduardo D. Sontag. Some structural complexity aspects of neural computation. In Proceedings of the Eigth Annual Structure in Complexity Theory Conference, San Diego, CA, USA, May 18-21, 1993, pages 253–265. IEEE Computer Society, 1993. doi:10.1109/SCT.1993.336521. Stephen Bellantoni and Stephen Cook. A new recursion-theoretic characterization of the poly-time functions. Computational Complexity, 2(2):97–110, 1992. doi:10.1007/bf01201998. Manon Blanc and Olivier Bournez. Simulation of turing machines with analytic discrete odes: Polynomialtime and space over the reals characterised with discrete ordinary differential equations. Journal of Logic and Analysis, 17, 2025. Accepted for publication. To appear. URL: http://logicandanalysis.org/ index.php/jla/article/view/505. Lenore Blum, Mike Shub, and Steve Smale. On a theory of computation and complexity over the real numbers; NP completeness, recursive functions and universal machines. Bulletin of the American Mathematical Society, 21(1):1–46, jul 1989. doi:10.1142/9789812792839\\\_0013. Olivier Bournez, Manuel Lameiras Campagnolo, Daniel Silva Graça, and Emmanuel Hainry. Polynomial differential equations compute all real computable functions on computable compact intervals. J. Complex., 23(3):317–335, 2007. doi:10.1016/j.jco.2006.12.005. Olivier Bournez and Arnaud Durand. Recursion schemes, discrete differential equations and characterization of polynomial time computations. In Peter Rossmanith, Pinar Heggernes, and Joost-Pieter Katoen, editors, 44th International Symposium on Mathematical Foundations of Computer Science, MFCS 2019, August 26-30, 2019, Aachen, Germany, volume 138 of LIPIcs, pages 23:1–23:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019. doi:10.4230/LIPIcs.MFCS.2019.23.
Olivier Bournez
10
Olivier Bournez and Arnaud Durand. A characterization of functions over the integers computable in polynomial time using discrete ordinary differential equations. Computational Complexity, 32(2):7, 2023. doi:10.1007/s00037-023-00240-1.
11
Olivier Bournez, Riccardo Gozzi, Daniel Silva Graça, and Amaury Pouly. A continuous characterization of PSPACE using polynomial ordinary differential equations. Journal of Complexity, 77:101755, august 2023. doi:10.1016/j.jco.2023.101755.
12
Olivier Bournez, Riccardo Gozzi, Daniel Silva Graça, and Amaury Pouly. A continuous characterization of PSPACE using polynomial ordinary differential equations. J. Complex., 77:101755, 2023. doi: 10.1016/j.jco.2023.101755.
13
Olivier Bournez, Daniel Silva Graça, and Amaury Pouly. Polynomial time corresponds to solutions of polynomial ordinary differential equations of polynomial length: The general purpose analog computer and computable analysis are two efficiently equivalent models of computations. In Ioannis Chatzigiannakis, Michael Mitzenmacher, Yuval Rabani, and Davide Sangiorgi, editors, 43rd International Colloquium on Automata, Languages, and Programming, ICALP 2016, Rome, Italy, July 11-15, 2016, volume 55 of LIPIcs, pages 109:1–109:15. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2016. ICALP’2016 Track B, Best Paper Award. doi:10.4230/LIPIcs.ICALP.2016.109.
14
Olivier Bournez, Daniel Silva Graça, and Amaury Pouly. Polynomial time corresponds to solutions of polynomial ordinary differential equations of polynomial length. Journal of the ACM, 64(6):38:1–38:76, 2017. doi:10.1145/3127496.
15
Olivier Bournez, Daniel Silva Graça, Amaury Pouly, and Ning Zhong. Computability and computational complexity of the evolution of nonlinear dynamical systems. In Paola Bonizzoni, Vasco Brattka, and Benedikt Löwe, editors, The Nature of Computation. Logic, Algorithms, Applications - 9th Conference on Computability in Europe, CiE 2013, Milan, Italy, July 1-5, 2013. Proceedings, volume 7921 of Lecture Notes in Computer Science, pages 12–21. Springer, 2013. doi:10.1007/978-3-642-39053-1\\_2.
16
Olivier Bournez and Emmanuel Hainry. Elementarily computable functions over the real numbers and r-sub-recursive functions. Theoretical Computer Science, 348(2-3):130–147, 2005. doi:10.1016/j. tcs.2005.09.010.
17
Olivier Bournez and Amaury Pouly. A survey on analog models of computation. In Vasco Brattka and Peter Hertling, editors, Handbook of Computability and Complexity in Analysis. Springer, 2021. doi:10.1007/978-3-030-59234-9\_6.
18
Michael S. Branicky. Universal computation and other capabilities of hybrid and continuous dynamical systems. Theoretical Computer Science, 138(1):67–100, 6~feb 1995. doi:10.1016/0304-3975(94) 00147-B.
19
Manuel L. Campagnolo. Computational complexity of real valued recursive functions and analog circuits. Thèse de doctorat, IST, Universidade Técnica de Lisboa, 2001.
20
Manuel L. Campagnolo, Cristopher Moore, and José Félix Costa. An analog characterization of the subrecursive functions. In P. Kornerup, editor, 4th Conference on Real Numbers and Computers, pages 91–109. Odense University Press, 2000.
21
Peter Clote and Evangelos Kranakis. Boolean Functions and Computation Models. Texts in Theoretical Computer Science. An EATCS Series. Springer, 2002. doi:10.1007/978-3-662-04943-3.
22
Alan Cobham. The intrinsic computational difficulty of functions. In Y. Bar-Hillel, editor, Proceedings of the International Conference on Logic, Methodology, and Philosophy of Science, pages 24–30. NorthHolland, Amsterdam, 1962.
23
Riccardo Gozzi. Analog Characterization of Complexity Classes. PhD thesis, Instituto Superior Técnico, Lisbon, Portugal and University of Algarve, Faro, Portugal, 2022.
24
Daniel Silva Graça, Manuel Lameiras Campagnolo, and Jorge Buescu. Robust simulations of turing machines with analytic maps and flows. In S. Barry Cooper, Benedikt Löwe, and Leen Torenvliet, editors, New Computational Paradigms, First Conference on Computability in Europe, CiE 2005, Amsterdam, The Netherlands, June 8-12, 2005, Proceedings, volume 3526 of Lecture Notes in Computer Science, pages 169–179. Springer-Verlag, 2005. URL: https://doi.org/10.1007/11494645\\_21, doi: 10.1007/11494645\\\_21.
XX:25
XX:26
Primitive Recursion without Composition
25 26 27 28 29
30
31
32 33 34
35 36 37 38 39 40
Daniel Silva Graça and José Félix Costa. Analog computers and recursive functions over the reals. J. Complex., 19(5):644–664, 2003. doi:10.1016/S0885-064X(03)00034-7. A. Grzegorczyk. Some classes of recursive functions. In Rozprawy Mate. IV. Warsaw, 1953. Andrzej Grzegorczyk. Computable functionals. Fundamenta Mathematicae, 42:168–202, 1955. doi: 10.4064/fm-42-1-168-202. Ker-I Ko. Complexity theory of real functions, volume 3 of Progress in theoretical computer science. Birkhäuser, Boston, 1991. doi:10.1007/978-1-4684-6802-1. Daniel Leivant and Jean-Yves Marion. Predicative functional recurrence and poly-space. In Michel Bidoit and Max Dauchet, editors, TAPSOFT’97: Theory and Practice of Software Development, 7th International Joint Conference CAAP/FASE, Lille, France, April 14-18, 1997, Proceedings, volume 1214 of Lecture Notes in Computer Science, pages 369–380. Springer, Apr 1997. doi:10.1007/BFb0030611. Wolfgang Maass. Bounds for the computational power and learning complexity of analog neural nets. In S. Rao Kosaraju, David S. Johnson, and Alok Aggarwal, editors, Proceedings of the Twenty-Fifth Annual ACM Symposium on Theory of Computing, May 16-18, 1993, San Diego, CA, USA, pages 335–344. ACM, 1993. doi:10.1145/167088.167193. Albert R. Meyer and Dennis M. Ritchie. The complexity of loop programs. In Proceedings of the 1967 22nd National Conference (ACM ’67), pages 465–469, New York, NY, USA, 1967. Association for Computing Machinery. doi:10.1145/800196.806014. Jerzy Mycka and José Félix Costa. A new conceptual framework for analog computation. Theor. Comput. Sci., 374(1-3):277–290, 2007. doi:10.1016/j.tcs.2007.01.005. P. Odifreddi. Classical Recursion Theory, volume 2. Elsevier, 1999. Amaury Pouly. Continuous models of computation: from computability to complexity. PhD thesis, Ecole Polytechnique and Unidersidade Do Algarve, Defended on July 6, 2015. 2015. https://pastel.archivesouvertes.fr/tel-01223284, Prix de Thèse de l’Ecole Polyechnique 2016, Ackermann Award 2017. Claude E. Shannon. Mathematical theory of the differential analyser. Journal of Mathematics and Physics MIT, 20:337–354, 1941. doi:10.1002/sapm1941201337. Hava T. Siegelmann. Neural networks and analog computation - beyond the Turing limit. Progress in theoretical computer science. Birkhäuser, 1999. Hava T. Siegelmann and Eduardo D. Sontag. Analog computation via neural networks. Theoretical Computer Science, 131(2):331–360, sep 1994. doi:10.1016/0304-3975(94)90178-3. Hava T. Siegelmann and Eduardo D. Sontag. On the computational power of neural nets. Journal of Computer and System Sciences, 50(1):132–150, feb 1995. doi:10.1006/jcss.1995.1013. David B. Thompson. Subrecursiveness: Machine-independent notions of computability in restricted time and storage. Mathematical Systems Theory, 6(1):3–15, 1972. doi:10.1007/BF01706069. Klaus Weihrauch. Computable Analysis - An Introduction. Texts in Theoretical Computer Science. An EATCS Series. Springer, 2000. doi:10.1007/978-3-642-56999-9.