docs(adr): ADR 009 leads with incremental formulas; exact f64 sums; ordered aggregates (improvise-felo.4)

- Reframed: incremental formula evaluation for streamed data (ADR 003) is the
  primary benefit, and ADR 008's grid is the second consumer.
- New clause 5a: aggregate over an expression, as a lifted map feeding the
  group aggregate, still evaluated by eval_expr.
- Clause 4: ordered state for LAST/FIRST (felo.7) and semi-additive measures
  (felo.8). New question 11.
- Clause 8: DBSP alone gives no floating-point benefit. An exact superaccumulator
  makes f64 SUM a true group. It is measured in retraction_drift: bit-equal after
  1M cycles, about 3.2 ns per update, 568 bytes per group. It is correct against
  math.fsum on 20,000 cases, and xsum 0.1.6 is rejected. The principle of
  deferring lossy operations, and a table of what each exact domain covers.
  Question 3 gains option (c) and recommends it.

No clause or question number moved.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdfEGXA7aJ6ftMsCc9PJze
This commit is contained in:
Claude (agent-felo4)
2026-09-24 21:53:29 -07:00
committed by Edward Langley
parent 0ab3ac7dbf
commit a2836f44f5
2 changed files with 544 additions and 70 deletions
+286 -44
View File
@@ -42,8 +42,38 @@ So the choice between Feldera's crate and our own is already made (see Alternati
what was checked). This ADR says what "the core of dbsp" is for improvise, where it
plugs in, and what it may not change.
## What it is for
**The primary benefit is incremental formula evaluation, so that data can stream.**
That is the user's own framing, both in the prompt that led to this ADR and in their
later question: "Isn't the benefit of adopting DBSP primarily incremental formula
updates?" Yes. When a cell changes, whether from a keystroke, a client's frame or a
live source (ADR 003, clause 10), only the aggregates that cell feeds are adjusted,
and only the formula cells whose inputs moved are re-evaluated. Nothing is recomputed
from scratch. The machinery is:
- the aggregates formulas read are **maintained**, not re-folded (clause 4);
- an aggregate over a per-record expression such as `SUM(Revenue * Price)` is a
maintained operator too (clause 5a);
- the circuit decides **which** formula cells moved, and the one existing evaluator
recomputes only those, in dependency order (clauses 5 and 6).
**The grid is the second consumer.** ADR 008 needs the same maintained state for other
purposes: occupancy for `prune_empty`, the occupancy delta for its frame observation,
and column totals over rows never scrolled to (clause 9). It gets them from the same
circuit, not from a second mechanism.
## Context
### The primary case: data that streams
A live-backed model (ADR 003) is a source whose rows change underneath the document.
Today the only way to reflect a change is the path every edit takes: recompute every
formula cell of the document (next subsection). That is tolerable when a human types
one cell at a time. It is not tolerable when a source delivers a stream of changes,
each touching a few records of a large table. Incremental evaluation makes the cost of
a change proportional to what the change touches, not to the document.
### Evaluation today is a from-scratch fold, rerun after every batch
- **Every batch recomputes every formula cell.** `App::apply_effects`
@@ -98,10 +128,14 @@ incremental evaluator could quietly disagree with the eager one.
## Decision
**A new, problem-agnostic workspace crate implements the DBSP core. improvise-core
defines a trait for the derived state it reads. Today's evaluator is one implementation
of that trait, and a circuit built on the new crate is the other. Both stay in the tree,
and a differential test holds them equal.** Eleven clauses.
**Formulas are evaluated incrementally. A change re-evaluates only the formula cells
whose inputs moved, over aggregates that are maintained rather than re-folded. The
machinery is a new, problem-agnostic workspace crate implementing the DBSP core.
improvise-core defines a trait for the derived state its evaluator reads. Today's
evaluator is one implementation of that trait, and a circuit built on the new crate is
the other. Both stay in the tree, and a differential test holds them equal.** Eleven
clauses, plus clause 5a. It is numbered 5a so that no clause number or question number
cited elsewhere moves.
### 1. The crate: `zset` (name is question 1)
@@ -265,7 +299,7 @@ comes from its `_Measure` coordinate through `measure_agg`, as `agg_for` does
| answer | state per group | cost per delta element | exact on f64? |
|---|---|---|---|
| `Count` | count of values, count of numbers | O(1) | yes (integers) |
| `Sum` | running sum | O(1) | **no** (clause 8) |
| `Sum` | running sum (an exact accumulator under question 3's recommendation) | O(1) | **no** as an f64; **yes** as the accumulator (clause 8) |
| `Avg` | derived: `Sum / Count` | O(1) | no, through `Sum` |
| `Min`, `Max` | multiset of the numbers | O(log d) | yes (order-free) |
| unanimity and `mixed` | multiset of all values | O(log d) | yes (order-free) |
@@ -294,6 +328,43 @@ one value, and `f64::min` skips a NaN. So the aggregator's *read* reproduces eag
rules explicitly and does not inherit `total_cmp`'s. The differential test generates
signed zeros and non-finite values (core.md §16) to pin it.
#### Ordered aggregates: `LAST` and `FIRST` need state ordered by an item, not a value
improvise-felo.7 asks for `LAST` and `FIRST` over a fold: the closing or opening value
of a stock price, say, across the months of a year. improvise-felo.8 asks for
semi-additive measures: `sum` across one category and `last` across another. Both are
open tickets with their own semantics to decide. What this ADR must say is whether its
per-group state can carry them. Today it cannot: the integrated state above is an
**unordered** multiset, keyed by the value itself.
**The structure.** An ordered aggregate keeps, per group, a `BTreeMap` keyed by the
**ordering category's item order**, with the tie-break felo.7 decides. Each entry holds
the values at that position with their weights. `LAST` reads the entry at the greatest
key, and `FIRST` the least.
- **It is the same shape as `Min`/`Max`.** Those also keep an integrated ordered map and
read an extreme; they order by the value, where `LAST`/`FIRST` order by the item's
position and read the value stored there. So one generic operator covers all four:
an integrated map keyed by `K`, read at an extreme. `Min`/`Max` instantiate it with
`K = value`, and `LAST`/`FIRST` with `K = (item position, tie-break)`.
- **Cost:** O(log d) per delta element, where d is the number of distinct positions in
the group, and one entry per contributing record. That is the same order as
`Min`/`Max`'s multisets.
- **What moves the order.** Adding an item in the middle of the ordering category, or
reordering its items, changes every stored position. So the ordering category's
item order must be part of the operator shape (clause 7), and a change to it rebuilds
the ordered aggregators that use it. The alternative is to key by a stable item
identity (ADR 007) and look the order up at read time. That makes a reorder free, but
a read then costs O(d). This is question 11.
- **Semi-additive measures (felo.8) are aggregates of aggregates.** For example,
`last` across `Date` and then `sum` across `Account`. The inner aggregate emits a
delta (retract the old value, insert the new), and the outer consumes it. So the
composition is incremental at O(log d) per stage. Which stage is inner is felo.8's
semantic decision, because sum-of-last and last-of-sum differ. The circuit takes
either order.
- **Exactness.** `LAST`/`FIRST` return a stored value unchanged, so they are exact on
f64 like `Min`/`Max`. A semi-additive `sum` stage still needs question 3's answer.
### 5. Formulas: the circuit finds what moved, and the evaluator recomputes it
**Stems are a `distinct`.** Today's stem set (`types.rs:986-997`) is every stored key
@@ -340,14 +411,61 @@ With the trait, both reads go to the maintained state.
not match falls back to raw aggregation (`types.rs:1047-1055`), so its input is its own
group and nothing else.
**An aggregate over an expression replaces `Expansion::members`.** `SUM(expr)` pools
`expr` evaluated at each member record (`agg_values`, `types.rs:89-139`). Incrementally,
that is `map(record ↦ expr at record)`, then `index_by(base key)`, then `Sum`. The
member records never have to be listed. **The boundary:** the map is sound only when
`expr` reads the record's own cells. An inner expression that reads another record, or
holds a nested aggregate, stays on the eager path for that formula, and the
differential test pins the boundary. Until `Sum` is exact, this path carries clause 8's
caveat. ADR 008 clause 5 adopts the same boundary.
**An aggregate over an expression** is the one place a formula needs an expression
*inside* an operator. It has its own clause, 5a.
### 5a. Aggregate over an expression: a lifted map feeding the group aggregate
`SUM(Revenue * Price)` does not sum a stored measure. It computes a value **per
record**, `Revenue * Price` at that record, and pools those values across the records
of a group. Today `agg_values` (`types.rs:89-139`) lists the member records of each
context and calls the record-level evaluator on each (`member_records`, then
`eval_at`). At a fold or a column total, the contexts come from `Expansion::members`
(`types.rs:201-213`), which materializes the cross product of the spanned categories'
items. That is the one O(grid) path in totals (ADR 008 clause 5).
**The design.** It is a three-stage pipeline over the record stream, and each stage is
an operator from clause 1:
1. **A lifted map over the per-record expression.** For each record in the cell
delta's `record_of` set, the map produces `(base key, expr at record)`. A record
whose cells changed is retracted at its old value and inserted at its new one, so
the map's output is itself a delta. `map` is linear, so this is its own incremental
form.
2. **`index_by(base key)`.** The base key is the record's key with the target category
removed, which is what `eval_agg` pools over (`types.rs:55-87`). An inner `WHERE`
(`SUM(expr WHERE Cat = Item)`) is a `filter` before the index. It is the same
rewrite `eval_agg` applies to its contexts (`types.rs:66-72`).
3. **The group aggregate** from clause 4: `Count`, `Sum`, `Min` or `Max` over the
mapped values, with `Avg` derived.
Member records never have to be listed, and no cross product is ever built. A context
in today's cross product that holds no records contributed nothing to the pool, so
pooling only the records that exist is the same answer.
**This is not compiling formulas into operators.** The operator is the map; the
expression inside it is still evaluated by `eval_expr` at the record, exactly as
`agg_values` calls it today. The circuit decides *which* records to re-evaluate (those
the delta names), and the one evaluator computes each value. Clause 2's rule stands:
there is one evaluator.
**The boundary.** The map is sound only when `expr` reads the record's own cells, so
that its value at a record changes only when that record changes. An inner expression
that reads another record, holds a nested aggregate, or refers to a formula whose value
depends on other records, is not a function of one record. That formula stays on the
eager path, and the differential test pins the boundary. ADR 008 clause 5 adopts the
same boundary.
**What stays eager until step 5, and until step 6.**
- Before step 5 (formulas), every aggregate over an expression is evaluated as today,
through `Expansion::members`.
- From step 5, the map and the index are maintained, so the cross product is gone.
`Count`, `Min` and `Max` over the mapped values are read from the maintained state,
since they are order-free and bit-exact on f64 (clause 8).
- `Sum` and `Avg` over the mapped values are **re-folded** from the group's maintained
values, in the store's record order, until step 6. That keeps them bit-equal to
eager, at a cost of O(group) per change rather than O(grid).
- From step 6 (on question 3's answer), `Sum` is maintained like any other.
### 6. Circular references: DBSP's recursion is not our semantics
@@ -425,7 +543,7 @@ stored cells, however many pages there are. Filter coordinates still filter; onl
coordinates key the slice. If the user answers the other way, a page change is an
operator-shape change like any other.
### 8. Floating point: retraction drifts, and exactness is ADR 005's
### 8. Floating point: retraction drifts, and exactness comes from the representation
f64 addition is not associative, so a sum maintained by `+v` and `−v` is not the sum
the eager fold computes. **Measured**, not asserted:
@@ -453,12 +571,12 @@ A Neumaier-compensated running sum recovers that case exactly. In the random cas
stays about as close to the exact sum as eager does, **but it is never bit-equal to
eager**, because eager folds in store order and is itself inexact.
Four conclusions follow:
Four conclusions follow about a **floating-point** running sum:
- **Bit-equality between an incremental and an eager f64 `Sum` is not attainable** by
any running-sum technique. The differential test can compare `Count`, `Min`, `Max`,
unanimity and occupancy bit for bit today, because they are order-free. It can compare
`Sum` and `Avg` only with a tolerance.
any f64 running-sum technique, naive or compensated. So the differential test can
compare `Count`, `Min`, `Max`, unanimity and occupancy bit for bit, because they are
order-free, but it can compare an f64 `Sum` and `Avg` only with a tolerance.
- **The gap reaches the epsilon.** At 18k cells, the incremental sum and the eager sum
differ by at least `FLOAT_EQ_EPSILON` at 13 of the 16 checkpoints. So
`IF(SUM(..) = x, ..)` can answer differently on the two paths. That is a behaviour
@@ -466,20 +584,115 @@ Four conclusions follow:
- **This is already a bug in eager, which the measurement found in passing.** With 3
decades of magnitude, the eager 18k-cell sum is up to 5.8e-10 from the exact decimal
sum, and above 1e-10 at 7 of its 8 checkpoints. So `=` against a decimal a user typed
can already be wrong today. That supports ADR 005 on its own
merits, and it is recorded here for improvise-l2b8.
- **Exact decimals dissolve both problems.** Under ADR 005, `Sum` is a true
`AbelianGroup`, retraction is exact, the order of folding stops mattering (gotcha 18's
whole class), and the differential test is bit-exact everywhere.
can already be wrong today. That is recorded here for improvise-l2b8.
- **Only an exact representation dissolves the problem.** The next two subsections give
two: an exact accumulator for f64 sums, and ADR 005's decimals.
**So ADR 005 is a prerequisite for the maintained `Sum`/`Avg` on the displayed path, but
not for the rest of this ADR** (question 3). Occupancy, `Count`, `Min`, `Max`,
unanimity, stems and invalidation can all land before it. Until ADR 005 lands, a `Sum`
or `Avg` the grid displays, including a column total that holds one, is read from the
eager fold. The maintained sum runs only inside the differential test, under a stated
tolerance. If the user wants maintained sums before ADR 005, the fallback is a
compensated sum that re-folds a group from its values when cancellation is detected.
The measurement above is the evidence for and against that fallback.
#### Does DBSP itself help with floating point? No.
The user asked whether rebuilding the evaluator on DBSP brings any correctness benefit
for floating point. On its own, **none**:
- **Only the Z-set weights are exact.** They are integers (clause 1). The *values* being
summed are whatever type the aggregator holds.
- **An f64 `Sum` maintained by retraction is worse than eager, not better.** The table
above is the evidence: it drifts further from the exact sum the longer it runs, and
the absorption case loses the answer entirely.
**The benefits come from two other places.**
- **The representation.** An exact accumulator (next subsection) or exact decimals
(ADR 005) make `Sum` a true abelian group. Then retraction is exact, incremental
equals eager bit for bit, summation order stops mattering (gotcha 18's order-dependent
`SUM` goes away as a class), and today's eager error beyond `FLOAT_EQ_EPSILON` goes
with it.
- **The differential test** between `Eager` and `Incremental` (clause 2) is itself a
correctness gain. Two independent evaluations of every formula cell catch evaluator
bugs in either one, which nothing in the tree does today.
What DBSP does contribute is **placement**. Every aggregate is a named operator with a
chosen algebra, so the circuit states exactly where values enter an exact domain and
where they leave it. That is the next subsection's principle.
#### An exact accumulator for f64 sums, and deferring the lossy operations
**A superaccumulator sums f64 values exactly.** It is a fixed-point integer wide enough
for the whole f64 range: in units of 2⁻¹⁰⁷⁴, a finite f64 is a 53-bit integer shifted
by at most 2045 bits. This is the idea behind Kulisch accumulators and Radford Neal's
xsum, and Shewchuk's expansions (behind Python's `math.fsum`) are a related technique.
Adding or retracting a value is integer addition, so the accumulator is a true abelian
group over the finite f64s. **Rounding happens once, when the value is read.** Then:
- `Sum` is the exact sum, rounded once. It is independent of order.
- `Avg` is the exact sum divided by the count, rounded once. It is correctly rounded
if the division is done on the exact value; dividing the rounded sum rounds twice,
and step 6 must say which it does.
- `Min` and `Max` are already exact.
- Non-finite values are counted (NaN, +∞, −∞), so they can be retracted too.
**Measured**, in the same example, by `ExactSum` (147 lines, including its rounding):
- **Maintained equals fresh, bit for bit.** At all 25 checkpoints above, an `ExactSum`
maintained through a million retract-and-insert cycles equals a fresh `ExactSum` over
the current values. That includes the absorption case, which reads exactly 10.00.
Against the exact *decimal* sum, it differs at 1 checkpoint of 25, by 9.1e-13. That
gap is the rounding of the typed decimals into f64, not the sum.
- **It is correct.** It was checked (not committed) against Python's `math.fsum` on
20,000 random cases. The cases include subnormals, values near the overflow bound,
and an insert-then-retract of an equal-sized noise set. There were 0 mismatches; the
1,596 cases where `fsum` overflows were checked against exact `Fraction` sums.
- **Cost, Apple M4 Max, best of 3 over 10M values:**
- about 3.2 ns per insert or retract, against about 2.3 ns for a plain f64 add in the
same loop;
- about 1.1 µs per read, which is the rounding, paid once per group read;
- **568 bytes per group**, against 8 for an f64.
The size is the real cost. A document with 100,000 `Sum` groups would hold about 57 MB
of accumulators. A sparse form that stores only the occupied span of limbs would be a
few limbs for typical data. **That is an estimate, not measured.**
**No existing crate was usable.** `xsum` 0.1.6 (MIT, no dependencies, and it builds for
`wasm32-unknown-unknown`) was checked and rejected on three counts:
- **It is wrong on a single negative value.** It rounds about a quarter of single
negative inputs to a neighbouring f64: 49,819 of 200,000 sampled values, all
negative, for example `-5890.88` read back as `-5890.880000000001`.
- It keeps one infinity and one NaN flag, not counts, so a non-finite value cannot be
retracted.
- Its accumulator is not `Clone`, and reading it takes `&mut self`.
`accurate` 0.4.1 has four required dependencies, plus optional `rayon`. `fsum` 0.1.4 was found and not
evaluated. So the accumulator is about 150 lines inside `zset`, with its own tests and
a differential test against an exact rational sum.
**Defer the lossy operations: the general principle.** Keep intermediate results in a
domain where the operation is exact, and apply the lossy operations at the edge:
- sums go in the accumulator;
- decimal add, subtract and multiply are exact under ADR 005;
- division, `POW` with a non-integer exponent, and rounding for display happen only
where a value leaves the exact domain.
This is honest only about where the gain comes from. **The representation earns it, not
DBSP.** What the circuit adds is that the edges are explicit: an aggregate is a named
operator with a chosen algebra, so "where does rounding happen" has one answer per
operator rather than one per call site.
#### What each exact domain covers
| | exact accumulator | ADR 005 decimals |
|---|---|---|
| `SUM`/`AVG` over stored f64 cells | yes: exact sum, one rounding | yes, if cells are decimals |
| per-cell arithmetic (`Revenue - Cost`, `* 1.1`) | **no**: still f64 | yes |
| `=` against a typed decimal (`FLOAT_EQ_EPSILON`) | only for sums | yes, exactly |
| typed values such as money | no | yes |
| a value's own rounding when typed (`0.1`) | no: the f64 is already inexact | yes |
| cost | 568 bytes per group, about 3.2 ns per update | bignum weight in the wasm bundle (improvise-l2b8's question 3) |
**So ADR 005 still matters**, for per-cell arithmetic and for typed values like money.
The accumulator covers only `SUM` and `AVG` over f64. But that is exactly the part this
ADR needs: it is what makes the maintained `Sum` bit-equal to eager, **without waiting
for ADR 005**, provided `Eager` uses the same accumulator. That is question 3's new
option, and its recommendation.
### 9. What ADR 008 gets from it
@@ -507,7 +720,8 @@ ADR 008 (revision 3, PR #308) assumes these, under the shared names below.
not a sum of row values. So the circuit keeps a per-(column tuple, raw measure)
aggregate, and the formula layer re-evaluates only the columns whose inputs moved,
through clause 2's trait. `Expansion::members` is replaced by clause 5. Totals that
hold a `Sum` or `Avg` wait for ADR 005 (clause 8).
hold a `Sum` or `Avg` wait for step 6, whose prerequisite is question 3's answer: the
exact accumulator (recommended) or ADR 005 (clause 8).
- **What ADR 008 therefore stops asking.** Its question 6 (totals: query and memoize,
or maintain) is struck there as decided by the user through this ADR. Its questions 1
and 2 stay open there.
@@ -627,14 +841,16 @@ displayed path until the differential test holds.
map that replaces `Expansion::members` (clause 5). Formula values that depend on a
`Sum` or `Avg` read it from the eager fold until step 6. **ADR 008's retirement of
`Expansion::members` depends on this step.**
6. **`Sum` and `Avg`, and column totals, after ADR 005** (or on question 3's
fallback). A column total is a formula over the pooled fold, so it needs step 5's
6. **`Sum` and `Avg`, and column totals**, on question 3's answer: the exact
accumulator in both backends (recommended), or after ADR 005. A column total is a formula over the pooled fold, so it needs step 5's
re-evaluation. **ADR 008 clause 5 (column totals) depends on this step.**
7. **Live sources as delta streams**, when ADR 003's source trait exists.
**Tickets.** improvise-felo.4 is this ADR. The steps above are not yet filed, which is
question 10. improvise-swre stays independent and is useful under either backend.
improvise-l2b8 (ADR 005) gains a dependant, and a measured argument, from clause 8.
improvise-l2b8 (ADR 005) gains a measured argument from clause 8: eager `SUM` is
already inexact beyond `FLOAT_EQ_EPSILON`. Under question 3's recommendation it is no
longer a prerequisite of step 6.
**Good.**
@@ -652,8 +868,12 @@ improvise-l2b8 (ADR 005) gains a dependant, and a measured argument, from clause
- Two backends, kept equal by a differential test that must itself be kept honest.
- Integrated state is O(cells) in memory, per distinct None axis, in addition to the
store.
- `Sum` and `Avg` cannot be maintained on the displayed path until ADR 005 lands, and
the measurement says so.
- An f64 `Sum` cannot be maintained on the displayed path at all. Maintaining it needs
an exact representation: the accumulator costs 568 bytes per group (clause 8), or
ADR 005's decimals.
- Under question 3's recommendation, `Eager`'s displayed sums change in the last bits,
because they become correctly rounded. That is a behaviour change, even though it is
a fix.
- Circularity is decided in two layers: the formula graph's components and the
per-cell bound.
- A mutating frame from a client already clones the workbook for rollback
@@ -677,13 +897,22 @@ the user decided", as ADRs 004 and 007 did.
*Recommendation:* a seam crate. improvise-core then never depends on `zset`, and
`Eager` stays buildable and testable alone. It is the `gp-tree`/`improvise-gp` shape.
3. **Floating point.** Should ADR 005 (exact decimals) be a prerequisite for maintaining
`Sum` and `Avg` on the displayed path, with the eager fold kept for them until then?
Or should a compensated f64 sum, re-folded on detected cancellation, ship first with
a tolerance in the differential test?
*Recommendation:* ADR 005 as the prerequisite. The measured drift reaches
`FLOAT_EQ_EPSILON`, so `=` could answer differently on the two paths. No f64 technique
makes the differential test bit-exact. Everything else in this ADR lands without it.
3. **Floating point.** What makes a maintained `Sum` and `Avg` safe to display?
- (a) ADR 005 (exact decimals) as a prerequisite, with the eager fold kept for sums
until then.
- (b) A compensated f64 sum, re-folded on detected cancellation, with a tolerance in
the differential test.
- (c) **An exact accumulator** (clause 8) for `Sum` and `Avg` over f64, used by
**both** backends, so that incremental equals eager bit for bit.
*Recommendation:* (c), in `zset`. It is measured bit-equal after a million
retract-and-insert cycles, and exact against `math.fsum` on 20,000 cases. It costs
about 3.2 ns per update and 568 bytes per group, and it does not wait for ADR 005.
(b) is never bit-exact, and its gap reaches `FLOAT_EQ_EPSILON`. (a) remains right for
per-cell arithmetic and typed values, and it composes with (c): under ADR 005 a
decimal sum is exact on its own. The cost to accept is that today's displayed sums
change in their last bits, because they become correctly rounded, and that
`Avg`'s one rounding must be stated in step 6.
4. **The sixteen-pass bound.** An acyclic chain needing all sixteen passes is reported
`circular` today (core.md §3). Should `Incremental` reproduce that?
@@ -727,3 +956,16 @@ the user decided", as ADRs 004 and 007 did.
steps 4, 5 and 6?
*Recommendation:* yes. The user's decision is recorded on felo.4, and the steps are
how it is met.
11. **Ordered aggregates' key** (clause 4, for improvise-felo.7 and improvise-felo.8).
Should `LAST`/`FIRST` key their per-group state by the ordering category's item
position, with that category's item order added to the operator shape so that a
reorder rebuilds them? Or should they key by stable item identity (ADR 007), and
look the order up at read time?
*Recommendation:* item position, with the order in the operator shape. Reads stay
O(log d), like `Min`/`Max`, and the same generic operator serves all four. An item
reorder is a rare structural edit, and clause 7 already rebuilds on those. This
extends the operator shape's `measure_agg` part to cover "and the item order of
each ordering category". ADR 008 clause 3 lists the parts as `view_shape`, group
membership and `measure_agg`. This reading changes what the `measure_agg` part
covers, not the list, so ADR 008's text and cross-references still hold.
@@ -1,5 +1,8 @@
//! How far a SUM maintained by insert/retract deltas drifts from the eager fold,
//! in f64. ADR 009 (improvise-felo.4), "Floating point".
//! in f64, and what an exact accumulator costs. ADR 009 (improvise-felo.4), clause 8.
use std::hint::black_box;
use std::time::Instant;
use improvise_core::format::format_f64;
use improvise_core::model::cell::{CellKey, CellValue};
@@ -90,31 +93,196 @@ impl Compensated {
}
}
const LIMB_BITS: u32 = 32;
/// Every finite f64 is `mant * 2^(shift - 1074)` with `shift <= 2045` and `mant < 2^53`.
const LIMBS: usize = 67;
/// An exact sum of f64s: a signed fixed-point integer in units of 2^-1074, wide enough
/// for every finite f64, with non-finite values counted so they can be retracted too.
#[derive(Clone)]
struct ExactSum {
limbs: [i64; LIMBS],
pending: u32,
nan: i64,
pos_inf: i64,
neg_inf: i64,
}
impl ExactSum {
fn new() -> Self {
ExactSum {
limbs: [0; LIMBS],
pending: 0,
nan: 0,
pos_inf: 0,
neg_inf: 0,
}
}
/// Add `x` with `weight` +1 (insert) or -1 (retract). Exact: no rounding happens here.
fn add(&mut self, x: f64, weight: i64) {
match (x.is_nan(), x.is_infinite()) {
(true, _) => self.nan += weight,
(false, true) if x > 0.0 => self.pos_inf += weight,
(false, true) => self.neg_inf += weight,
(false, false) => self.add_finite(x, weight),
}
}
fn add_finite(&mut self, x: f64, weight: i64) {
let bits = x.to_bits();
let biased = ((bits >> 52) & 0x7ff) as u32;
let fraction = bits & ((1 << 52) - 1);
let (mant, shift) = match biased {
0 => (fraction, 0),
e => (fraction | (1 << 52), e - 1),
};
let sign = if x.is_sign_negative() {
-weight
} else {
weight
};
let wide = u128::from(mant) << (shift % LIMB_BITS);
let k = (shift / LIMB_BITS) as usize;
for i in 0..3 {
let part = ((wide >> (LIMB_BITS * i as u32)) & 0xffff_ffff) as i64;
self.limbs[k + i] += sign * part;
}
self.pending += 1;
if self.pending == 1 << 30 {
propagate(&mut self.limbs);
self.pending = 0;
}
}
/// The exact sum rounded to the nearest f64, ties to even: the one rounding.
fn value(&self) -> f64 {
if self.nan > 0 || (self.pos_inf > 0 && self.neg_inf > 0) {
return f64::NAN;
}
if self.pos_inf > 0 {
return f64::INFINITY;
}
if self.neg_inf > 0 {
return f64::NEG_INFINITY;
}
let mut limbs = self.limbs;
propagate(&mut limbs);
let negative = limbs[LIMBS - 1] < 0;
if negative {
limbs.iter_mut().for_each(|l| *l = -*l);
propagate(&mut limbs);
}
let magnitude = round_to_f64(&digits(&limbs));
match negative {
true => -magnitude,
false => magnitude,
}
}
}
/// Carry each limb into the next, leaving every limb but the top in [0, 2^32).
fn propagate(limbs: &mut [i64; LIMBS]) {
for k in 0..LIMBS - 1 {
let carry = limbs[k] >> LIMB_BITS;
limbs[k] -= carry << LIMB_BITS;
limbs[k + 1] += carry;
}
}
/// A non-negative normalized accumulator as 32-bit digits, the top limb split out.
fn digits(limbs: &[i64; LIMBS]) -> Vec<u64> {
let mut out: Vec<u64> = limbs.iter().map(|&l| l as u64).collect();
let mut top = out.pop().unwrap_or(0);
loop {
out.push(top & 0xffff_ffff);
top >>= LIMB_BITS;
if top == 0 {
break;
}
}
out
}
fn bit(d: &[u64], pos: usize) -> u64 {
(d[pos / 32] >> (pos % 32)) & 1
}
/// Round an exact multiple of 2^-1074 to the nearest f64, ties to even.
fn round_to_f64(d: &[u64]) -> f64 {
let Some(k) = d.iter().rposition(|&x| x != 0) else {
return 0.0;
};
let msb = 32 * k + 63 - d[k].leading_zeros() as usize;
if msb < 53 {
let int = d[0] | (d[1] << 32);
return int as f64 * f64::from_bits(1);
}
let mut mant = (0..53).fold(0u64, |m, i| m | (bit(d, msb - 52 + i) << i));
let round = bit(d, msb - 53) == 1;
let sticky = (0..msb - 53).any(|p| bit(d, p) == 1);
let mut top = msb;
if round && (sticky || mant & 1 == 1) {
mant += 1;
if mant == 1 << 53 {
mant >>= 1;
top += 1;
}
}
let exp = top as i64 - 52 - 1074;
if exp > 971 {
return f64::INFINITY;
}
let scale = match exp >= -1022 {
true => f64::from_bits(((exp + 1023) as u64) << 52),
false => f64::from_bits(1u64 << (exp + 1074)),
};
mant as f64 * scale
}
fn exact_over(values: &[i64]) -> f64 {
let mut acc = ExactSum::new();
values.iter().for_each(|&v| acc.add(typed(v), 1));
acc.value()
}
fn ulps(a: f64, b: f64) -> u64 {
(a.to_bits() as i64).abs_diff(b.to_bits() as i64)
}
struct Sums {
exact: f64,
eager: f64,
naive: f64,
comp: f64,
xs_maintained: f64,
xs_fresh: f64,
}
/// Each running sum against the eager fold (the differential oracle) and against the
/// exact decimal sum. `flip`: the gap reaches `FLOAT_EQ_EPSILON`, so `=` can answer differently.
fn report(label: &str, cycles: usize, exact: f64, eager: f64, naive: f64, comp: f64) {
fn report(label: &str, cycles: usize, s: &Sums) {
let shown = |v: f64| format_f64(v, true, 2);
println!(
"{label:<22} {cycles:>8} {:>10.2e} {:>6} {:>10.2e} {:>6} {:>5} {:>10.2e} {:>10.2e} {:>10.2e}",
(naive - eager).abs(),
ulps(naive, eager),
(comp - eager).abs(),
ulps(comp, eager),
(naive - eager).abs() >= 1e-10,
(eager - exact).abs(),
(naive - exact).abs(),
(comp - exact).abs(),
"{label:<22} {cycles:>8} {:>10.2e} {:>6} {:>10.2e} {:>6} {:>5} {:>10.2e} {:>10.2e} {:>10.2e} {:>7} {:>10.2e}",
(s.naive - s.eager).abs(),
ulps(s.naive, s.eager),
(s.comp - s.eager).abs(),
ulps(s.comp, s.eager),
(s.naive - s.eager).abs() >= 1e-10,
(s.eager - s.exact).abs(),
(s.naive - s.exact).abs(),
(s.comp - s.exact).abs(),
s.xs_maintained.to_bits() == s.xs_fresh.to_bits(),
(s.xs_maintained - s.exact).abs(),
);
if shown(naive) != shown(eager) || shown(eager) != shown(exact) {
if shown(s.naive) != shown(s.eager) || shown(s.eager) != shown(s.exact) {
println!(
" shown at ,.2: naive {} eager {} exact {}",
shown(naive),
shown(eager),
shown(exact)
" shown at ,.2: naive {} eager {} exact {} maintained-exact-sum {}",
shown(s.naive),
shown(s.eager),
shown(s.exact),
shown(s.xs_maintained)
);
}
}
@@ -126,6 +294,8 @@ fn random_updates(seed: u64, cells: usize, decades: u32, checkpoints: &[usize])
let mut naive = eager(&m);
let mut comp = Compensated::default();
comp.add(naive);
let mut xs = ExactSum::new();
values.iter().for_each(|&v| xs.add(typed(v), 1));
let mut done = 0;
for &at in checkpoints {
while done < at {
@@ -134,12 +304,21 @@ fn random_updates(seed: u64, cells: usize, decades: u32, checkpoints: &[usize])
naive = naive - typed(values[row]) + typed(new);
comp.add(-typed(values[row]));
comp.add(typed(new));
xs.add(typed(values[row]), -1);
xs.add(typed(new), 1);
values[row] = new;
m.set_cell(key(row), CellValue::Number(typed(new)));
done += 1;
}
let label = format!("{cells} cells, {decades} dec");
report(&label, at, exact(&values), eager(&m), naive, comp.value());
let sums = Sums {
exact: exact(&values),
eager: eager(&m),
naive,
comp: comp.value(),
xs_maintained: xs.value(),
xs_fresh: exact_over(&values),
};
report(&format!("{cells} cells, {decades} dec"), at, &sums);
}
}
@@ -150,9 +329,11 @@ fn absorbed_by_a_large_value() {
m.category_mut("Row").expect("just added").add_item("r0");
let mut naive = 0.0;
let mut comp = Compensated::default();
let mut xs = ExactSum::new();
m.set_cell(key(0), CellValue::Number(1e17));
naive += 1e17;
comp.add(1e17);
xs.add(1e17, 1);
for row in 1..=1000 {
m.category_mut("Row")
.expect("just added")
@@ -160,23 +341,71 @@ fn absorbed_by_a_large_value() {
m.set_cell(key(row), CellValue::Number(0.01));
naive += 0.01;
comp.add(0.01);
xs.add(0.01, 1);
}
m.clear_cell(&key(0));
naive -= 1e17;
comp.add(-1e17);
report(
"1e17, 1000 x 0.01, -1e17",
1001,
10.0,
eager(&m),
xs.add(1e17, -1);
let mut fresh = ExactSum::new();
(0..1000).for_each(|_| fresh.add(0.01, 1));
let sums = Sums {
exact: 10.0,
eager: eager(&m),
naive,
comp.value(),
comp: comp.value(),
xs_maintained: xs.value(),
xs_fresh: fresh.value(),
};
report("1e17, 1000 x 0.01, -1e17", 1001, &sums);
}
/// Nanoseconds per operation over `n` values, best of three runs.
fn per_op(n: usize, mut op: impl FnMut() -> f64) -> f64 {
(0..3)
.map(|_| {
let t = Instant::now();
black_box(op());
t.elapsed().as_secs_f64() * 1e9 / n as f64
})
.fold(f64::INFINITY, f64::min)
}
fn cost() {
let n = 10_000_000;
let mut rng = Rng(7);
let values: Vec<f64> = (0..n).map(|_| typed(decimal(&mut rng, 9))).collect();
let plain = per_op(n, || {
values
.iter()
.fold(0.0, |s, &v| black_box(s + v) - black_box(v) + v)
});
let exact_add = per_op(n, || {
let mut xs = ExactSum::new();
values.iter().for_each(|&v| {
xs.add(black_box(v), 1);
xs.add(black_box(v), -1);
xs.add(v, 1);
});
xs.limbs[0] as f64
});
let mut xs = ExactSum::new();
values.iter().for_each(|&v| xs.add(v, 1));
let read = per_op(1, || xs.value());
println!();
println!("cost, {n} values, best of 3 (per op = one insert + one retract + one insert):");
println!(" plain f64 add/sub {plain:>8.2} ns per op");
println!(" exact accumulator {exact_add:>8.2} ns per op");
println!(" exact read (round) {read:>8.0} ns per read");
println!(
" state per group {} bytes, against 8 for an f64",
std::mem::size_of::<ExactSum>()
);
}
fn main() {
println!(
"{:<22} {:>8} {:>10} {:>6} {:>10} {:>6} {:>5} {:>10} {:>10} {:>10}",
"{:<22} {:>8} {:>10} {:>6} {:>10} {:>6} {:>5} {:>10} {:>10} {:>10} {:>7} {:>10}",
"case",
"cycles",
"naive-eag",
@@ -186,7 +415,9 @@ fn main() {
"flip",
"eag-exact",
"nai-exact",
"cmp-exact"
"cmp-exact",
"xs=fresh",
"xs-exact"
);
let checkpoints = [1_000, 10_000, 100_000, 1_000_000];
for seed in [1, 2] {
@@ -195,4 +426,5 @@ fn main() {
random_updates(seed, 18_000, 9, &checkpoints);
}
absorbed_by_a_large_value();
cost();
}