REAL and Numerical Assurance — Chain Walkthrough

A first-principles walkthrough of the numerical assurance chain using a first-order lowpass filter as the running example.

All numbers in this document are produced by the tool itself and pinned by golden tests in tests/test_examples_t2_golden.py. If a number here ever disagrees with the test, the test is authoritative.

Supporting files:

examples/t2/lowpass.lola            — main example program
examples/t2/lowpass_report.py       — runs the full numerical assurance chain, prints all warrants
examples/t2/external_pid_contract.py — PID external real contract + two implementation warrants
examples/t2/pid_rust.rs             — illustrative Rust f64 implementation
examples/t2/pid_st.st               — illustrative ST LREAL implementation

---

Part I Why REAL means ℝ

LoLa's REAL type is the mathematical real numbers ℝ — not IEEE 754 floats.

This is not a naming convention. It is a proof-theoretic choice with consequences throughout the compiler:

• Z3 is the prover. When LoLa proves INVARIANT s >= -1.0, it calls Z3 on the mathematical statement ∀ s,u ∈ ℝ : s∈[-1,1] ∧ u∈[-1,1] ⟹ 0.98s + 0.02u ∈ [-1,1]. Z3 reasons in ℝ, not in Float64.

• Float64 is a representation, not a type definition. Choosing profile=FLOAT64_REPR tells the compiler "the hardware computes this state variable in 64-bit floating-point". This is a representation boundary, not a redefinition of what REAL means.

• Error is the gap between ℝ and Float64. The entire numerical assurance chain exists to bound this gap with a traceable, certified number.

The principle: a Float implementation can only approximate a mathematically-defined REAL program. It can never redefine it.

---

Part II A first REAL program

The running example is a first-order lowpass filter derived from the continuous ODE

ẋ = −2s + 2u

with pole at −2 rad/s. Forward Euler discretisation at T_s = 10 ms gives

s[k+1]  =  (1 + h·(−2))·s[k] + h·2·u[k]  =  0.98·s[k] + 0.02·u[k]

The LoLa source is in examples/t2/lowpass.lola:

FUNCTION_BLOCK LowpassFilter
(*
  Discrete first-order lowpass filter derived from
  the continuous ODE  ẋ = -2·s + 2·u  (pole at -2 rad/s)
  via Forward Euler with step h = 1/100 s:

      s[k+1] = (1 + h·(-2))·s[k] + h·2·u[k]
             = 0.98·s[k] + 0.02·u[k]

  T2-claim: for all k,  |s_k^machine  -  s_k^exact|  ≤  Δ
  where Δ is the CombinedImplementationTrajectoryBound computed
  from Float64 roundoff + physical sensor error + discretisation
  error + timing jitter.  All bounds are proved or declared;
  none are left implicit.
*)
VAR_INPUT
    u : REAL;       (* normalised input; sensor reading in [-1, 1] *)
END_VAR
VAR_OUTPUT
    s : REAL;       (* filter state / output *)
END_VAR
IMPLEMENTATION
    s: SET 0.98 * s + 0.02 * u OTHERWISE;
    REQUIRE u >= -1.0;
    REQUIRE u <= 1.0;
    INVARIANT s >= -1.0;
    INVARIANT s <= 1.0;
END_IMPLEMENTATION
END_FUNCTION_BLOCK

FUNCTION_BLOCK declares a stateful component with inputs, outputs, and an update rule.

s: SET expr OTHERWISE is a conditional update: s takes the value of expr at every step (the OTHERWISE clause covers the "no condition restricts this" case). In LoLa, SET expressions use REAL = ℝ arithmetic.

REQUIRE declares an input constraint. Z3 treats it as a precondition: "assume u ∈ [-1, 1] when proving the INVARIANT".

INVARIANT is a state invariant. Z3 proves it inductively: if s ∈ [-1, 1] and u ∈ [-1, 1], then 0.98·s + 0.02·u ∈ [-1, 1] in ℝ. The proof is exact and needs no margin — it works in ℝ.

Compiling this alone (no numerical assurance options) already runs the Z3 proof:

from lola.compiler import compile_source
from lola import profiles

result = compile_source(src, profile=profiles.FLOAT64_REPR)
iw = result.program.assurance.invariant_warrants
# iw["s"].status == "proved"

---

Part III From exact mathematics to Float64

Z3 proved the INVARIANT in ℝ. But the PLC runs in Float64. Do the Float64 computations stay close to the ℝ result?

profile=FLOAT64_REPR is the answer. It does two things:

1. Declares the representation. Every state variable s is mapped to a 64-bit IEEE 754 double. This is stored in the assurance report as a first-class annotation.

2. Triggers the Gappa roundoff proof. For each state variable, the compiler produces a Gappa script that asks: "given that the Float64 arithmetic computes fl(0.98)·fl(s) + fl(0.02)·fl(u) instead of the exact 0.98·s + 0.02·u, how large is the per-step error?"

The Lipschitz bound is computed algebraically:

L  =  |∂SET/∂s|  =  0.98  =  49/50

proved, because Z3 can differentiate the linear SET expression and Z3 arithmetic is exact. L < 1 means the discrete map is contractive: errors shrink geometrically rather than accumulate.

---

Part IV Roundoff and state accumulation

Gappa proves the per-step machine roundoff (StateTransitionBound):

ε_m  =  281 / 112589990684262400  ≈  2.50 × 10⁻¹⁵  per step

This is a verified, rigorous upper bound — not a heuristic estimate. Gappa constructs a formal certificate of the bound.

Error recurrence. If E_k is the accumulated error after k steps, then

E_{k+1}  ≤  L · E_k  +  ε_m

Starting at E_0 = 0 (at-rest initial condition), this recurrence gives:

E_k  ≤  ε_m · (1 + L + L² + … + L^{k-1})  =  ε_m · (1 − L^k) / (1 − L)

For L = 49/50 and k → ∞ (the asymptotic, worst-case bound):

E_∞  =  ε_m / (1 − L)  =  ε_m / (1/50)  =  50 · ε_m

The TubeClosure + TrajectoryBound gives the proved value:

E_∞ (machine only)  =  345 / 2251799813685248  ≈  1.53 × 10⁻¹³

(Exact values are pinned to the toolchain's prover, Gappa 1.8.0: the adaptive search certifies the first bound the prover accepts, so a different Gappa version yields a slightly different — equally sound — number.) The machine-only bound sits around 10⁻¹³ — far tighter than anything a physical sensor can provide. Parts V–VII add the real engineering uncertainties.

---

Part V Physical uncertainty

The sensor reading of u is not perfect. The physical measurement error is declared as:

physical_errors={"u": Fraction(1, 100)}   # |δu| ≤ 0.01  (1 % of range)

LoLa creates a PhysicalInputBound: a declared, assumed warrant that the physical sensor error on u is at most 1/100.

The InputSensitivityBound is then proved algebraically:

|∂SET/∂u|  =  0.02  =  1/50     [proved]

so a 1/100 error on u contributes at most

ε_p  =  1/100 · 1/50  =  1/5000  =  0.0002  per step

The StepErrorBudget adds machine roundoff and physical input error:

machine_roundoff      ≈  1.17 × 10⁻¹⁵   [proved — inside the wider tube]
physical_input_error  =  1/5000 = 0.0002 [proved-under-assumptions]
───────────────────────────────────────────
SEB total             ≈  0.0002

Physical uncertainty dominates machine roundoff by seven orders of magnitude. The machine is not the problem; the sensor is.

---

Part VI Sampling, discretization, and jitter

The model s[k+1] = 0.98·s[k] + 0.02·u[k] is a discrete approximation of the continuous ODE ẋ = −2s + 2u. LoLa formalises this relationship with a SamplingWarrant:

sampling_models={"s": SamplingSpec(
    sample_period=Fraction(1, 100),          # T_s = 10 ms
    method="forward-euler",
    continuous_rhs=RealBin("+",
        RealBin("*", RealConst(Fraction(-2)), Ref("s")),
        RealBin("*", RealConst(Fraction(2)),  Ref("u")),
    ),
)}

This warrant says: "the discrete LoLa program is the Forward Euler approximation of this specific continuous ODE". It is declared, not proved — LoLa does not prove that the ODE models the physical plant. That is the engineer's responsibility.

Discretisation error. Forward Euler has a local truncation error of h²/2 · |ẍ|_max. Given the input derivative bound |u̇| ≤ 0 (piecewise constant, declared):

ε_trunc  =  121 / 244500  ≈  4.95 × 10⁻⁴   [proved-under-assumptions]

Timing jitter. Real PLCs do not sample at exactly T_s = 10 ms. Declare a jitter bound J = 1 ms:

jitter_bounds={"s": Fraction(1, 1000)}

The jitter step error bounds the effect of timing variation on one step:

ε_jitter  =  2 / 489  ≈  4.09 × 10⁻³   [proved-under-assumptions]

The effective single-step timing+discretisation error is:

ε_eff  =  ε_trunc + ε_jitter  =  1121 / 244500  ≈  4.58 × 10⁻³

The worst-case discrete Lipschitz at the extremes of the jitter interval h ∈ [9 ms, 11 ms] is:

L_Φ^max  =  max(|1 − 2·0.009|, |1 − 2·0.011|)  =  491 / 500  =  0.982

---

Part VII The unified engineering error report

With physical errors, sampling, and jitter all declared, the compiler assembles the CombinedImplementationTrajectoryBound:

ε_combined  =  SEB.total + ε_eff  ≈  0.004785   [proved-under-assumptions]
L_eff       =  max(L_machine, L_Φ^max)  =  491/500 = 0.982

Because L_eff < 1, the combined system is still contractive. The asymptotic trajectory bound is:

E_∞  =  ε_combined / (1 − L_eff)  ≈  0.2658
E_100 ≈  0.2226   (after 100 steps from rest)

Interpretation: starting from an exact initial condition, the Float64 PLC implementation running with 1 % sensor noise and 1 ms jitter will have an accumulated state error of at most 0.266 — in normalised units, about 26.6 % of the full scale. This is the price of a 1 ms jitter with a 10 ms sample period (10 % relative jitter); reducing J tightens the bound sharply.

The EngineeringErrorReport assembles all of this into one machine-readable artefact with a stable digest. Running PYTHONPATH=. python examples/t2/lowpass_report.py prints:

══════════════════════════════════════════════════════════════════════════════
  Engineering Error Budget — s
══════════════════════════════════════════════════════════════════════════════
  Representation : float64
  Overall status : proved-under-assumptions
  Warrant ceiling: proved-under-assumptions

┌─ Semantics ──────────────────────────────────────────────────────────────────
  LoLa SET (discrete) : ((0.98 * s) + (0.02 * u))
  Continuous ODE rhs  : Add(Mul(Integer(-1), Integer(2), Symbol('s', real=True)),
                            Mul(Integer(2), Symbol('u', real=True)))

┌─ Timing ─────────────────────────────────────────────────────────────────────
  Sample period T_s   : 0.010000  (1/100) s
  Jitter bound J      : 0.001000  (1/1000) s
  h ∈ [0.009000  (9/1000), 0.011000  (11/1000)] s

┌─ Per-step error budget ──────────────────────────────────────────────────────
  Machine roundoff ε_m           0.000000  (3297/2814749767106560000)  [proved]
  Physical input  ε_p            0.000200  (1/5000)  [proved-under-assumptions]
  Discretization  ε_d            0.000495  (121/244500)  [proved-under-assumptions]
  Jitter          ε_j            0.004090  (2/489)  [proved-under-assumptions]
  ────────────────────────────────────────────────────────────────────────────
  Combined ε_combined                      0.004785

┌─ Trajectory bound ───────────────────────────────────────────────────────────
  L_machine   : 0.980000  (49/50)
  L_Φ^max     : 0.982000  (491/500)
  L_eff       : 0.982000  (491/500)  (max of above)
  Regime      : contractive
  Horizon N   : 100
  E_N         : 0.222600
  E_∞         : 0.265826

┌─ Declared assumptions ───────────────────────────────────────────────────────
  [SamplingWarrant]  Nominal sample period T_s = 1/100 s
  [InputDerivativeBound]  |u_dot| <= 0 (declared)
  [InputRangeWarrant]  u in [-1, 1] (from REQUIRE)
  [JitterBound]  Timing jitter |T_k - T_s| <= J = 1/1000 s
  [PhysicalInputBound]  |delta_u| <= 1/100 (declared)

┌─ Proof chain ────────────────────────────────────────────────────────────────
  [CombinedStepBudget]  proved  ceil=proved-under-assumptions
  [SemanticCompositionWarrant]  proved  ceil=proved-under-assumptions
  [LipschitzBound]  proved
  [CombinedImplementationTrajectoryBound]  proved  ceil=proved-under-assumptions
  [JitterStepBound]  proved  ceil=proved-under-assumptions
  [StepErrorBudget]  proved  ceil=proved-under-assumptions

┌─ Report digest ──────────────────────────────────────────────────────────────
  77d8ebfba73e73748ae444c8ff35b17836aa48211d52fc053296e05d7b03a265
══════════════════════════════════════════════════════════════════════════════

The digest 77d8ebf... is a SHA-256 over the complete report payload. It changes if any input (source, SamplingSpec, jitter bound, sensor error, or horizon) changes. This makes the report a traceable artefact: a CI log entry or a datasheet footnote can reference this digest to pin the exact claim.

---

Part VIII External implementations

Everything in Parts I–VII concerns REAL programs that LoLa's compiler can directly analyse — the FUNCTION_BLOCK source is present and the SET expressions are known.

Some components cannot be expressed in LoLa: a hardware codec, a vendor-supplied PID library, a DSP filter running in firmware. For these, LoLa provides ExternalRealContract and ExternalImplementationWarrant.

The principle

A Float implementation can never define the mathematical meaning of a LoLa REAL component. It can only approximate a mathematically-defined ExternalRealContract.

An ExternalRealContract holds the exact REAL (= ℝ) semantics of a component as LoLa IR expressions — the same RealBin/RealConst/Ref nodes used inside FUNCTION_BLOCKs. The contract is always status="assumed": the engineer declares the intended mathematical behaviour. LoLa does not prove that the hardware or library actually computes it.

An ExternalImplementationWarrant binds one concrete implementation (e.g. a specific Rust crate at a specific version) to an ExternalRealContract. It declares the per-output and per-state error bounds relative to the contract's exact REAL values, and the representation (float64 or lreal64). It is status="assumed" by default — the error bounds are declared, not proved by Gappa (though a future extension could add Gappa or Kani proofs).

The proof DAG is:

ExternalRealContract[assumed]
        ↑  (assumption_dep)
ExternalImplementationWarrant[assumed]

The ERC is a leaf: it has no LoLa source to verify against. The EIW's assumption_dep_warrant_keys records the ERC's digest, so the ProofRegistry can enforce that the warrant references a specific, sealed contract.

---

Part IX PID preview: one mathematical contract, two implementations

The PID controller is a natural example for ExternalRealContract because:

• It is stateful (integral accumulator, previous error).
• Its mathematical semantics are well-defined and compact.
• Different implementations (Rust f64 on an embedded target vs ST LREAL on an OpenPLC runtime) can both be certified against the same sealed contract.

The contract

The controller computes (with Kp = 2, Ki = 1/10, Kd = 1/100):

e[k]            =  setpoint − measurement
integral[k+1]   =  integral[k] + (Ki · dt) · e[k]
control[k]      =  Kp·e[k] + integral[k] + Kd·(e[k] − e_prev[k]) / dt
e_prev[k+1]     =  e[k]

In examples/t2/external_pid_contract.py, these are written as exact REAL IR expressions and sealed into an ExternalRealContract:

from fractions import Fraction
from lola.ir import RealBin, RealConst, Ref
from lola.contract_ranges import Range
from lola.external_real_contract import make_external_real_contract

_e = RealBin("-", Ref("setpoint"), Ref("measurement"))
_i_next = RealBin("+", Ref("integral"),
    RealBin("*", RealBin("*", RealConst(Fraction(1,10)), Ref("dt")), _e))
_deriv = RealBin("*", RealConst(Fraction(1,100)),
    RealBin("/", RealBin("-", _e, Ref("e_prev")), Ref("dt")))
_control = RealBin("+",
    RealBin("+", RealBin("*", RealConst(Fraction(2)), _e), Ref("integral")),
    _deriv)

contract = make_external_real_contract(
    component_name="PidController",
    kind="stateful-block",
    input_names=["setpoint", "measurement", "dt"],
    output_names=["control"],
    state_names=["integral", "e_prev"],
    exact_outputs={"control": _control},
    exact_next_state={"integral": _i_next, "e_prev": _e},
    input_ranges={
        "setpoint":    Range(Fraction(-10), Fraction(10)),
        "measurement": Range(Fraction(-10), Fraction(10)),
        "dt":          Range(Fraction(1, 1000), Fraction(1, 10)),
    },
    state_invariants={
        "integral": Range(Fraction(-100), Fraction(100)),
        "e_prev":   Range(Fraction(-20), Fraction(20)),
    },
)
# contract.digest == "a16f5dcb5d7fcb970d7ca8b12e2601e628936f2edb9b2e18e930999fdfffd8bc"

The digest a16f5dcb... identifies this exact mathematical specification. Any change to Kp, Ki, Kd, or the formula changes the digest.

Two implementations, one contract

examples/t2/pid_rust.rs implements the controller in Rust using f64. The ExternalImplementationWarrant declares:

w_rust = make_external_implementation_warrant(
    contract=contract,
    backend="rust",
    implementation_id="pid_rust_f64_v1",
    representation="float64",
    output_error_bounds={"control": Fraction(1, 10000)},   # ε ≤ 0.0001
    state_transition_error_bounds={
        "integral": Fraction(1, 100000),
        "e_prev": Fraction(0),
    },
    state_lipschitz_bounds={"integral": Fraction(1, 5), "e_prev": Fraction(1)},
    status="assumed",
)
# w_rust.digest == "7f6d5374babd4db989cdbde1ba058464902907e743eecbb883edc1e3a6b44e9a"

examples/t2/pid_st.st implements the same controller as an IEC 61131-3 FUNCTION_BLOCK using LREAL. A separate warrant declares looser bounds:

w_st = make_external_implementation_warrant(
    contract=contract,
    backend="structured-text",
    implementation_id="pid_st_lreal64_v1",
    representation="lreal64",
    output_error_bounds={"control": Fraction(1, 1000)},    # ε ≤ 0.001
    state_transition_error_bounds={
        "integral": Fraction(1, 10000),
        "e_prev": Fraction(0),
    },
    state_lipschitz_bounds={"integral": Fraction(1, 5), "e_prev": Fraction(1)},
    status="assumed",
)
# w_st.digest == "391eae6f573652eeb4f1130339db3c0811c0ea099615729cace1a128cdb1ba66"

Both warrants reference the same sealed contract (mathematical_contract_digest == "a16f5dcb..."), but with different implementation-specific bounds:

Backend	Representation	ε_out[control]
Rust f64	float64	1/10000 = 0.0001
ST LREAL	lreal64	1/1000 = 0.001

The Rust bound is 10× tighter because the Rust implementation avoids intermediate LREAL rounding modes and accumulates fewer inter-operation conversions.

What this buys

The contract digest is proof that both implementations commit to the same mathematical PID definition. A system integrator can:

1. Choose the implementation that meets their error budget.
2. Verify that the bound declared in the EIW is consistent with their requirements.
3. Record the contract digest in the system's assurance case to pin exactly which mathematical PID definition was certified.

The error bounds in the EIW are declared, not Gappa-proved. A future status="proved" path (using Kani or Gappa with Rust or LREAL backends) would replace the assumption with a verified bound — without changing the mathematical contract.

---

Part X What LoLa proves — and what it does not

Proved (Z3, Gappa, or algebraic chain)

Claim	Method	Status
INVARIANT s ∈ [−1, 1] holds at every step	Z3 (ℝ)	proved
Lipschitz constant L = 49/50	algebraic	proved
Per-step machine roundoff ε_m ≈ 5.06 × 10⁻¹⁶	Gappa	proved
Input sensitivity	∂SET/∂u	= 1/50
Discretisation correspondence: SET = s + h·f(s,u) exactly	algebraic	proved
Discretisation step error ε_trunc = 121/244500	algebraic	proved-under-assumptions
Jitter step error ε_jitter = 2/489	algebraic	proved-under-assumptions
Combined trajectory bound E_∞ ≈ 0.2658	algebraic chain	proved-under-assumptions

Assumed (declared by the engineer)

Claim	Warrant	Impact
Sensor error	δu	≤ 1/100
Sample period T_s = 1/100 s	SamplingWarrant	drives ε_trunc and L_Φ^max
Continuous ODE is ẋ = −2s + 2u	SamplingWarrant	defines the reference
Input is piecewise constant (	u̇	= 0)
Timing jitter	T_k − T_s	≤ 1/1000 s
PID ExternalRealContract semantics	ExternalRealContract	pins the contract
PID implementation error bounds	ExternalImplementationWarrant	per-implementation claim

Not in scope (by design)

Functional correctness of external implementations. The ExternalImplementationWarrant says "our Rust PID outputs are within 0.0001 of the contract's ℝ value". It does not say "the Rust code has no bugs". A Kani proof of the Rust code would upgrade the warrant to status="proved".

Physical plant stability. LoLa proved that the discrete filter is contractive (L < 1). It did not prove that the continuous closed-loop system is stable. Continuous stability analysis is a separate activity (e.g. Lyapunov, Nyquist) that takes the LoLa DiscreteMapLipschitz as one input.

Scheduling and OS timing. The JitterBound is declared. LoLa does not prove that the scheduler guarantees |T_k − T_s| ≤ J. That is a realtime systems concern, properly addressed by WCET analysis or measurement.

Numerical conditioning of sequences. LoLa bounds the error after N steps from a single exact initial condition. It does not address ensemble behaviour, noise spectral density, or aliasing — those are signal-processing concerns that use the LoLa bounds as inputs, not outputs.

---

End of LoLa by Example.

Next: docs/architecture.md gives the full proof-registry design. docs/adr-numerical-semantics.md records the decision to use ℝ rather than a Float semantics as the primary type.