diff --git a/Cargo.lock b/Cargo.lock index 77a3c5b8..51a21dda 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -97,6 +97,12 @@ dependencies = [ "num-traits", ] +[[package]] +name = "arrayvec" +version = "0.7.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3fb67a6e08acf24fdeccbac2cb6ac4305825bd1f117462e0e6f2f193345ad56" + [[package]] name = "atomic" version = "0.6.1" @@ -178,6 +184,30 @@ dependencies = [ "hybrid-array", ] +[[package]] +name = "borsh" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "553c5d846a6ba5150c65e3b1b8ec073bcf1abc20f9b7220de384a4443ea4e20a" +dependencies = [ + "borsh-derive", + "bytes", + "cfg_aliases", +] + +[[package]] +name = "borsh-derive" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "12cdfe656708a01f89b451a7d36466e6fe6c414de0aa18fc54f864f6f9ca9f56" +dependencies = [ + "once_cell", + "proc-macro-crate", + "proc-macro2", + "quote", + "syn 3.0.6", +] + [[package]] name = "bumpalo" version = "3.20.3" @@ -705,6 +735,17 @@ dependencies = [ "version_check", ] +[[package]] +name = "getrandom" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + [[package]] name = "getrandom" version = "0.3.4" @@ -1216,6 +1257,16 @@ dependencies = [ "minimal-lexical", ] +[[package]] +name = "num-bigint" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367" +dependencies = [ + "num-integer", + "num-traits", +] + [[package]] name = "num-conv" version = "0.2.2" @@ -1233,6 +1284,26 @@ dependencies = [ "syn 2.0.119", ] +[[package]] +name = "num-integer" +version = "0.1.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ce2d95d4b3734dc35aa2f45e1aa22cd416814592a4f9d9205e11affd5b8e10b" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + [[package]] name = "num-traits" version = "0.2.19" @@ -1251,6 +1322,18 @@ dependencies = [ "libc", ] +[[package]] +name = "numtower" +version = "0.1.0" +dependencies = [ + "num-bigint", + "num-integer", + "num-rational", + "num-traits", + "proptest", + "rust_decimal", +] + [[package]] name = "once_cell" version = "1.21.4" @@ -1449,6 +1532,15 @@ dependencies = [ "zerocopy", ] +[[package]] +name = "proc-macro-crate" +version = "3.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e67ba7e9b2b56446f1d419b1d807906278ffa1a658a8a5d8a39dcb1f5a78614f" +dependencies = [ + "toml_edit", +] + [[package]] name = "proc-macro2" version = "1.0.107" @@ -1469,7 +1561,7 @@ dependencies = [ "bitflags 2.13.2", "num-traits", "rand 0.9.5", - "rand_chacha", + "rand_chacha 0.9.0", "rand_xorshift", "regex-syntax", "rusty-fork", @@ -1510,6 +1602,8 @@ version = "0.8.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e058c7de0b26af77780c769414d6257830bb240f3c38477dbc2c16e5f54d6d4c" dependencies = [ + "libc", + "rand_chacha 0.3.1", "rand_core 0.6.4", ] @@ -1519,7 +1613,7 @@ version = "0.9.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" dependencies = [ - "rand_chacha", + "rand_chacha 0.9.0", "rand_core 0.9.5", ] @@ -1534,6 +1628,16 @@ dependencies = [ "rand_core 0.10.1", ] +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core 0.6.4", +] + [[package]] name = "rand_chacha" version = "0.9.0" @@ -1549,6 +1653,9 @@ name = "rand_core" version = "0.6.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom 0.2.17", +] [[package]] name = "rand_core" @@ -1713,6 +1820,23 @@ version = "0.8.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" +[[package]] +name = "rust_decimal" +version = "1.43.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7653272e75dcac41dc199fbea6f5797633994fafd339943c06c9af16bf29cd3a" +dependencies = [ + "arrayvec", + "borsh", + "bytes", + "num-traits", + "rand 0.8.8", + "rand 0.9.5", + "serde", + "serde_json", + "wasm-bindgen", +] + [[package]] name = "rustc_version" version = "0.4.1" @@ -2161,6 +2285,36 @@ dependencies = [ "tungstenite", ] +[[package]] +name = "toml_datetime" +version = "1.1.2+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b86d767906c6c42421dcba507eb9d203e779497710a47782a224bb871653053" +dependencies = [ + "serde_core", +] + +[[package]] +name = "toml_edit" +version = "0.25.17+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3641d5bbb5349a79e1020a242d251efbc546ad8048d133958323ce9c40a9c9c" +dependencies = [ + "indexmap", + "toml_datetime", + "toml_parser", + "winnow", +] + +[[package]] +name = "toml_parser" +version = "1.1.5+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baa693a8032d7e1cada7d0041e96126df243179ff061456783ac7f12bda4744c" +dependencies = [ + "winnow", +] + [[package]] name = "tungstenite" version = "0.30.0" @@ -2525,6 +2679,15 @@ dependencies = [ "windows-link", ] +[[package]] +name = "winnow" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23b97319f7b8343df12cc98938e5c3eb436064524c8d2b4e30a1d3a36eecdf81" +dependencies = [ + "memchr", +] + [[package]] name = "wit-bindgen" version = "0.57.1" diff --git a/Cargo.toml b/Cargo.toml index 1a5467a7..a07fe043 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,5 +1,5 @@ [workspace] -members = [".", "crates/improvise-formula", "crates/improvise-core", "crates/improvise-protocol", "crates/improvise-io", "crates/improvise-command", "crates/improvise-wasm-client", "crates/improvise-ws-server", "crates/improvise-worker-server", "crates/gp-tree", "crates/improvise-gp"] +members = [".", "crates/improvise-formula", "crates/improvise-core", "crates/improvise-protocol", "crates/improvise-io", "crates/improvise-command", "crates/improvise-wasm-client", "crates/improvise-ws-server", "crates/improvise-worker-server", "crates/gp-tree", "crates/improvise-gp", "crates/numtower"] [package] name = "improvise" diff --git a/context/code/numtower.md b/context/code/numtower.md new file mode 100644 index 00000000..ea27d290 --- /dev/null +++ b/context/code/numtower.md @@ -0,0 +1,328 @@ +# `numtower`: the reasoning + +`crates/numtower` is ADR 005 migration step 1 (improvise-59nw.1): the numeric tower, +first and alone. Nothing in improvise depends on it yet, so it changes no behaviour. +The decisions are the ADR's (`context/adr/005-exact-decimal-arithmetic.md`). This file +records how the crate reads them, the places the ADR left a choice open or said two +things, and the measurements. The evidence for each test is in +`context/test-evidence/numtower.md`. + +## The shape + +| type | what it is | ADR | +|---|---|---| +| `Num` | the stored value: an exact `Decimal`, an exact `Ratio`, or an `f64` | clauses 1, 2 | +| `NumRef` | a borrowed view of a `Num`'s arm, for exhaustive matching | — | +| `Interim` | a value inside one evaluation; exact intermediates widen to a bigint rational | question 13 (b) | +| `NumSum` | a retractable, order-free `SUM` over any `Num` | clause 8, question 16 | +| `ExactSum` | ADR 009's f64 superaccumulator | question 15 (a) | +| `parse_number`, `is_number`, `Typing` | clause 3's number language, gate and parse | clause 3, question 12 | +| `Display` / `Num::from_canonical` | the canonical text: `12.50`, `10/3`, `#i0.1`, `inf`, `-inf`, `nan` | clauses 10, 11, 14 | +| `Num::round_to_places` | display rounding, the mode a parameter | clause 12 "Format" | + +**`Num` is a struct with private arms, not the ADR's `pub enum`.** The ADR's sketch +(clause 2) shows `pub enum Num { Exact(Decimal), Fraction(Box>), +Inexact(f64) }`. A public enum would let any caller build values the tower never +produces: `Fraction(1/2)` (which must be the decimal `0.5`), an unreduced fraction, or an +exact `-0`. Each would break exact equality or the fingerprint of clause 9. So the arms +are private, every constructor normalizes, and callers match on `NumRef`. That is +`context/design-principles.md` §3, smart constructors. It is also what clause 14's table +asks for, since improvise calls "`Num`'s accessors for exactness, mantissa, scale and +fraction parts". The `Box` lives in one constructor, as `AGENTS.md` asks. `Num` is 24 +bytes, as the ADR measured, and `a_num_is_24_bytes` pins it. + +**The normal form.** An exact zero has no sign. A fraction is reduced and has a positive +denominator. A value whose reduced denominator is `2^p·5^q`, with `k = max(p, q) <= 28` +and a mantissa that fits 96 bits, is always the `Decimal`, never a `Fraction`. So +`Fraction(n/1)` exists only when `n` does not fit 96 bits, as with `10^35/1`. Every `i128` +is a valid numerator, `i128::MIN` included, because the arm is `Ratio`. That one +value cannot be negated in `i128`, so `Operand::parts` keeps it out of the `i128` fast +path and its arithmetic takes the bigint path (Copilot on PR #337). + +## Arithmetic + +### Why there are two types: `Num` and `Interim` + +Question 13 (b) bounds only the **stored** result: "every intermediate inside one formula +evaluation is unbounded up to the ceiling". The bounded `Num` cannot be such an +intermediate, so `Interim` is the widened type. It is either a stored `Num` or an exact +`BigQ` (a reduced `BigInt` fraction) that does not fit `i128`. `Interim::store` applies +the stored bound. The operators on `&Num` are one-operation evaluations: each builds the +`Interim`, applies the operation and stores the result. So `(a * a) / a` returns +`a` exactly when computed through `Interim`, and is inexact when computed through `Num` +(`an_intermediate_beyond_i128_comes_back_exact`). That difference is (b)'s whole point. +Step 3's evaluators should evaluate through `Interim` and store once per cell. + +**So `Interim` compares, and knows zero, without storing.** `IF(a*a/a = a, …)` and the +division guard both run before the stored bound. If `=` had to store its operands first, +its answer would depend on how the formula is spelled, which is what question 13 (b) +rules out. `Interim::cmp_exact` is clause 6's exact order over every pairing: two stored +values go through `Num::cmp_exact`, and anything wide compares as a bigint rational, +against an f64 by its exact binary value. `Interim::is_zero` is the one zero predicate. A +wide value is never zero, because zero always fits a stored arm. The comparison through +storing is the rival (`an_interim_compares_exactly_without_storing_catches_storing_first`): +`a` and `a ± 1/(3·10^30)` store to one f64 and compare equal. The review of `302594bf` +asked for this or for a note naming it as a step-3 item; the crate does the first. + +### Contagion, and how many roundings a demotion has + +An inexact operand makes the result inexact. Each exact operand is converted by the +correctly rounded conversion, and then the f64 operation runs. That is two roundings, +which is clause 4's "computed in f64". + +**An exact ∘ exact result that does not fit is rounded once.** It is the correctly rounded +f64 of the exact value. Clause 4 still says a demoted sum or product is "rounded twice, +once per operand and once in the operation". That text predates option (d). Option (d)'s +"Operations" and question 13 (b) both say that a result failing the stored bound is "the +correctly rounded f64 of the exact fraction", and under (b) the exact value is always in +hand. The crate follows (d), because (d) is decided and later. Both texts cannot hold for +one operation, so this is a reading, and it is flagged. + +### Clause 4: the hint, the recompute, the stored scale + +`Decimal + Decimal` and `Decimal × Decimal` first try `rust_decimal`'s checked operation. +The result is trusted only when its scale equals the nominal one (`max(sa, sb)` for `+` +and `−`, `sa + sb` for `×`). A lower scale is the **hint**, and so is a `None`. On the +hint the exact integer result is recomputed in a 256-bit integer, and its trailing zeros +are stripped. If the stripped mantissa fits 96 bits at a scale of at most 28, the result +is a decimal at the **stored scale**, `min(nominal, 28)`, padded back with zeros while the +mantissa fits. Otherwise the result goes to the fraction arm, then to a wide `Interim`. +The ADR's own example holds: `1.2500000000000000 × 4.0000000000000000` is `5` at scale 28 +(`a_recomputed_result_keeps_the_stored_scale_clause_4_names`). A `None` counts as the +hint, so that an overflow `rust_decimal` reports is never mistaken for an exact result. + +### Division: the direct test, one implementation + +Clause 5 is normative: write the quotient as a reduced `n / (2^p·5^q)`, then it is a +decimal iff `k <= 28` and the padded mantissa fits. The crate does this the same way for +every exact pairing, not only for decimals. It reduces the quotient (in `i128` when that +cannot overflow, and in a bigint otherwise), then normalizes. `small_ratio` and +`decimal_of` are the one normalization, so a decimal reached through any path has the +same representation. Multiplying back is never used. + +**A quotient carries its minimal scale** (`10.00 / 2` is `5`, `1 / 4` is `0.25`). The ADR +does not say what scale a quotient has. A quotient's scale is not a function of its +operands' scales (clause 5), and a result that passed through the fraction arm carries +its minimal scale (option (d)), so the minimal scale is the consistent choice. It is a +reading, and it is flagged. + +**The ADR's argument against multiplying back is weaker than it says.** Clause 5 rejects +it as incomplete, on `1e-28 / 0.5`: `q × 0.5` "needs scale 29, so the check flags it". +That holds under revision 6's rule, which demotes on the hint +(`division_decides_fit_directly_catches_multiplying_back` keeps that rival). But under +revision 7's recompute the product `1e-28` is found exact after its zeros are stripped. +So multiplying back by `Num`'s own `*` accepts the ADR's counterexample. In 20 runs of the +division property (5,120 random cases plus the three seeds), it agreed with the direct +test every time. That was an experiment and is not committed. The direct test stays the +rule, because the ADR decided it and because it needs no second division. + +### `POW` + +- **The exponent.** It is exact when it is an exact **integer value**, so `2.0` counts. +- **The size, decided before any multiplication.** The result's numerator has at least + `|n|·(bits(num) − 1)` bits, and likewise its denominator. If either lower bound passes + the ceiling, the result is `powf`, inexact, with nothing allocated. So + `POW(1/3, 10^9)` costs nothing. Inside the bound the result is computed exactly, and it + can still exceed the ceiling (`3^6000` is 9,510 bits). It then demotes, which is the + ceiling's stated exception. +- **The scale.** A positive power of a decimal takes clause 4's stored scale with nominal + `s·n`, the scale repeated multiplication would give, so `1.50^2` is `2.2500`. A + negative power is a division and has its minimal scale, **unit bases included**: + `1.0^-2` is `1` and `(-1.0)^-3` is `-1`, as `2.0^-2` is `0.25`. Until the review of + `302594bf` the unit-base arm padded every power to `s·|n|`, so `1.0^-2` read `1.00` + (Copilot on PR #337). The power property now checks the scale of every decimal result + against this rule, with unit bases among its seeds. +- **The edge values.** `0^0` is `1`. `0^-n` is `div/0`. `±1^n` is `±1` for any `n`, + including one too large to compute. +- **A demoted power takes its sign from the exact exponent.** `powf` gets the base's + magnitude, and the sign is set from `base < 0 && n odd`. An exponent above 2^53 rounds + to an even f64, so handing `powf` the signed base made `(-2)^(2^53 + 1)` read `+inf` + (Copilot on PR #337). + +### The ceiling + +`CEILING_BITS = 8,192`. It applies to the bit length of an intermediate's numerator *or* +denominator. The ADR states the ceiling as a width, not as which of the two it measures. +Past it, the intermediate becomes the correctly rounded f64 of its exact value. + +## Comparison + +- **Exact with exact.** Cross-multiplication in 256 bits. A decimal is `m / 10^s` with + `10^28 < 2^94`, so every product fits. +- **Exact with inexact.** The ADR's method: convert the exact operand, and the rounded + value decides unless it lands on the f64 itself. A tie is resolved against the f64's + exact binary value, in a bigint. +- **That fast path is sound only because the conversion is correctly rounded.** Injecting + `rust_decimal`'s `to_f64` failed `exact_mixed_comparison_is_transitive` as well as the + conversion test. A future integer conversion (clause 4 allows one) must keep that. + +`MIN` and `MAX` skip NaN. Equal values break their tie as follows. Only the first rule is +the ADR's (question 11); the other two are this crate's, and they make the result +independent of argument order: +- the exact value is preferred to the inexact one (question 11); +- of two equal decimals, the larger scale is kept, for the same reason: it keeps more of + what was typed; +- of `0.0` and `-0.0`, `MIN` takes `-0.0` and `MAX` takes `0.0`, as IEEE 754's + `minimum` and `maximum` do. + +## The number language + +**The gate is literally `f64::from_str` reading a finite value.** The parse is the +crate's own lexer, and it never goes through an f64 on an exact path. The property test +checks that the two agree, and that the exact value read rounds to the f64 +`f64::from_str` reads. + +- **A lexeme's scale.** It is the scale the lexeme implies (`1.50` → 2, `1.5e3` → 0). When + the implied scale is above 28 but the value fits a decimal (`0.1` written with 29 + places), clause 4's stored-scale rule applies: strip the zeros, then pad toward 28. + The ADR only says "with the scale the lexeme implies", which cannot exceed 28. +- **Zero.** It is exact, unsigned, at the implied scale clamped to 0–28, so + `0e99999999999999999999` is `0`. +- **An exponent too long to read** saturates, because no exact arm can reach it. A + significand of more than 180 digits, or a leading digit outside 10^-130 to 10^40, goes + straight to f64. Those bounds are wider than the fraction arm can reach (a reduced + denominator `10^k / 2^a` is at least `2^k`, so `k < 127`). +- **The forced `n/d` arm.** Its value is a decimal, then a fraction, then **the nearest + f64**. The ADR says `n/d` "is a fraction or, once normalized, a decimal", and does not + say what a forced `n/d` too large for `i128` becomes. The crate applies clause 3's order + to it, as it does to every other spelled number. This is a reading, and it is flagged. + **Except that the nearest f64 must be finite.** A forced `n/d` whose value rounds to + infinity (`n` of 400 digits, `/1`) is text, exactly as the gate makes `1e999` text. + Clause 3: "a finite lexeme must never read as infinity". One that underflows is the + inexact zero, as `1e-999` is. `is_number` and `parse_number` share the one helper, + `forced_fraction`, so the gate and the parse cannot disagree on it. + +**The canonical text is `Display`.** An inexact finite value is `#i` followed by Rust's +`{:?}` spelling, which is the shortest that reads back (`#i1e300`, `#i0.1`). +`from_canonical` reads all six spellings back, and reads an exact spelling only as an +exact value. Clause 14's list puts "`#i…`" in the crate, while its table says "the `#i` +prefix and the grammar are improvise's". Those two contradict each other. The crate +holds the prefix, because the wire and the file both need the one spelling, and step 7 +may move the prefix out if the grammar wants it. + +## Rounding to places + +Every arm goes through one exact path: the value as a bigint fraction, +`n·10^places / d`, then rounded by the caller's mode (`HalfAwayFromZero` or `HalfEven`). +An inexact value is first read exactly from its shortest round-trip spelling, which is +improvise-felo.9's fix, so a typed-then-computed `1.005` shows `1.01`. A result that +rounds to zero has no sign (`-0.001` at two places is `0.00`). The ADR does not decide +that, and today's f64 path would print `-0.00`. + +## Sums + +`NumSum` keeps five things, and every one is an abelian group, so the state is a function +of the multiset of terms and not of their order. + +| part | holds | width | +|---|---|---| +| the decimal partition | every decimal, and every fraction whose denominator is `2^p·5^q` with `k <= 28`, as an integer at scale 28 | 320 bits | +| the groups | per other reduced denominator, the summed numerator | 256 bits each, in a `BTreeMap` | +| `group_floats` | an `ExactSum` of each group's correctly rounded value | ADR 009's | +| `inexact` | an `ExactSum` of the f64 terms | ADR 009's | +| `inexact_terms` | the net weight per f64 bit pattern, entries at zero dropped | a `BTreeMap` | + +**The read.** It is exact when no inexact term is present and at most K = 8 groups +remain. The groups and the decimal partition are then combined **over the lcm of their +denominators, and reduced once**, and the stored bound applies to the result. Past K the +read is the f64 accumulators: `group_floats`, the decimal partition converted once, and +the inexact terms. That reading rounds each group, then the total. With inexact terms and +at most K groups, the exact total is converted once and added to the inexact terms. A +group whose numerator returns to zero is dropped, so "how many groups" is a function of +state. + +**"Is an inexact term present" is asked per value, never of the net weight.** The first +version kept one net count. Inserting `1.0` and retracting `2.0` leaves the count at zero +while the accumulator holds −1.0, so the read took the exact branch and dropped it +(ChatGPT Codex, P1 on PR #337). Presence is now a map from f64 bit pattern to net weight, +dropped at zero, the same rule the groups follow. It costs O(distinct inexact values), +the same order as the groups. The K cutoff overrides representability: nine denominators whose exact total fits +still read inexact (`more_than_k_denominators_read_inexact_even_when_the_total_fits`). +That is question 13's second exception, and the ADR's flagged reading. + +**`ExactSum` is ADR 009's prototype**, from `crates/improvise-core/examples/retraction_drift.rs` +on PR #309. It has three changes: +- it splits a weight into base-2^30 digits, as ADR 009 says the `numtower` version must, + instead of looping per chunk; +- it holds 70 limbs, not 67, because a digit's shift adds up to 60 bits; +- it reads through the crate's one rational-to-f64 rounding. + +It counts NaN and the infinities, so they can be retracted. A count that is nonzero in +either direction is treated as present. + +**`ExactSum`'s `==` compares sums, not representations.** A derived equality compared the +raw limbs and the `pending` carry counter, so `0.5 + 0.5 != 1.0`, an insert then retract +was not equal to an empty sum, and the same weights added in another order compared +unequal (Copilot on PR #337). For a type whose whole point is order-freedom, and which +ADR 009's group newtype will wrap, that equality is a trap. It now compares the +propagated limbs, which are canonical (every limb but the top in `[0, 2^32)`), plus the +three special-value counts, and ignores `pending`. `NumSum` derives no equality: compare +its `read()`. + +**For ADR 009's seam crate.** `NumSum` and `ExactSum` both have the shape a Z-set group +needs: `add(&Num, i64)` / `add(f64, i64)`, `insert`, `retract`, `read`/`value`, `Clone`. +`ExactSum` also has `absorb`, to merge two accumulators. Neither implements a group trait +yet, because `zset` does not exist yet. The newtype of question 15 is where that goes. + +## The wide integers + +`Fixed` is a two's-complement integer of `N` 64-bit limbs that never allocates. It is +used for the clause 4 recompute (256 bits), exact comparison (256 bits), the decimal +partition (320 bits) and the group numerators (256 bits). Clause 14 also names "256-bit +helpers for the division reduction". Division reduces in `i128` when the products cannot +overflow, and otherwise in a bigint. It does not have a hand-rolled 256-bit gcd. That +path is the slow path, and it is measured below. + +## Dependencies + +The normal edges are `rust_decimal` 1.43.0 (std only), `num-rational` 0.4.2 (no default +features), `num-bigint` 0.4.8, `num-integer` 0.1.47 and `num-traits` 0.2.19. The last two +are direct dependencies, used for `gcd`/`lcm` and the numeric traits. They were already +in the tree through `num-rational` and `num-bigint`, so they add no crate. The ADR names +only the first three. There is no serde, and `cargo tree --edges normal` shows no +`borsh`. `BigRational` is a **dev**-dependency feature only, for the oracle, so the +oracle and the implementation do not share a rational type. + +## Measured + +On an Apple M4 Max (64 GB), rustc 1.98.1, `cargo run --release -p numtower --example +bench`. Each range is two runs. + +| what | figure | +|---|---| +| SUM, 18k and 90k typed two-place values, ns per element | f64 0.53; `rust_decimal` fold 2.7–3.0; `Num` `+` fold 7.8–8.2; `NumSum` insert and read 7.9–8.1 | +| per operation over 18,000 typed pairs, ns | `Num` add 3.1–3.3, mul 3.0–3.2, div 146–153; `rust_decimal` checked add 1.4–1.5, mul 1.2, div 24.7–25.3 | +| `(Rev − Cost) / Rev`, ns per record | 151–153 | +| `to_f64` of a decimal (the text path), ns | 36.7–37.1 | +| `round_to_places(2)`, ns | 370–382 | +| `NumSum` read, K distinct 126-bit denominators, numerators of about 190 bits | K = 2: 5.2–5.4 µs; K = 8: 41.4–42.0 µs | +| `ExactSum`, 1M values | 1.78–1.81 ns per insert, 0.29 µs per read | + +Here is how these sit against the ADR: +- **The grouped read.** At K = 8 it is inside the ADR's 52–127 µs. The ADR measured + 250-bit numerators, and these are about 190 bits, because `NumSum` can only reach a + group's numerator through its public `add`. +- **`ExactSum`.** It agrees with ADR 009's 2.0 ns per update and 0.58 µs per read. +- **Division.** At about 150 ns it is above the ADR's *estimate* of 60–110 ns. Most + quotients of typed values do not terminate, so each one pays two `i128` gcds and the + `Box` of a `Fraction`. +- **The rest.** The ADR's whole-recompute measurement at 18k and 90k cells waits for + step 3, as the ADR says. + +**The wasm32 cost** comes from a probe outside the repository, built exactly as the +ADR's appendix builds one: a `cdylib`, `opt-level = "s"`, LTO, `wasm-opt -Oz +--enable-bulk-memory --enable-nontrapping-float-to-int`, `gzip -9`. The baseline already +parses and prints an f64. + +| probe | after `wasm-opt -Oz` | gzipped | +|---|---|---| +| f64 only | 48,331 B | 29,162 B | +| the ADR's five operations (parse, `+`, `−`, `/`, print) through `numtower` | +89,375 B | **+36,714 B** | +| plus `NumSum`, `cmp_exact` and `round_to_places` | +119,277 B | **+48,221 B** | + +The narrow figure confirms the ADR's estimate of "up to about +36 KB together". The full +surface is about +47 KB gzipped. That is what the bundle pays once step 4 sums through +the crate and step 6 rounds through it. It is more than the +15.4 KB the user answered +question 8 against, and the ADR says question 8 "reopens" if the figure changes the +answer. Neither figure is the real bundle, which shares allocator and formatting code +with everything else. diff --git a/context/test-evidence/numtower.md b/context/test-evidence/numtower.md new file mode 100644 index 00000000..c57949a3 --- /dev/null +++ b/context/test-evidence/numtower.md @@ -0,0 +1,468 @@ +# Test evidence: `numtower` + +What each test guards, which rival it was run against, and where that rival failed. The +entries are keyed by exact test name, and a rename orphans its entry, which +`scripts/check-comment-runs.py --check-evidence` catches. The reasoning these tests defend +is in `context/code/numtower.md`. The decisions are ADR 005's. + +**How the rivals are kept.** Every property in `crates/numtower/tests/properties.rs` is a +checker that takes the implementation as a parameter. The property's own test runs it on +the real code. A companion test, named `_catches_`, runs the same checker +on the rival, with a deterministic seed, and asserts two things: +- the run fails; +- the failure message carries **the property's own assertion text**, so a rival that + fails in setup, panics, or fails some other assertion does not count (`assert_caught` + in `tests/oracle/mod.rs`). + +So the rival is watched failing on every run of the suite, not once by hand. Two rivals +were also injected into the crate's source, and the run was recorded. Those runs are noted +under their entries. + +**The oracle** is `num-rational`'s `BigRational`, a dev-dependency feature. The crate +itself never uses `BigRational`, so the oracle shares no rational type with the code under +test. "Correctly rounded" is checked in the review's interval form (`is_nearest`). The +exact value must lie inside the round-to-nearest-even interval of the returned f64, with +the endpoints included only when the f64's mantissa is even. Comparing against +`BigRational::to_f64` would have been circular in the way ADR 005's option (d) section +describes. + +--- + +## `crates/numtower/tests/properties.rs` + +### an_exact_result_is_never_a_rounded_one + +Every exact pairing of arms and every operation (`+ − × ÷`) is checked against the oracle. +When the result claims to be exact, its value must be the exact value. The generator +spans 96-bit mantissas, scales 0–28 and `i128` fractions. A fixed case pins ADR 005's own +example, `1000 + 0.3333333333333333333333333333`. `rust_decimal` returns that at scale +25, rounded, and the crate must not. + +**Injected in the source.** In `arith.rs::decimal_sum`, the scale filter was deleted, so +`checked_add`'s result was trusted. Three tests failed: this one, +`contagion_is_total`, and `a_recomputed_result_keeps_the_stored_scale_clause_4_names`. +The revert was followed by a full rerun. + +### an_exact_result_is_never_a_rounded_one_catches_a_silently_rounding_add + +The rival is `rust_decimal`'s own `checked_add`, trusted as exact. That is the one +outcome ADR 005 clause 4 exists to prevent. + +### contagion_is_total + +Every pairing of the three arms, NaN and the infinities included, through `+ − × ÷`: +- an inexact operand makes the result inexact; +- an exact pair gives an exact result **iff** the oracle value fits a stored arm, and + otherwise the correctly rounded f64 of that value; +- division refuses iff the divisor is exactly zero, at either level. + +"Total" means no pairing is missing. A pairing the code never handled would panic, and the +runner reports a panic as a failure. + +### contagion_is_total_catches_an_integral_inexact_kept_exact + +The rival is a "harmless" promotion: an integral f64 is treated as exact, so `2 + 3.0` is +exact. It is caught on the seed pair `(2, Inexact(3.0))`. That pair is in the generator on +purpose. Integral f64s are rare among random floats, so without the seed the rival was +caught only by luck. + +### pow_is_exact_for_an_exact_base_and_integer_exponent_when_it_fits + +Bases from every arm, exponents −6 to 6. An exact base gives an exact result iff the +oracle's power fits, and otherwise the correctly rounded value. `0^-n` refuses. + +### pow_decides_its_size_before_computing_and_honours_the_ceiling + +- `POW(1/3, 10^9)` returns at once, inexact. Its lower bound is past the ceiling, so + nothing is allocated. Without the size check the test does not finish. +- `(1/3)^6000` is computed, but is 9,510 bits and demotes. That is the ceiling's + exception. +- `1.05^28` is exact and `1.05^29` is not, the option (d) table's k = 29. +- Also pinned: `(-1.0)^(10^20)`, `0^0`, the scale of `1.50^2`, and that a non-integer + exponent is `powf`. + +### the_f64_conversion_agrees_with_parsing_the_decimal_text + +The generator spans **96-bit mantissas and scales 0–28**, as ADR 005 requires. The oracle +text is built digit by digit, not by `rust_decimal`'s `Display`, which is what the crate +itself calls. So the test does not compare the crate's path with itself. The review's +case `63294.607176721429540940386554` is pinned. + +**Injected in the source.** `num.rs::decimal_to_f64` was replaced with `rust_decimal`'s +`ToPrimitive::to_f64`, and three tests failed: +- this one, minimal case `-0.10806611601966961` at scale 17; +- `the_number_languages_gate_and_parse_agree`, on `92240.01133710762224`; +- `exact_mixed_comparison_is_transitive`. + +The third was not expected. The comparison's fast path decides by the rounded value, and +that is sound only if the rounding is correct. + +### the_f64_conversion_agrees_with_parsing_the_decimal_text_catches_rust_decimal_to_f64 + +The rival is `rust_decimal`'s `to_f64`, which is not correctly rounded (ADR 005, B2). The +test also asserts that the rival misses on the review's case, so the pinned case is known +to discriminate. + +### rust_decimal_to_f64_agrees_on_fifteen_digit_mantissas_so_a_narrow_generator_is_vacuous + +This test **passes the rival** on purpose. With mantissas of at most 15 digits, +`rust_decimal`'s `to_f64` agrees with the text every time, at this test's 256 cases as +in ADR 005's 200,000 of 200,000. So a test whose generator stays in that range proves +nothing about the conversion, and this test is why the generator above spans 96 bits. It +ran at 4,096 cases until Copilot pointed out, on PR #337, that the repository's budget is +256 (`context/design-principles.md` §7). The 200,000 is the ADR's measurement, and this +test only pins its direction. + +If `rust_decimal` ever fixes its conversion, this test still passes. What fails then is +`the_f64_conversion_agrees_with_parsing_the_decimal_text_catches_rust_decimal_to_f64`, +whose rival would agree. That is the signal to re-read both. An earlier version of this +entry named this test instead, which was wrong. + +### a_fractions_conversion_lies_inside_the_round_to_nearest_even_interval + +Reduced fractions with numerator and denominator widths from 20 to 126 bits, both signs: +the review's independent form. It also checks the product of two such fractions, which +reaches the bigint path and is stored. That covers `ratio_to_f64` on values wider than +`i128`. + +### a_fractions_conversion_lies_inside_the_round_to_nearest_even_interval_catches_an_f64_division + +The rival is `n as f64 / d as f64`, which rounds three times. + +### exact_mixed_comparison_is_transitive + +The generator is a **crowd**: values crowded around one f64 `x`. +- `x` itself, and its two neighbours; +- the shortest decimal that spells `x`, and that decimal ± `k·10^-25`; +- a third of the shortest decimal. + +Every triple is checked for transitivity, and every pair for antisymmetry. Random values +would almost never put two exacts and an inexact on one f64, which is exactly where the +rival breaks. + +### exact_mixed_comparison_is_transitive_catches_comparing_through_f64 + +The rival is revision 1's rule: convert the exact operand to f64, then compare. It fails +as ADR 005 clause 6 says it does. `Exact(s) = Inexact(x)` and `Inexact(x) = Exact(s + +10^-25)`, while the two exacts differ. The marker is `comparison: not transitive`, so it +is the transitivity assertion that fails, not antisymmetry. + +### comparison_agrees_with_the_oracle_at_every_pairing_of_arms + +The differential half of comparison: the crowd, plus mixed random pools of decimals, +fractions and finite f64s, against `BigRational`'s order. NaN is unordered, and an +infinity orders past every exact value. + +### min_and_max_prefer_the_exact_value_and_do_not_depend_on_order + +Question 11's tie-break, plus the two the crate adds: the larger scale between equal +decimals, and IEEE's `-0.0` for `MIN` and `0.0` for `MAX`. Each pair is asserted in both +argument orders, the order-independence property the tie-breaks exist for. + +### the_accumulator_is_independent_of_order + +Random **weighted** term lists, from every arm, are summed in one order and in a shuffled +order. The results must agree in **arm, value and scale**. The terms include fractions +over a pool of fourteen non-decimal denominators, so groups receive several terms and +lists often hold more than K of them, and arbitrary fractions too. ADR 005 clause 8's +seed case, `[0.30000000000000004, 1e12, −1e12]`, is one of the generator's arms and is +also pinned exactly. It caught the review's F12 (a group's old f64 not retracted), where +the stale values make the past-K read depend on the order of updates. + +**On weights** (the review of `302594bf`, B2). Until then every `NumSum` property called +`insert` and `retract` only, so every weight was ±1. That was the cause the answer to +Codex's P1 had itself named, and the generators had not changed. Two faults passed the +whole suite: F7 (the decimal partition ignored the weight's magnitude) and F8 (a group +multiplied by `weight.signum()`). Weights now come from `weight()`: +- the units 1, −1 and 0; +- −3 to 3; +- magnitudes up to 2^40; +- `±i64::MAX` and 2^30 − 1. + +`i64::MIN` is left out so a weight always negates. + +### the_accumulator_is_independent_of_order_catches_decimal_addition + +The rival is clause 8's own rival, a fold with `Decimal`'s `+`, demoting to the tower +where it overflows. On the seed case, one order rounds `0.30000000000000004 + 1e12` to +scale 16 and the other does not. + +### the_accumulator_retracts_exactly + +Adding weighted `xs` and `ys`, then adding each `y` again with its weight negated in +reverse order, must equal a fresh sum of `xs` in arm, value and scale. A group whose +numerator returns to zero must be dropped, or a retracted denominator would still count +toward K. With weights it catches F7 and F12 as well. + +### the_accumulator_reads_the_exact_sum_or_its_stated_roundings + +The differential test of the read against a `BigRational` oracle, over weighted terms. It +states the value of all three branches: +- **No inexact term, at most K live groups.** The read is exact iff the exact total fits, + and is then that total. Otherwise it is the total's correctly rounded f64. +- **Inexact terms present.** The exact total is rounded once, then the weighted f64s are + added and the result is rounded once more. +- **Past K groups.** The nearest f64 of Σ round(live group) + round(decimal partition) + + Σ w·f. + +**Until the review of `302594bf` (B1) the past-K branch asserted only that the read was +inexact, and this entry claimed the test stated every branch. Both were wrong.** Three +faults passed the whole suite: +- F12: a group's old f64 was not retracted from `group_floats`; +- F14: the past-K read dropped the decimal partition; +- F15: the past-K read dropped every inexact term. + +The test now has the reviewer's oracle, at the repository's 256-case budget. 34ab0fc3 ran +it at 512, and Copilot pointed out the budget. At 256, every fault in the table below +still failed by name, and F7, F8, F12, F14 and F15 failed in four of four passes. Half the lists come from the +fourteen-denominator pool, up to 40 terms long, so groups receive several terms and K is +passed often. **It also counts how many cases reach each branch and fails if any branch +is reached fewer than 20 times**, so a generator change cannot silently starve one. + +The oracle counts **live** groups, those whose summed numerator is nonzero. The first +version counted denominators ever seen, and that disagrees with the state-based count +whenever terms cancel. Inexact presence is likewise per bit pattern, by net weight. + +**Every fault in the review's script, re-run after the fix.** Each was injected into the +source and `cargo test -p numtower` was run on a confirmed rebuild. Each failed by name: + +| fault | failed | +|---|---| +| F1 decimal sum's scale filter dropped | `a_recomputed_result_keeps_the_stored_scale_clause_4_names`, `an_exact_result_is_never_a_rounded_one`, `contagion_is_total` | +| F2 `rust_decimal`'s `to_f64` | `exact_mixed_comparison_is_transitive`, `comparison_agrees_with_the_oracle_at_every_pairing_of_arms`, `the_f64_conversion_agrees_with_parsing_the_decimal_text`, `the_number_languages_gate_and_parse_agree` | +| F3 mixed comparison through f64 | `exact_mixed_comparison_is_transitive`, `comparison_agrees_with_the_oracle_at_every_pairing_of_arms`, `an_interim_compares_exactly_without_storing` | +| F4 one net inexact count | `the_accumulator_keeps_inexact_terms_whose_weights_cancel`, `the_accumulator_retracts_exactly`, this test | +| F5 product, not lcm | `sum::tests::the_grouped_read_combines_over_the_lcm_and_reduces_once` | +| F6 `ExactSum` rounded per limb | `exact_sum_is_the_exact_weighted_total_rounded_once`, `more_than_k_denominators_read_inexact_even_when_the_total_fits`, `the_accumulator_keeps_inexact_terms_whose_weights_cancel`, this test | +| F7 decimal weight magnitude ignored | this test, `the_accumulator_retracts_exactly` | +| F8 group weight as its sign | this test | +| F9 decimal product's scale filter dropped | `a_typed_value_and_the_same_value_computed_are_equally_exact`, `an_exact_result_is_never_a_rounded_one`, `contagion_is_total` | +| F10 `ExactSum`'s low weight digit only | `exact_sum_is_the_exact_weighted_total_rounded_once`, this test | +| F11 `POW` sign lost | `a_huge_odd_exponent_keeps_the_sign_of_a_negative_base` | +| F12 old group f64 kept | this test, `the_accumulator_is_independent_of_order`, `the_accumulator_retracts_exactly` | +| F13 mixed tie not checked exactly | `comparison_agrees_with_the_oracle_at_every_pairing_of_arms`, `exact_mixed_comparison_is_transitive`, `an_interim_compares_exactly_without_storing` | +| F14 past K, decimal partition dropped | this test | +| F15 past K, inexact terms dropped | this test | + +- **F13's text now occurs twice.** `Interim`'s comparison added a second copy of the + replaced line, so the review's script refuses to inject it. It was injected at its + original site, in `cmp.rs::mixed`. +- **The five former survivors were run twice and failed both times.** F12's set varied + between runs: `the_accumulator_retracts_exactly` failed in one run and not the other, + and this test failed in both. + +### more_than_k_denominators_read_inexact_even_when_the_total_fits + +`1/3 + 1/7 + … + 1/31` is nine distinct denominators. The exact total fits an `i128` +fraction, yet the read is inexact. That is question 13's second exception, the K cutoff +overriding representability. Eight of them read as an exact `Fraction`. The read past K is +asserted to be **each group rounded, then the total rounded once**. The first version +asserted one rounding of the exact total, and it failed: clause 8's maintained f64 +accumulator of group values rounds twice by design. + +### exact_sum_is_the_exact_weighted_total_rounded_once + +ADR 009's `ExactSum` against the oracle. Weights are drawn from small values, from the +whole `i64` range and up to 2^40. A whole-range weight exercises all three base-2^30 +digits, which ADR 009's chunk loop could not reach in reasonable time. Values include +subnormals. + +### division_decides_fit_directly + +ADR 005 clause 5's three seed cases, as both property inputs and pinned assertions: +- `1 / 2^40` and `1 / 5^30` are fractions; +- `1e-28 / 0.5` is the exact decimal `2e-28`. + +Random decimal pairs must also be an exact decimal iff the quotient fits one, an exact +fraction iff it fits only that, and otherwise inexact. The test also pins `AVG`'s example +(`70/3`), `1/3 × 3 = 1`, and division by an inexact `-0.0`. + +### division_decides_fit_directly_catches_multiplying_back + +The rival is multiplying back with **revision 6's** scale test (demote on the hint). It +fails on `1e-28 / 0.5`, as ADR 005 says. An experiment, not committed, ran multiplying +back with `Num`'s own `*`, which carries revision 7's recompute. It agreed with the direct +test on all 5,120 random cases and the seeds. The ADR's argument that multiplying back is +incomplete holds only under the earlier rule (`context/code/numtower.md`, "Division"). + +### the_number_languages_gate_and_parse_agree + +The generated lexemes have four sources: +- structured lexemes: an optional sign, up to 40 digits, an optional fraction, and an + optional exponent of up to 4 digits, so the empty parts give `1.`, `.5` and `.`; +- junk over the alphabet of digits, `eE+-._/aifnIN` and space; +- `n/d` lexemes; +- the named spellings: `inf`, `NaN`, `Infinity`, `1_000`, `+.5`, `e5`. + +For each lexeme, the parse succeeds iff `f64::from_str` reads a finite value. The exact +value parsed must round to the f64 `f64::from_str` reads. `is_number` and `parse_number` +must agree in both typings. + +### the_number_languages_gate_and_parse_agree_catches_from_str_exact + +The rival is `rust_decimal`'s `from_str_exact`. ADR 005 clause 3 names it as the parse +whose language differs from the gate's: it rejects `1.5e3` and accepts `1_000`. + +### the_number_language_seed_cases + +Every seed case ADR 005's migration step 1 lists: +- `1.5e3` is the exact `1500`; +- `1_000` is text; +- `-0.0` is the exact `0.0`, with no sign; +- `3/4` is text when the type is inferred, and `0.75` when it is forced; +- `1e999`, `-1e999`, `nan`, `inf` and `Infinity` are text; +- `1e-999` is an inexact zero; +- `1e-30` is the exact `1/10^30`; +- `0.12345678901234567890123456789` is an exact fraction, not a rounded decimal; +- `2^120` is an exact `n/1`; +- `1e-40` is inexact. + +### a_typed_value_and_the_same_value_computed_are_equally_exact + +For random decimal pairs (and the seed `1e-15 × 1e-15`), the exact product is spelled as +`e-` and parsed. The parse must have the arm and the value of the computed +product. An earlier draft built the text only when the reduced denominator was a power of +ten, and rejected most cases. A reduced decimal product's denominator is `2^a·5^b`, so the +text now scales by the smallest power of ten that makes the product an integer. + +### a_typed_value_and_the_same_value_computed_are_equally_exact_catches_demoting_on_the_hint + +The rival is revision 6's multiplication, which demotes whenever the returned scale drops. +`1e-15 × 1e-15` then becomes an f64, while the typed `1e-30` is the exact fraction. That +is the asymmetry revision 6's N2 was about. + +### a_recomputed_result_keeps_the_stored_scale_clause_4_names + +ADR 005 clause 4's example: `1.2500000000000000 × 4.0000000000000000` is `5` at scale 28, +not scale 0 and not 32. The test also pins that a common-path sum keeps its returned scale +(`12.50 + 1.5 = 14.00`), that `1000 + 0.333…` goes to the fraction arm, and that an exact +zero result has no sign. + +### an_intermediate_beyond_i128_comes_back_exact + +The point of question 13 (b). `a = 10^30/7`; `a × a` does not fit `i128`, yet `(a × a) / +a` through `Interim` is `a`, exactly. The same product stored through `Num` is inexact +and correctly rounded. + +### the_canonical_text_reads_back_as_the_same_number + +Every arm's `Display` text reads back through `from_canonical` as the same arm, value and +scale. The test pins the four spellings ADR 005 clause 11 tabulates, and that canonical +exact text never reads back as an inexact value. + +### the_accumulator_keeps_inexact_terms_whose_weights_cancel + +ChatGPT Codex, P1 on PR #337, against `edbb50c4`. `NumSum` decided "is an inexact term +present" from the **net** inexact weight. Adding `1.0` with weight 1 and `2.0` with weight +−1 nets to zero while the accumulator holds −1.0. So `read()` took the exact branch and +returned the exact partition alone, dropping the inexact term silently. + +This test **failed at `edbb50c4`**, on its first assertion: the read of `[1.0 ×1, 2.0 ×−1, +0.5]` claimed to be exact. Two things made the bug invisible to the existing properties: +- every generated term was inserted with weight +1; +- the retraction test retracts the very values it inserted, so per-value presence and net + presence agree there by construction. + +The test also asserts that retracting both inexact terms returns the exact `0.5`, so the +fix cannot simply make every sum that ever saw an f64 inexact. + +### a_huge_odd_exponent_keeps_the_sign_of_a_negative_base + +Copilot on PR #337, against `edbb50c4`. When the size check demotes a power to `powf`, the +exponent becomes an f64. Above 2^53 it rounds to an even value, so `(-2)^(2^53 + 1)` +returned `+inf` where the sign must be negative, and the negative exponent gave `+0.0`. +This test **failed at `edbb50c4`** on its first assertion. It also checks a fraction base, +and that an even exponent stays positive, so a fix that always negates cannot pass. + +### i128_min_is_an_exact_fraction_numerator + +Copilot on PR #337, against `edbb50c4`. `Num::fit` rejected `i128::MIN` as a numerator +because `to_i128` filtered it for both halves, so `Num::ratio(i128::MIN, 1)` became inexact +although `Ratio` holds it. This test **failed at `edbb50c4`** on its first +assertion. The oracle's `fits` had the same symmetric-range assumption, so no property +could see the bug. Both are fixed, and `fraction()` now generates `i128::MIN / d`. + +Allowing the value exposes a second hazard: it cannot be negated in `i128`, and the old +`negate` computed `-n`, which panics in a debug build. So the test also negates +`i128::MIN / 3` (an inexact, correctly rounded result), takes it through `Interim` twice +and back, divides it, compares it and reads its canonical text. + +### a_negative_power_has_its_minimal_scale + +Copilot on PR #337, and B3 of the review of `302594bf`. The unit-base arm of `POW` +padded every result to `s·|n|`, so `1.0^-2` was `1.00` and `(-1.0)^-3` was `-1.000`, +while `2.0^-2` was `0.25`. `context/code/numtower.md` already said "a negative power … +has its minimal scale", so the code contradicted the crate's stated rule. This test +**failed at `302594bf`** on its first assertion. It also pins that a positive power keeps +`s·n` (`(-1.0)^3 = -1.000`), so a fix that stripped every scale cannot pass. + +The broader guard is the scale assertion now in +`pow_is_exact_for_an_exact_base_and_integer_exponent_when_it_fits`, which had checked +value and exactness but never scale. Its oracle, `power_scale`, is the rule restated: +- start from the minimal scale; +- pad toward `min(s·k, 28)` for a decimal base with `k > 0` only, while the mantissa + fits 96 bits. + +Unit bases (`1.0`, `-1.0`, `1.000`, `-1`, `0.00`) are among its seeds. It **failed at +`302594bf`** with the minimal case `1.0^-1 = 1.0`. The injection G1 (a negative power +keeps the base's scale) fails both tests. + +### exact_sums_are_equal_when_their_values_are + +Copilot on PR #337, and B3 of the review. `ExactSum` derived `PartialEq` over its raw +limbs and its `pending` carry counter. This test **failed at `302594bf`** on its first +assertion, `0.5 + 0.5 == 1.0`. It also asserts: +- insert then retract equals an empty sum; +- weights `(2^30 − 1, 2)` added in either order are equal; +- different values are unequal, and a counted NaN is part of the value, so an `==` that + always answers true cannot pass. + +Injection G2, equality back on the representation, fails it. + +### a_forced_fraction_never_reads_as_infinity + +A should-fix from the review of `302594bf`. A forced `n/d` with a 400-digit `n` read as +`Inexact(inf)`, though clause 3 says a finite lexeme never reads as infinity. This test +**failed at `302594bf`**. It also asserts: +- `is_number` agrees with the parse; +- a negative one is text too; +- an underflowing `1/9…9` is still the inexact zero, so a fix that refused every inexact + forced fraction cannot pass. + +Injection G3 fails it. + +### an_interim_compares_exactly_without_storing + +New API, from a should-fix in the review: `Interim` had no comparison and no zero test +while step 3 is told to evaluate through it. For two wide products `a`, it compares `a` +against: +- another wide product; +- `a ± 1/(3·10^30)`; +- `a` itself; +- the f64 nearest `a`, and that f64's upper neighbour. + +Each comparison is made in both directions, against the `BigRational` order. It also +asserts `a − a` is zero and `a` is not. Injection G4, storing before comparing, fails it, +and so do F3 and F13, because the stored arms delegate to `Num::cmp_exact`. + +### an_interim_compares_exactly_without_storing_catches_storing_first + +The rival is the comparison the review warned step 3 would otherwise have to make: store +both operands, then compare. `a` and `a ± 1/(3·10^30)` store to one f64, so the rival +calls them equal. + +--- + +## `crates/numtower/src/sum.rs` + +### sum::tests::the_grouped_read_combines_over_the_lcm_and_reduces_once + +The combine must return the **lcm** as its denominator, and the reduced result must equal +the `BigRational` sum. Group numerators are up to `i128 × i64`, and the denominators share +factors (3, 77, 2^20, 3^40), so the lcm and the product differ. + +### sum::tests::the_grouped_read_combines_over_the_lcm_and_reduces_once_catches_the_product + +The rival combines over the product of the denominators. Its value is right, so only the +lcm assertion can see it. That is why the property asserts the denominator before the +value. diff --git a/crates/numtower/Cargo.toml b/crates/numtower/Cargo.toml new file mode 100644 index 00000000..54ca4e27 --- /dev/null +++ b/crates/numtower/Cargo.toml @@ -0,0 +1,25 @@ +[package] +name = "numtower" +version = "0.1.0" +edition = "2024" +license = "Apache-2.0" +description = "A three-level numeric tower: exact decimals, bounded exact fractions and inexact f64, with Scheme's contagion rules and exact, order-free sums." +keywords = ["decimal", "rational", "exact", "numeric-tower"] +categories = ["mathematics"] + +# ADR 005 clause 14: rust_decimal (std only), num-rational without its default +# features (for Ratio), and num-bigint always. No serde. +[dependencies] +rust_decimal = { version = "1.43", default-features = false, features = ["std"] } +num-rational = { version = "0.4.2", default-features = false } +num-bigint = { version = "0.4", default-features = false, features = ["std"] } +num-integer = { version = "0.1", default-features = false } +num-traits = { version = "0.2", default-features = false, features = ["std"] } + +[dev-dependencies] +num-rational = { version = "0.4.2", features = ["num-bigint-std"] } +proptest = "1" + +[[example]] +name = "bench" +test = false diff --git a/crates/numtower/examples/bench.rs b/crates/numtower/examples/bench.rs new file mode 100644 index 00000000..809a3be0 --- /dev/null +++ b/crates/numtower/examples/bench.rs @@ -0,0 +1,187 @@ +//! ADR 005 migration step 1's benchmark: `cargo run --release -p numtower --example bench`. +//! A microbenchmark, not the evaluator; what it does and does not show is in `context/code/numtower.md`. + +use std::hint::black_box; +use std::time::Instant; + +use numtower::{Decimal, ExactSum, Num, NumSum, Rounding, Typing, parse_number}; + +struct SplitMix64(u64); + +impl SplitMix64 { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = self.0; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } +} + +/// Typed two-place values from 0.01 to 999,999.99, as the ADR's appendix draws them. +fn typed_texts(seed: u64, n: usize) -> Vec { + let mut rng = SplitMix64(seed); + (0..n) + .map(|_| { + let cents = 1 + rng.next() % 99_999_999; + format!("{}.{:02}", cents / 100, cents % 100) + }) + .collect() +} + +fn nums(texts: &[String]) -> Vec { + texts + .iter() + .map(|t| parse_number(t, Typing::Inferred).unwrap()) + .collect() +} + +/// The best of `runs`, in nanoseconds per item. +fn best(runs: usize, items: usize, mut f: impl FnMut()) -> f64 { + (0..runs) + .map(|_| { + let t = Instant::now(); + f(); + t.elapsed().as_nanos() as f64 / items as f64 + }) + .fold(f64::INFINITY, f64::min) +} + +fn sums() { + println!("SUM over typed two-place values, ns per element, best of 20"); + for n in [18_000, 90_000] { + let texts = typed_texts(1, n); + let floats: Vec = texts.iter().map(|t| t.parse().unwrap()).collect(); + let decimals: Vec = texts.iter().map(|t| t.parse().unwrap()).collect(); + let values = nums(&texts); + let f64_fold = best(20, n, || { + black_box(black_box(&floats).iter().sum::()); + }); + let decimal_fold = best(20, n, || { + let s = black_box(&decimals) + .iter() + .fold(Decimal::ZERO, |a, b| a.checked_add(*b).unwrap()); + black_box(s); + }); + let num_fold = best(20, n, || { + let zero = Num::exact(Decimal::ZERO); + black_box(black_box(&values).iter().fold(zero, |a, b| &a + b)); + }); + let accumulator = best(20, n, || { + let mut s = NumSum::new(); + black_box(&values).iter().for_each(|v| s.insert(v)); + black_box(s.read()); + }); + println!( + " {n:>6} cells: f64 {f64_fold:.2}, rust_decimal {decimal_fold:.2}, Num + fold {num_fold:.2}, NumSum {accumulator:.2}" + ); + } +} + +fn operations() { + println!("Per operation over 18,000 typed pairs, ns, best of 20"); + let (a, b) = ( + nums(&typed_texts(21, 18_000)), + nums(&typed_texts(22, 18_000)), + ); + let decimal = + |v: &[Num]| -> Vec { v.iter().map(|n| n.to_string().parse().unwrap()).collect() }; + let (da, db) = (decimal(&a), decimal(&b)); + let pairs = a.len(); + let rd_add = best(20, pairs, || { + da.iter().zip(&db).for_each(|(x, y)| { + black_box(x.checked_add(*y)); + }) + }); + let rd_mul = best(20, pairs, || { + da.iter().zip(&db).for_each(|(x, y)| { + black_box(x.checked_mul(*y)); + }) + }); + let rd_div = best(20, pairs, || { + da.iter().zip(&db).for_each(|(x, y)| { + black_box(x.checked_div(*y)); + }) + }); + let add = best(20, pairs, || { + a.iter().zip(&b).for_each(|(x, y)| drop(black_box(x + y))) + }); + let mul = best(20, pairs, || { + a.iter().zip(&b).for_each(|(x, y)| drop(black_box(x * y))) + }); + let div = best(20, pairs, || { + a.iter() + .zip(&b) + .for_each(|(x, y)| drop(black_box(x.checked_div(y)))) + }); + let margin = best(20, pairs, || { + a.iter() + .zip(&b) + .for_each(|(rev, cost)| drop(black_box((rev - cost).checked_div(rev)))); + }); + let to_f64 = best(20, pairs, || { + a.iter().for_each(|x| { + black_box(x.to_f64()); + }) + }); + let places = best(20, pairs, || { + a.iter() + .for_each(|x| drop(black_box(x.round_to_places(2, Rounding::HalfAwayFromZero)))); + }); + println!(" rust_decimal checked: add {rd_add:.1}, mul {rd_mul:.1}, div {rd_div:.1}"); + println!(" Num: add {add:.1}, mul {mul:.1}, div {div:.1}"); + println!( + " Num (Rev - Cost) / Rev {margin:.1}, to_f64 {to_f64:.1}, round_to_places(2) {places:.1}" + ); +} + +/// K groups of near-`i128` denominators: the ADR's worst case for the exact read. +fn grouped_read() { + println!("NumSum read over K distinct 126-bit denominators, us per read, best of 200"); + for k in [2usize, 8] { + let mut rng = SplitMix64(0x5EED); + let mut sum = NumSum::new(); + for _ in 0..k { + let den = ((rng.next() as i128) << 62 | rng.next() as i128 >> 2) | 1 | (1 << 125); + let num = (rng.next() as i128) << 60 ^ rng.next() as i128; + let term = Num::ratio(num, den).unwrap(); + sum.add(&term, i64::MAX); + sum.add(&term, i64::MAX); + } + assert_eq!(sum.groups(), k); + let read = best(200, 1, || drop(black_box(sum.read()))); + println!( + " K = {k}: {:.1} us (exact: {})", + read / 1000.0, + sum.read().is_exact() + ); + } +} + +fn exact_sum() { + println!("ExactSum (f64 superaccumulator), best of 20"); + let floats: Vec = typed_texts(1, 1_000_000) + .iter() + .map(|t| t.parse().unwrap()) + .collect(); + let mut acc = ExactSum::new(); + let insert = best(20, floats.len(), || { + acc = ExactSum::new(); + floats.iter().for_each(|x| acc.add(*x, 1)); + }); + let read = best(20, 1, || { + black_box(acc.value()); + }); + println!( + " {insert:.2} ns per insert, {:.2} us per read", + read / 1000.0 + ); +} + +fn main() { + println!("numtower bench (release); name the machine when you quote these"); + sums(); + operations(); + grouped_read(); + exact_sum(); +} diff --git a/crates/numtower/src/arith.rs b/crates/numtower/src/arith.rs new file mode 100644 index 00000000..1b19a38b --- /dev/null +++ b/crates/numtower/src/arith.rs @@ -0,0 +1,445 @@ +//! Arithmetic and contagion (ADR 005 clauses 4 and 5, question 13 (b)). +//! The reasoning for each rule is in `context/code/numtower.md`, "Arithmetic". + +use std::ops::{Add, Mul, Neg, Sub}; + +use num_bigint::{BigInt, BigUint}; +use num_integer::Integer; +use num_traits::{One, Signed, ToPrimitive, Zero}; +use rust_decimal::Decimal; + +use crate::fixed::I256; +use crate::num::{MANTISSA_LIMIT, MAX_SCALE, Num, NumRef, decimal, pow10, small_ratio}; +use crate::ratio::BigQ; + +/// The widest an exact intermediate may grow, in bits of numerator or denominator. +pub const CEILING_BITS: u64 = 8192; + +/// A value inside one evaluation: exact intermediates are unbounded up to the ceiling. +#[derive(Clone, Debug)] +pub struct Interim(Inner); + +#[derive(Clone, Debug)] +enum Inner { + Stored(Num), + Wide(Box), +} + +/// The divisor was exactly zero, exact or inexact (clause 5). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct DivisionByZero; + +impl From for Interim { + fn from(n: Num) -> Interim { + Interim(Inner::Stored(n)) + } +} + +impl Interim { + fn wide(q: BigQ) -> Interim { + Interim(Inner::Wide(Box::new(q))) + } + + /// The stored bound of question 13 (b): a wide value is the nearest f64. + pub fn store(self) -> Num { + match self.0 { + Inner::Stored(n) => n, + Inner::Wide(q) => Num::inexact(q.to_f64()), + } + } + + pub fn is_exact(&self) -> bool { + match &self.0 { + Inner::Stored(n) => n.is_exact(), + Inner::Wide(_) => true, + } + } + + fn operand(&self) -> Operand<'_> { + match &self.0 { + Inner::Stored(n) => Operand::of(n), + Inner::Wide(q) => Operand::Big(q), + } + } + + pub fn checked_div(&self, o: &Interim) -> Result { + divide(self.operand(), o.operand()) + } + + pub fn pow(&self, o: &Interim) -> Result { + power(self.operand(), o.operand()) + } + + /// Exactly zero; a wide value never is, because zero fits a stored arm. + pub fn is_zero(&self) -> bool { + match &self.0 { + Inner::Stored(n) => n.is_zero(), + Inner::Wide(_) => false, + } + } + + /// Clause 6's exact order, before the stored bound; `None` only when either is NaN. + pub fn cmp_exact(&self, o: &Interim) -> Option { + match (&self.0, &o.0) { + (Inner::Stored(a), Inner::Stored(b)) => a.cmp_exact(b), + _ => match (self.operand().big(), o.operand().big()) { + (Some(a), Some(b)) => Some(a.cmp_value(&b)), + (Some(a), None) => crate::cmp::exact_against_f64(&a, o.operand().float()), + (None, Some(b)) => { + crate::cmp::exact_against_f64(&b, self.operand().float()).map(|c| c.reverse()) + } + (None, None) => self.operand().float().partial_cmp(&o.operand().float()), + }, + } + } +} + +impl Num { + pub fn checked_div(&self, o: &Num) -> Result { + divide(Operand::of(self), Operand::of(o)).map(Interim::store) + } + + pub fn pow(&self, o: &Num) -> Result { + power(Operand::of(self), Operand::of(o)).map(Interim::store) + } +} + +macro_rules! binary_ops { + ($($trait:ident $method:ident $f:ident;)*) => {$( + impl $trait<&Interim> for &Interim { + type Output = Interim; + fn $method(self, o: &Interim) -> Interim { + $f(self.operand(), o.operand()) + } + } + impl $trait<&Num> for &Num { + type Output = Num; + fn $method(self, o: &Num) -> Num { + $f(Operand::of(self), Operand::of(o)).store() + } + } + )*}; +} + +binary_ops! { + Add add plus; + Sub sub minus; + Mul mul times; +} + +impl Neg for &Interim { + type Output = Interim; + fn neg(self) -> Interim { + negate(self.operand()) + } +} + +impl Neg for &Num { + type Output = Num; + fn neg(self) -> Num { + negate(Operand::of(self)).store() + } +} + +#[derive(Clone, Copy)] +enum Operand<'a> { + Dec(&'a Decimal), + Small(i128, i128), + Big(&'a BigQ), + Float(f64), +} + +impl<'a> Operand<'a> { + fn of(n: &'a Num) -> Operand<'a> { + match n.kind() { + NumRef::Exact(d) => Operand::Dec(d), + NumRef::Fraction(r) => Operand::Small(*r.numer(), *r.denom()), + NumRef::Inexact(x) => Operand::Float(x), + } + } + + /// The `i128` fast path's operands; `i128::MIN` cannot be negated, so it takes the bigint path. + fn parts(self) -> Option<(i128, i128)> { + match self { + Operand::Dec(d) => Some((d.mantissa(), pow10(d.scale()))), + Operand::Small(n, d) => (n != i128::MIN).then_some((n, d)), + Operand::Big(_) | Operand::Float(_) => None, + } + } + + fn big(self) -> Option { + match self { + Operand::Dec(d) => Some(BigQ::from_i128(d.mantissa(), pow10(d.scale()))), + Operand::Small(n, d) => Some(BigQ::from_i128(n, d)), + Operand::Big(q) => Some(q.clone()), + Operand::Float(_) => None, + } + } + + fn float(self) -> f64 { + match self { + Operand::Dec(d) => crate::num::decimal_to_f64(d), + Operand::Small(n, d) => crate::ratio::ratio_to_f64(&BigInt::from(n), &BigInt::from(d)), + Operand::Big(q) => q.to_f64(), + Operand::Float(x) => x, + } + } + + fn is_float(self) -> bool { + matches!(self, Operand::Float(_)) + } + + fn is_zero(self) -> bool { + match self { + Operand::Dec(d) => d.is_zero(), + Operand::Small(..) | Operand::Big(_) => false, + Operand::Float(x) => x == 0.0, + } + } +} + +fn inexact(x: f64) -> Interim { + Interim::from(Num::inexact(x)) +} + +/// An exact result in the first arm it fits; wide while under the ceiling, else demoted. +fn settle_wide(q: BigQ) -> Interim { + match Num::fit(&q) { + Some(n) => Interim::from(n), + None if q.bits() > CEILING_BITS => inexact(q.to_f64()), + None => Interim::wide(q), + } +} + +fn exact_or_wide( + a: Operand, + b: Operand, + small: impl Fn(i128, i128, i128, i128) -> Option<(i128, i128)>, + big: impl Fn(&BigQ, &BigQ) -> BigQ, +) -> Interim { + let fast = a + .parts() + .zip(b.parts()) + .and_then(|((n1, d1), (n2, d2))| small(n1, d1, n2, d2)); + fast.and_then(|(n, d)| small_ratio(n, d)) + .map(Interim::from) + .unwrap_or_else(|| match a.big().zip(b.big()) { + Some((x, y)) => settle_wide(big(&x, &y)), + None => unreachable!("both operands are exact here"), + }) +} + +fn plus(a: Operand, b: Operand) -> Interim { + match (a, b) { + (x, y) if x.is_float() || y.is_float() => inexact(x.float() + y.float()), + (Operand::Dec(x), Operand::Dec(y)) => decimal_sum(x, y), + _ => exact_or_wide( + a, + b, + |n1, d1, n2, d2| { + let n = n1.checked_mul(d2)?.checked_add(n2.checked_mul(d1)?)?; + Some((n, d1.checked_mul(d2)?)) + }, + BigQ::add, + ), + } +} + +fn minus(a: Operand, b: Operand) -> Interim { + plus(a, negate(b).operand()) +} + +fn times(a: Operand, b: Operand) -> Interim { + match (a, b) { + (x, y) if x.is_float() || y.is_float() => inexact(x.float() * y.float()), + (Operand::Dec(x), Operand::Dec(y)) => decimal_product(x, y), + _ => exact_or_wide( + a, + b, + |n1, d1, n2, d2| { + let (g1, g2) = (n1.gcd(&d2).max(1), n2.gcd(&d1).max(1)); + Some(( + (n1 / g1).checked_mul(n2 / g2)?, + (d1 / g2).checked_mul(d2 / g1)?, + )) + }, + BigQ::mul, + ), + } +} + +fn negate(a: Operand) -> Interim { + match a { + Operand::Dec(d) => Interim::from(Num::exact(-*d)), + Operand::Float(x) => inexact(-x), + Operand::Small(n, d) => n + .checked_neg() + .and_then(|m| small_ratio(m, d)) + .map(Interim::from) + .unwrap_or_else(|| settle_wide(BigQ::from_i128(n, d).neg())), + Operand::Big(q) => settle_wide(q.neg()), + } +} + +fn divide(a: Operand, b: Operand) -> Result { + match (a, b) { + (_, y) if y.is_zero() => Err(DivisionByZero), + (x, y) if x.is_float() || y.is_float() => Ok(inexact(x.float() / y.float())), + _ => Ok(exact_or_wide( + a, + b, + |n1, d1, n2, d2| { + let (g1, g2) = (n1.gcd(&n2).max(1), d1.gcd(&d2).max(1)); + Some(( + (n1 / g1).checked_mul(d2 / g2)?, + (d1 / g2).checked_mul(n2 / g1)?, + )) + }, + |x, y| { + x.div(y) + .unwrap_or_else(|| unreachable!("the divisor is not zero")) + }, + )), + } +} + +/// Clause 4's hint: a returned scale below the nominal one may mean rounding. +fn decimal_sum(x: &Decimal, y: &Decimal) -> Interim { + let nominal = x.scale().max(y.scale()); + match x.checked_add(*y).filter(|r| r.scale() == nominal) { + Some(r) => Interim::from(Num::exact(r)), + None => { + let align = |d: &Decimal| { + I256::from_i128(d.mantissa()) + .checked_mul_u128(pow10(nominal - d.scale()) as u128) + .expect("an aligned mantissa is below 2^190") + }; + let sum = align(x).checked_add(&align(y)).expect("below 2^191"); + settle_decimal_integer(sum, nominal, nominal) + } + } +} + +fn decimal_product(x: &Decimal, y: &Decimal) -> Interim { + let nominal = x.scale() + y.scale(); + match x.checked_mul(*y).filter(|r| r.scale() == nominal) { + Some(r) => Interim::from(Num::exact(r)), + None => { + settle_decimal_integer(I256::mul_i128(x.mantissa(), y.mantissa()), nominal, nominal) + } + } +} + +/// The recompute of clause 4: strip trailing zeros, fit 96 bits at scale <= 28, else a fraction. +fn settle_decimal_integer(m: I256, scale: u32, nominal: u32) -> Interim { + let (mut m, mut s) = (m, scale); + while s > 0 { + let (q, r) = m.div_rem_u64(10); + if r != 0 { + break; + } + m = q; + s -= 1; + } + let fits = m + .to_i128() + .filter(|v| v.unsigned_abs() < MANTISSA_LIMIT as u128); + match fits { + Some(v) if s <= MAX_SCALE => Interim::from(padded(v, s, nominal.min(MAX_SCALE))), + _ => settle_wide(BigQ::new(m.to_bigint(), BigInt::from(10).pow(s))), + } +} + +/// The stored scale of clause 4: pad toward `target` while the mantissa fits 96 bits. +pub(crate) fn padded(mantissa: i128, scale: u32, target: u32) -> Num { + let (mut m, mut s) = (mantissa, scale); + while s < target && (m * 10).unsigned_abs() < MANTISSA_LIMIT as u128 { + m *= 10; + s += 1; + } + decimal(m, s).expect("a padded mantissa still fits") +} + +/// The exponent of an exact integer-valued operand. +fn integer_value(b: Operand) -> Option { + match b { + Operand::Dec(d) if d.fract().is_zero() => { + Some(BigInt::from(d.mantissa() / pow10(d.scale()))) + } + Operand::Small(n, 1) => Some(BigInt::from(n)), + Operand::Big(q) if q.den.is_one() => Some(q.num.clone()), + _ => None, + } +} + +fn power(a: Operand, b: Operand) -> Result { + match (a.is_float(), integer_value(b)) { + (false, Some(n)) => exact_power(a, &n), + _ => Ok(inexact(a.float().powf(b.float()))), + } +} + +/// Clause 5's `POW`: the size is checked before anything is computed. +fn exact_power(a: Operand, n: &BigInt) -> Result { + let q = a.big().expect("an exact base"); + let base_scale = match a { + Operand::Dec(d) => d.scale(), + _ => 0, + }; + let k = n.magnitude(); + let target = match n.is_negative() { + true => 0, + false => k + .to_u32() + .and_then(|k| k.checked_mul(base_scale)) + .map_or(MAX_SCALE, |s| s.min(MAX_SCALE)), + }; + let unit = q.den.is_one() && q.num.magnitude().is_one(); + match (q.num.is_zero(), n.is_negative(), unit) { + (_, _, _) if n.is_zero() => Ok(Interim::from(padded(1, 0, 0))), + (true, true, _) => Err(DivisionByZero), + (true, false, _) => Ok(Interim::from(padded(0, 0, target))), + (false, _, true) => { + let negative = q.num.is_negative() && k.is_odd(); + Ok(Interim::from(padded( + if negative { -1 } else { 1 }, + 0, + target, + ))) + } + (false, negative, false) => { + let lower = |x: &BigInt| BigUint::from(x.bits().saturating_sub(1)) * k; + let too_big = lower(&q.num).max(lower(&q.den)) > BigUint::from(CEILING_BITS); + let odd_negative = q.num.is_negative() && k.is_odd(); + match (too_big, k.to_u32()) { + (false, Some(k)) => Ok(exact_integer_power(&q, k, negative, target)), + _ => Ok(inexact(signed_powf(a.float(), n, odd_negative))), + } + } + } +} + +/// `powf` on the magnitude, the sign taken from the exact exponent's parity. +fn signed_powf(base: f64, n: &BigInt, odd_negative: bool) -> f64 { + let magnitude = base.abs().powf(n.to_f64().unwrap_or(f64::INFINITY)); + match odd_negative { + true => -magnitude, + false => magnitude, + } +} + +fn exact_integer_power(q: &BigQ, k: u32, negative: bool, target: u32) -> Interim { + let (num, den) = (q.num.pow(k), q.den.pow(k)); + let value = match negative { + true => BigQ::new(den, num), + false => BigQ { num, den }, + }; + let settled = settle_wide(value); + match (&settled.0, negative) { + (Inner::Stored(n), false) => match n.kind() { + NumRef::Exact(d) => Interim::from(padded(d.mantissa(), d.scale(), target)), + _ => settled, + }, + _ => settled, + } +} diff --git a/crates/numtower/src/cmp.rs b/crates/numtower/src/cmp.rs new file mode 100644 index 00000000..81cac9b7 --- /dev/null +++ b/crates/numtower/src/cmp.rs @@ -0,0 +1,107 @@ +//! Exact comparison at every pairing of arms (clause 6), and the `MIN`/`MAX` +//! tie-breaks (clause 8). Why: `context/code/numtower.md`, "Comparison". + +use std::cmp::Ordering; + +use crate::fixed::I256; +use crate::num::{Num, NumRef}; +use crate::ratio::BigQ; + +impl Num { + /// The order of the two values, exactly; `None` only when either is NaN. + pub fn cmp_exact(&self, o: &Num) -> Option { + match (self.exact_parts(), o.exact_parts()) { + (Some(a), Some(b)) => Some(cross(a, b)), + (Some(_), None) => mixed(self, o.to_f64()), + (None, Some(_)) => mixed(o, self.to_f64()).map(Ordering::reverse), + (None, None) => self.to_f64().partial_cmp(&o.to_f64()), + } + } + + /// The smaller, NaN skipped; ties keep the most information (clause 8, question 11). + pub fn min_of<'a>(&'a self, o: &'a Num) -> &'a Num { + pick(self, o, Ordering::Less) + } + + pub fn max_of<'a>(&'a self, o: &'a Num) -> &'a Num { + pick(self, o, Ordering::Greater) + } +} + +impl PartialEq for Num { + fn eq(&self, o: &Num) -> bool { + self.cmp_exact(o) == Some(Ordering::Equal) + } +} + +impl PartialOrd for Num { + fn partial_cmp(&self, o: &Num) -> Option { + self.cmp_exact(o) + } +} + +/// `a.0/a.1` against `b.0/b.1`, denominators positive, in 256 bits. +fn cross(a: (i128, i128), b: (i128, i128)) -> Ordering { + I256::mul_i128(a.0, b.1).cmp(&I256::mul_i128(b.0, a.1)) +} + +/// An exact value against an f64: the rounded value decides unless it ties. +fn mixed(exact: &Num, x: f64) -> Option { + match (x.is_nan(), x.is_infinite()) { + (true, _) => None, + (false, true) => Some(if x > 0.0 { + Ordering::Less + } else { + Ordering::Greater + }), + (false, false) => match exact.to_f64().partial_cmp(&x)? { + Ordering::Equal => { + let q = exact.exact_q().expect("an exact operand"); + Some(q.cmp_value(&BigQ::from_f64(x))) + } + decided => Some(decided), + }, + } +} + +/// A wide exact value against an f64, with no fast path: the slow, total form of `mixed`. +pub(crate) fn exact_against_f64(q: &BigQ, x: f64) -> Option { + match (x.is_nan(), x.is_infinite(), x > 0.0) { + (true, _, _) => None, + (false, true, true) => Some(Ordering::Less), + (false, true, false) => Some(Ordering::Greater), + (false, false, _) => Some(q.cmp_value(&BigQ::from_f64(x))), + } +} + +fn pick<'a>(a: &'a Num, b: &'a Num, want: Ordering) -> &'a Num { + match (a.cmp_exact(b), is_nan(a), is_nan(b)) { + (_, true, _) => b, + (_, false, true) => a, + (Some(o), _, _) if o == want => a, + (Some(Ordering::Equal), _, _) => tie(a, b, want), + _ => b, + } +} + +fn is_nan(n: &Num) -> bool { + matches!(n.kind(), NumRef::Inexact(x) if x.is_nan()) +} + +/// Equal values: exact over inexact, the larger scale, and IEEE's signed zero. +fn tie<'a>(a: &'a Num, b: &'a Num, want: Ordering) -> &'a Num { + match (a.kind(), b.kind()) { + (NumRef::Inexact(_), NumRef::Exact(_) | NumRef::Fraction(_)) => b, + (NumRef::Exact(_) | NumRef::Fraction(_), NumRef::Inexact(_)) => a, + (NumRef::Exact(x), NumRef::Exact(y)) if y.scale() > x.scale() => b, + (NumRef::Inexact(x), NumRef::Inexact(y)) => { + let negative_wanted = want == Ordering::Less; + match y.is_sign_negative() == negative_wanted && x.is_sign_negative() != negative_wanted + { + true => b, + false => a, + } + } + _ => a, + } +} diff --git a/crates/numtower/src/fixed.rs b/crates/numtower/src/fixed.rs new file mode 100644 index 00000000..ba551181 --- /dev/null +++ b/crates/numtower/src/fixed.rs @@ -0,0 +1,234 @@ +//! Fixed-width two's-complement integers for the hot paths: no allocation. +//! Why they are hand-rolled: `context/code/numtower.md`, "The wide integers". + +use std::cmp::Ordering; + +use num_bigint::{BigInt, Sign}; + +/// A signed integer of `64 * N` bits, little-endian limbs, two's complement. +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub(crate) struct Fixed { + limbs: [u64; N], +} + +/// The 256-bit integer clauses 4–6 use, and the 320-bit decimal partition of clause 8. +pub(crate) type I256 = Fixed<4>; +pub(crate) type I320 = Fixed<5>; + +impl Fixed { + pub(crate) const ZERO: Self = Fixed { limbs: [0; N] }; + + pub(crate) fn from_i128(x: i128) -> Self { + let fill = if x < 0 { u64::MAX } else { 0 }; + let mut limbs = [fill; N]; + limbs[0] = x as u64; + limbs[1] = (x >> 64) as u64; + Fixed { limbs } + } + + pub(crate) fn is_negative(&self) -> bool { + self.limbs[N - 1] >> 63 == 1 + } + + pub(crate) fn is_zero(&self) -> bool { + self.limbs.iter().all(|&l| l == 0) + } + + fn wrapping_neg(&self) -> Self { + let mut out = [0u64; N]; + let mut carry = 1u64; + for (o, &l) in out.iter_mut().zip(&self.limbs) { + let (v, c) = (!l).overflowing_add(carry); + *o = v; + carry = c as u64; + } + Fixed { limbs: out } + } + + fn wrapping_add(&self, other: &Self) -> Self { + let mut out = [0u64; N]; + let mut carry = false; + for (o, (a, b)) in out.iter_mut().zip(self.limbs.iter().zip(&other.limbs)) { + let (s, c1) = a.overflowing_add(*b); + let (s, c2) = s.overflowing_add(carry as u64); + *o = s; + carry = c1 || c2; + } + Fixed { limbs: out } + } + + pub(crate) fn checked_add(&self, other: &Self) -> Option { + let sum = self.wrapping_add(other); + let same = self.is_negative() == other.is_negative(); + (!same || sum.is_negative() == self.is_negative()).then_some(sum) + } + + pub(crate) fn checked_neg(&self) -> Option { + let neg = self.wrapping_neg(); + (self.is_zero() || neg.is_negative() != self.is_negative()).then_some(neg) + } + + /// The magnitude as unsigned limbs, with the sign split out. + fn magnitude(&self) -> (bool, [u64; N]) { + match self.is_negative() { + true => (true, self.wrapping_neg().limbs), + false => (false, self.limbs), + } + } + + fn from_magnitude(negative: bool, limbs: [u64; N]) -> Option { + let unsigned = Fixed { limbs }; + match (negative, unsigned.is_negative()) { + (false, false) => Some(unsigned), + (true, _) => { + let neg = unsigned.wrapping_neg(); + (neg.is_negative() || neg.is_zero()).then_some(neg) + } + (false, true) => None, + } + } + + pub(crate) fn checked_mul_u128(&self, k: u128) -> Option { + let (negative, mag) = self.magnitude(); + let mut out = [0u64; N]; + for (j, kk) in [k as u64, (k >> 64) as u64].into_iter().enumerate() { + if kk == 0 { + continue; + } + let mut carry = 0u128; + for (i, &limb) in mag.iter().enumerate() { + let idx = i + j; + let prod = limb as u128 * kk as u128 + carry; + match idx < N { + true => { + let s = out[idx] as u128 + (prod as u64) as u128; + out[idx] = s as u64; + carry = (prod >> 64) + (s >> 64); + } + false if prod != 0 => return None, + false => carry = 0, + } + } + if carry != 0 { + return None; + } + } + Self::from_magnitude(negative, out) + } + + /// Divides the magnitude by `k`, returning the quotient (sign kept) and the remainder. + pub(crate) fn div_rem_u64(&self, k: u64) -> (Self, u64) { + let (negative, mag) = self.magnitude(); + let mut out = [0u64; N]; + let mut rem = 0u128; + for i in (0..N).rev() { + let cur = (rem << 64) | mag[i] as u128; + out[i] = (cur / k as u128) as u64; + rem = cur % k as u128; + } + let q = Self::from_magnitude(negative, out).unwrap_or(Self::ZERO); + (q, rem as u64) + } + + pub(crate) fn to_i128(self) -> Option { + let fill = if self.is_negative() { u64::MAX } else { 0 }; + let hi_ok = self.limbs[2..].iter().all(|&l| l == fill); + let v = (self.limbs[0] as u128 | ((self.limbs[1] as u128) << 64)) as i128; + (hi_ok && (v < 0) == self.is_negative()).then_some(v) + } + + pub(crate) fn to_bigint(self) -> BigInt { + let (negative, mag) = self.magnitude(); + let digits: Vec = mag + .iter() + .flat_map(|&l| [l as u32, (l >> 32) as u32]) + .collect(); + let sign = if negative { Sign::Minus } else { Sign::Plus }; + BigInt::from_slice(sign, &digits) + } + + /// The product of two `i128`s; needs `N >= 4`. + pub(crate) fn mul_i128(a: i128, b: i128) -> Self { + let negative = (a < 0) != (b < 0); + let (x, y) = (a.unsigned_abs(), b.unsigned_abs()); + let (x0, x1) = (x as u64 as u128, x >> 64); + let (y0, y1) = (y as u64 as u128, y >> 64); + let p00 = x0 * y0; + let p01 = x0 * y1; + let p10 = x1 * y0; + let p11 = x1 * y1; + let mid = (p00 >> 64) + (p01 as u64 as u128) + (p10 as u64 as u128); + let mut limbs = [0u64; N]; + limbs[0] = p00 as u64; + limbs[1] = mid as u64; + let hi = (mid >> 64) + (p01 >> 64) + (p10 >> 64) + p11; + limbs[2] = hi as u64; + limbs[3] = (hi >> 64) as u64; + let unsigned = Fixed { limbs }; + match negative { + true => unsigned.wrapping_neg(), + false => unsigned, + } + } +} + +impl Default for Fixed { + fn default() -> Self { + Self::ZERO + } +} + +impl Ord for Fixed { + fn cmp(&self, other: &Self) -> Ordering { + match (self.is_negative(), other.is_negative()) { + (true, false) => Ordering::Less, + (false, true) => Ordering::Greater, + _ => self.limbs.iter().rev().cmp(other.limbs.iter().rev()), + } + } +} + +impl PartialOrd for Fixed { + fn partial_cmp(&self, other: &Self) -> Option { + Some(self.cmp(other)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn the_256_bit_product_of_two_i128s_matches_bigint() { + let cases = [ + (i128::MAX, i128::MAX), + (i128::MAX, -i128::MAX), + (-1, 1), + (0, i128::MIN), + (123_456_789_012_345_678_901_234, -98_765_432_109_876_543_210), + ]; + for (a, b) in cases { + let want = BigInt::from(a) * BigInt::from(b); + assert_eq!(I256::mul_i128(a, b).to_bigint(), want, "{a} * {b}"); + } + } + + #[test] + fn checked_arithmetic_refuses_what_overflows_and_keeps_what_fits() { + let max = I256::from_i128(i128::MAX) + .checked_mul_u128(u128::MAX) + .unwrap(); + assert_eq!( + max.to_bigint(), + BigInt::from(i128::MAX) * BigInt::from(u128::MAX) + ); + assert!(max.checked_mul_u128(u128::MAX).is_none()); + let neg = max.checked_neg().unwrap(); + assert_eq!(neg.to_bigint(), -max.to_bigint()); + let (q, r) = neg.div_rem_u64(10); + assert_eq!(q.to_bigint() * 10 - BigInt::from(r), neg.to_bigint()); + assert_eq!(I256::from_i128(-5).to_i128(), Some(-5)); + assert_eq!(max.to_i128(), None); + assert!(I256::from_i128(-1) < I256::ZERO); + } +} diff --git a/crates/numtower/src/lang.rs b/crates/numtower/src/lang.rs new file mode 100644 index 00000000..b5d00960 --- /dev/null +++ b/crates/numtower/src/lang.rs @@ -0,0 +1,192 @@ +//! Clause 3's one number language: the gate, the exact parse, the forced `n/d` arm, +//! and the canonical text. Why each rule: `context/code/numtower.md`, "The number language". + +use std::fmt; + +use num_bigint::BigInt; +use num_traits::{One, Zero}; + +use crate::arith::padded; +use crate::num::{MAX_SCALE, Num, NumRef}; +use crate::ratio::BigQ; + +/// Whether the type of the text is being guessed or was forced (question 12). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Typing { + Inferred, + Forced, +} + +/// The gate: the strings `f64::from_str` accepts and reads as a finite value. +pub fn is_number(text: &str, typing: Typing) -> bool { + gate(text) || (typing == Typing::Forced && forced_fraction(text).is_some()) +} + +/// A number if the text is one, read exactly; `None` means the text is text. +pub fn parse_number(text: &str, typing: Typing) -> Option { + match (gate(text), typing) { + (true, _) => Lexeme::read(text).map(|lx| lx.value(text)), + (false, Typing::Forced) => forced_fraction(text), + (false, Typing::Inferred) => None, + } +} + +/// A forced `n/d`; one that rounds to infinity is text, as a finite lexeme never reads as one. +fn forced_fraction(text: &str) -> Option { + let n = Num::settle(&fraction_lexeme(text)?); + n.to_f64().is_finite().then_some(n) +} + +fn gate(text: &str) -> bool { + text.parse::().is_ok_and(f64::is_finite) +} + +/// `[+-]digits/digits` with a nonzero denominator, as an exact value. +fn fraction_lexeme(text: &str) -> Option { + let (n, d) = text.split_once('/')?; + let digits = |s: &str| !s.is_empty() && s.bytes().all(|b| b.is_ascii_digit()); + let unsigned = n.strip_prefix(['+', '-']).unwrap_or(n); + let ok = digits(unsigned) && digits(d); + let num: BigInt = ok.then(|| n.trim_start_matches('+').parse().ok())??; + let den: BigInt = d.parse().ok()?; + (!den.is_zero()).then(|| BigQ::new(num, den)) +} + +/// A lexeme of the gate's language: sign, significant digits, and the power of ten. +#[derive(Debug, PartialEq, Eq)] +pub(crate) struct Lexeme { + pub(crate) negative: bool, + pub(crate) digits: String, + pub(crate) exp10: i64, +} + +/// Exponents beyond this cannot reach either exact arm, and saturate. +const EXP_SATURATE: i64 = 1 << 40; + +impl Lexeme { + /// `[+-]? (d+ (. d*)? | . d+) ([eE] [+-]? d+)?`, and nothing else. + pub(crate) fn read(text: &str) -> Option { + let (negative, rest) = match text.as_bytes().first()? { + b'-' => (true, &text[1..]), + b'+' => (false, &text[1..]), + _ => (false, text), + }; + let (mantissa, exponent) = match rest.find(['e', 'E']) { + Some(i) => (&rest[..i], Some(&rest[i + 1..])), + None => (rest, None), + }; + let (int, frac) = mantissa.split_once('.').unwrap_or((mantissa, "")); + let all_digits = |s: &str| s.bytes().all(|b| b.is_ascii_digit()); + let mantissa_ok = + all_digits(int) && all_digits(frac) && !(int.is_empty() && frac.is_empty()); + let exp = match exponent { + Some(e) => read_exponent(e)?, + None => 0, + }; + let digits = format!("{int}{frac}"); + mantissa_ok.then(|| Lexeme { + negative, + digits, + exp10: exp.saturating_sub(frac.len() as i64), + }) + } + + /// Clause 3's order: a decimal at the lexeme's scale, then a fraction, then the nearest f64. + fn value(&self, text: &str) -> Num { + let significant = self.digits.trim_start_matches('0'); + let stripped = significant.trim_end_matches('0'); + let implied_scale = (-self.exp10).clamp(0, MAX_SCALE as i64) as u32; + match stripped.is_empty() { + true => padded(0, 0, implied_scale), + false => { + let shift = self + .exp10 + .saturating_add((significant.len() - stripped.len()) as i64); + self.exact(stripped, shift, implied_scale) + .unwrap_or_else(|| Num::inexact(text.parse().unwrap_or(f64::NAN))) + } + } + } + + fn exact(&self, stripped: &str, shift: i64, implied_scale: u32) -> Option { + let leading = stripped.len() as i64 + shift; + let reachable = stripped.len() <= 180 && (-130..=40).contains(&leading) && shift >= -400; + let q = reachable.then(|| self.rational(stripped, shift))?; + Num::fit(&q).map(|n| match n.kind() { + NumRef::Exact(d) => padded(d.mantissa(), d.scale(), implied_scale), + _ => n, + }) + } + + pub(crate) fn rational(&self, stripped: &str, shift: i64) -> BigQ { + let m: BigInt = stripped.parse().unwrap_or_default(); + let m = if self.negative { -m } else { m }; + let ten = BigInt::from(10); + match shift >= 0 { + true => BigQ::new(m * ten.pow(shift as u32), BigInt::one()), + false => BigQ::new(m, ten.pow((-shift) as u32)), + } + } + + /// The exact value, whatever its size: for display rounding. + pub(crate) fn exact_value(&self) -> BigQ { + let stripped = self.digits.trim_start_matches('0'); + match stripped.is_empty() { + true => BigQ::new(BigInt::zero(), BigInt::one()), + false => self.rational(stripped, self.exp10), + } + } +} + +fn read_exponent(e: &str) -> Option { + let (negative, digits) = match e.as_bytes().first()? { + b'-' => (true, &e[1..]), + b'+' => (false, &e[1..]), + _ => (false, e), + }; + let ok = !digits.is_empty() && digits.bytes().all(|b| b.is_ascii_digit()); + let significant = digits.trim_start_matches('0'); + let magnitude = match significant.len() > 12 { + true => EXP_SATURATE, + false => significant.parse::().unwrap_or(0).min(EXP_SATURATE), + }; + ok.then_some(if negative { -magnitude } else { magnitude }) +} + +/// The shortest text that reads back as `x`; `None` for a non-finite value. +pub fn shortest_round_trip(x: f64) -> Option { + x.is_finite().then(|| format!("{x:?}")) +} + +const INEXACT_MARK: &str = "#i"; + +impl Num { + /// Reads every spelling `Display` writes: a decimal, `n/d`, `#i…`, `inf`, `-inf`, `nan`. + pub fn from_canonical(text: &str) -> Option { + match text { + "inf" => Some(Num::inexact(f64::INFINITY)), + "-inf" => Some(Num::inexact(f64::NEG_INFINITY)), + "nan" => Some(Num::inexact(f64::NAN)), + _ => match text.strip_prefix(INEXACT_MARK) { + Some(rest) if gate(rest) => rest.parse().ok().map(Num::inexact), + Some(_) => None, + None => parse_number(text, Typing::Forced).filter(Num::is_exact), + }, + } + } +} + +/// The canonical text: a decimal with its scale, `n/d`, `#i` and the shortest spelling, or a token. +impl fmt::Display for Num { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self.kind() { + NumRef::Exact(d) => write!(f, "{d}"), + NumRef::Fraction(r) => write!(f, "{}/{}", r.numer(), r.denom()), + NumRef::Inexact(x) if x.is_nan() => f.write_str("nan"), + NumRef::Inexact(x) if x.is_infinite() => { + f.write_str(if x > 0.0 { "inf" } else { "-inf" }) + } + NumRef::Inexact(x) => write!(f, "{INEXACT_MARK}{x:?}"), + } + } +} diff --git a/crates/numtower/src/lib.rs b/crates/numtower/src/lib.rs new file mode 100644 index 00000000..458d3e9f --- /dev/null +++ b/crates/numtower/src/lib.rs @@ -0,0 +1,20 @@ +//! A three-level numeric tower: exact decimal, exact bounded fraction, inexact f64. +//! The decisions are ADR 005; the reasoning is `context/code/numtower.md`. + +mod arith; +mod cmp; +mod fixed; +mod lang; +mod num; +mod ratio; +mod round; +mod sum; + +pub use arith::{CEILING_BITS, DivisionByZero, Interim}; +pub use lang::{Typing, is_number, parse_number, shortest_round_trip}; +pub use num::{Num, NumRef}; +pub use round::Rounding; +pub use sum::{ExactSum, K, NumSum}; + +pub use num_rational::Ratio; +pub use rust_decimal::Decimal; diff --git a/crates/numtower/src/num.rs b/crates/numtower/src/num.rs new file mode 100644 index 00000000..1f0448f1 --- /dev/null +++ b/crates/numtower/src/num.rs @@ -0,0 +1,198 @@ +//! The stored number: an exact decimal, an exact bounded fraction, or an f64. +//! The arms are private so a value is always normalized: `context/code/numtower.md`. + +use num_bigint::BigInt; +use num_integer::Integer; +use num_rational::Ratio; +use num_traits::{One, ToPrimitive}; +use rust_decimal::Decimal; + +use crate::ratio::{BigQ, ratio_to_f64}; + +/// The most significant digits a `Decimal` holds, and its mantissa bound. +pub(crate) const MAX_SCALE: u32 = 28; +pub(crate) const MANTISSA_LIMIT: i128 = 1 << 96; + +/// A number of the tower (ADR 005 clause 2), always in its one normal form. +#[derive(Clone, Debug)] +pub struct Num(Repr); + +#[derive(Clone, Debug)] +enum Repr { + Exact(Decimal), + Fraction(Box>), + Inexact(f64), +} + +/// A borrowed view of a `Num`'s arm, for exhaustive matching by callers. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum NumRef<'a> { + Exact(&'a Decimal), + Fraction(&'a Ratio), + Inexact(f64), +} + +impl Num { + /// An exact decimal. Zero loses its sign: the sign of zero is an f64 artefact. + pub fn exact(d: Decimal) -> Num { + let mut d = d; + if d.is_zero() { + d.set_sign_positive(true); + } + Num(Repr::Exact(d)) + } + + pub fn inexact(x: f64) -> Num { + Num(Repr::Inexact(x)) + } + + /// The exact value `num / den`, normalized; `None` when `den` is zero. + pub fn ratio(num: i128, den: i128) -> Option { + (den != 0).then(|| Num::settle(&BigQ::from_i128(num, den))) + } + + fn fraction(num: i128, den: i128) -> Num { + Num(Repr::Fraction(Box::new(Ratio::new_raw(num, den)))) + } + + pub fn kind(&self) -> NumRef<'_> { + match &self.0 { + Repr::Exact(d) => NumRef::Exact(d), + Repr::Fraction(r) => NumRef::Fraction(r), + Repr::Inexact(x) => NumRef::Inexact(*x), + } + } + + pub fn is_exact(&self) -> bool { + !matches!(self.0, Repr::Inexact(_)) + } + + /// Exactly zero, at either level; the one zero predicate of clause 6. + pub fn is_zero(&self) -> bool { + match &self.0 { + Repr::Exact(d) => d.is_zero(), + Repr::Fraction(_) => false, + Repr::Inexact(x) => *x == 0.0, + } + } + + /// The correctly rounded f64 (clause 4): a decimal through its text, a fraction exactly. + pub fn to_f64(&self) -> f64 { + match &self.0 { + Repr::Exact(d) => decimal_to_f64(d), + Repr::Fraction(r) => ratio_to_f64(&BigInt::from(*r.numer()), &BigInt::from(*r.denom())), + Repr::Inexact(x) => *x, + } + } + + /// The exact value as a rational; `None` for an inexact number. + pub(crate) fn exact_parts(&self) -> Option<(i128, i128)> { + match &self.0 { + Repr::Exact(d) => Some((d.mantissa(), pow10(d.scale()))), + Repr::Fraction(r) => Some((*r.numer(), *r.denom())), + Repr::Inexact(_) => None, + } + } + + pub(crate) fn exact_q(&self) -> Option { + self.exact_parts().map(|(n, d)| BigQ::from_i128(n, d)) + } + + /// The stored bound: a decimal if it fits, else an `i128` fraction, else the nearest f64. + pub(crate) fn settle(q: &BigQ) -> Num { + Num::fit(q).unwrap_or_else(|| Num::inexact(q.to_f64())) + } + + /// `q` in an exact arm, at its minimal scale, if it fits one. + pub(crate) fn fit(q: &BigQ) -> Option { + decimal_of(q).or_else(|| { + let num = q.num.to_i128()?; + let den = q.den.to_i128()?; + Some(Num::fraction(num, den)) + }) + } +} + +/// `num / den` normalized without allocating; `None` hands the case to the bigint path. +pub(crate) fn small_ratio(num: i128, den: i128) -> Option { + let representable = num != i128::MIN && den != i128::MIN && den != 0; + representable.then(|| { + let g = num.gcd(&den); + let (n, d) = match den < 0 { + true => (-num / g, -den / g), + false => (num / g, den / g), + }; + let p = d.trailing_zeros(); + let mut rest = d >> p; + let mut q5 = 0; + while rest % 5 == 0 { + rest /= 5; + q5 += 1; + } + let k = p.max(q5); + let as_decimal = (rest == 1 && k <= MAX_SCALE) + .then(|| 2i128.pow(k - p) * 5i128.pow(k - q5)) + .and_then(|pad| n.checked_mul(pad)) + .and_then(|m| decimal(m, k)); + as_decimal.unwrap_or_else(|| Num::fraction(n, d)) + }) +} + +pub(crate) fn pow10(k: u32) -> i128 { + 10i128.pow(k) +} + +pub(crate) fn decimal_to_f64(d: &Decimal) -> f64 { + d.to_string().parse().unwrap_or(f64::NAN) +} + +/// A decimal from a mantissa and scale, if the mantissa fits 96 bits. +pub(crate) fn decimal(mantissa: i128, scale: u32) -> Option { + let fits = mantissa.unsigned_abs() < MANTISSA_LIMIT as u128 && scale <= MAX_SCALE; + fits.then(|| Num::exact(Decimal::from_i128_with_scale(mantissa, scale))) +} + +/// The direct test of clause 5: a reduced `n / (2^p 5^q)` is a decimal iff it fits. +fn decimal_of(q: &BigQ) -> Option { + let p = q.den.trailing_zeros().unwrap_or(0); + let mut rest = &q.den >> p as usize; + let five = BigInt::from(5); + let mut q5 = 0u64; + while rest.is_multiple_of(&five) && q5 <= MAX_SCALE as u64 { + rest /= &five; + q5 += 1; + } + let k = p.max(q5); + if !rest.is_one() || k > MAX_SCALE as u64 { + return None; + } + let pad = (BigInt::one() << (k - p) as usize) * five.pow((k - q5) as u32); + let mantissa = (&q.num * pad).to_i128()?; + decimal(mantissa, k as u32) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn a_num_is_24_bytes() { + assert_eq!(std::mem::size_of::(), 24); + } + + #[test] + fn a_fraction_that_is_a_fitting_decimal_is_stored_as_the_decimal() { + let half = Num::ratio(1, 2).unwrap(); + assert_eq!(half.kind(), NumRef::Exact(&Decimal::new(5, 1))); + let third = Num::ratio(-10, 3).unwrap(); + assert_eq!(third.kind(), NumRef::Fraction(&Ratio::new_raw(-10, 3))); + let tiny = Num::ratio(1, 1 << 40).unwrap(); + assert!(matches!(tiny.kind(), NumRef::Fraction(_)), "{tiny:?}"); + let mut minus_zero = Decimal::new(0, 2); + minus_zero.set_sign_negative(true); + assert!(minus_zero.is_sign_negative()); + let zero = Num::exact(minus_zero); + assert!(!zero.to_f64().is_sign_negative()); + assert!(Num::ratio(1, 0).is_none()); + } +} diff --git a/crates/numtower/src/ratio.rs b/crates/numtower/src/ratio.rs new file mode 100644 index 00000000..0c1a2d5c --- /dev/null +++ b/crates/numtower/src/ratio.rs @@ -0,0 +1,172 @@ +//! Exact rationals in bigints: the widened intermediate of question 13 (b), and +//! the one correctly rounded conversion of a rational to f64 (clauses 4 and 5). + +use std::cmp::Ordering; + +use num_bigint::{BigInt, BigUint, Sign}; +use num_integer::Integer; +use num_traits::{One, Signed, ToPrimitive, Zero}; + +/// A reduced fraction with a positive denominator. +#[derive(Clone, Debug, PartialEq, Eq)] +pub(crate) struct BigQ { + pub(crate) num: BigInt, + pub(crate) den: BigInt, +} + +impl BigQ { + pub(crate) fn new(num: BigInt, den: BigInt) -> BigQ { + let (num, den) = match den.is_negative() { + true => (-num, -den), + false => (num, den), + }; + let g = num.gcd(&den); + match g.is_zero() || g.is_one() { + true => BigQ { num, den }, + false => BigQ { + num: num / &g, + den: den / &g, + }, + } + } + + pub(crate) fn from_i128(num: i128, den: i128) -> BigQ { + BigQ::new(BigInt::from(num), BigInt::from(den)) + } + + pub(crate) fn add(&self, o: &BigQ) -> BigQ { + BigQ::new(&self.num * &o.den + &o.num * &self.den, &self.den * &o.den) + } + + pub(crate) fn neg(&self) -> BigQ { + BigQ { + num: -&self.num, + den: self.den.clone(), + } + } + + pub(crate) fn mul(&self, o: &BigQ) -> BigQ { + BigQ::new(&self.num * &o.num, &self.den * &o.den) + } + + /// `None` when `o` is zero. + pub(crate) fn div(&self, o: &BigQ) -> Option { + (!o.num.is_zero()).then(|| BigQ::new(&self.num * &o.den, &self.den * &o.num)) + } + + pub(crate) fn cmp_value(&self, o: &BigQ) -> Ordering { + (&self.num * &o.den).cmp(&(&o.num * &self.den)) + } + + pub(crate) fn bits(&self) -> u64 { + self.num.bits().max(self.den.bits()) + } + + pub(crate) fn to_f64(&self) -> f64 { + ratio_to_f64(&self.num, &self.den) + } + + /// The exact value of a finite f64. + pub(crate) fn from_f64(x: f64) -> BigQ { + let bits = x.to_bits(); + let biased = ((bits >> 52) & 0x7ff) as i64; + let fraction = bits & ((1u64 << 52) - 1); + let (mant, exp) = match biased { + 0 => (fraction, -1074), + e => (fraction | (1 << 52), e - 1075), + }; + let signed = match x.is_sign_negative() { + true => -BigInt::from(mant), + false => BigInt::from(mant), + }; + match exp >= 0 { + true => BigQ::new(signed << exp as usize, BigInt::one()), + false => BigQ::new(signed, BigInt::one() << (-exp) as usize), + } + } +} + +/// `num / den` rounded to the nearest f64, ties to even, subnormals and overflow included. +pub(crate) fn ratio_to_f64(num: &BigInt, den: &BigInt) -> f64 { + let negative = (num.sign() == Sign::Minus) != (den.sign() == Sign::Minus); + let a = num.magnitude(); + let b = den.magnitude(); + let magnitude = match a.is_zero() { + true => 0.0, + false => magnitude_to_f64(a, b), + }; + match negative { + true => -magnitude, + false => magnitude, + } +} + +fn magnitude_to_f64(a: &BigUint, b: &BigUint) -> f64 { + let shift = 54 - (a.bits() as i64 - b.bits() as i64); + let (n, d) = match shift >= 0 { + true => (a << shift as usize, b.clone()), + false => (a.clone(), b << (-shift) as usize), + }; + let (q, r) = n.div_rem(&d); + let qbits = q.bits() as i64; + let drop = (qbits - 53).max(shift - 1074); + if drop > qbits + 1 { + return 0.0; + } + let kept = (&q >> drop as usize).to_u64().unwrap_or(0); + let round = q.bit(drop as u64 - 1); + let sticky = !r.is_zero() || q.trailing_zeros().is_some_and(|t| (t as i64) < drop - 1); + let up = round && (sticky || kept & 1 == 1); + ldexp(kept + up as u64, drop - shift) +} + +/// `m * 2^e` for `m <= 2^53`, exact unless it overflows to infinity. +fn ldexp(m: u64, e: i64) -> f64 { + let top = e + 63 - m.leading_zeros() as i64; + match (m == 0, top > 1023) { + (true, _) => 0.0, + (false, true) => f64::INFINITY, + (false, false) => { + let half = e / 2; + m as f64 * pow2(half) * pow2(e - half) + } + } +} + +fn pow2(e: i64) -> f64 { + match e { + e if e >= -1022 => f64::from_bits(((e + 1023) as u64) << 52), + e => f64::from_bits(1u64 << (e + 1074)), + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn rounding_to_f64_reaches_subnormals_ties_and_overflow() { + let one = BigInt::one(); + let min_sub = f64::from_bits(1); + assert_eq!(ratio_to_f64(&one, &(BigInt::one() << 1074)), min_sub); + assert_eq!(ratio_to_f64(&one, &(BigInt::one() << 1075)), 0.0); + let three = BigInt::from(3); + assert_eq!(ratio_to_f64(&three, &(BigInt::one() << 1076)), min_sub); + assert_eq!(ratio_to_f64(&one, &BigInt::from(3)), 1.0 / 3.0); + let tie = (BigInt::one() << 53) + 1; + assert_eq!(ratio_to_f64(&tie, &one), 9007199254740992.0); + let big = BigInt::one() << 1024; + assert_eq!(ratio_to_f64(&big, &one), f64::INFINITY); + assert_eq!(ratio_to_f64(&-big, &one), f64::NEG_INFINITY); + let max = BigQ::from_f64(f64::MAX); + assert_eq!(max.to_f64(), f64::MAX); + } + + #[test] + fn every_f64_converts_to_its_exact_value_and_back() { + for x in [0.1, -2.5e-310, 1e300, f64::MIN_POSITIVE, 123456.789, -7.0] { + assert_eq!(BigQ::from_f64(x).to_f64().to_bits(), x.to_bits(), "{x}"); + } + assert!(BigQ::from_f64(-0.0).num.is_zero()); + } +} diff --git a/crates/numtower/src/round.rs b/crates/numtower/src/round.rs new file mode 100644 index 00000000..c1a17633 --- /dev/null +++ b/crates/numtower/src/round.rs @@ -0,0 +1,107 @@ +//! Rounding any value to `places` for display, exactly (clause 12, "Format"). +//! An inexact value is rounded over its shortest round-trip spelling (improvise-felo.9). + +use num_bigint::BigInt; +use num_integer::Integer; +use num_traits::Signed; + +use crate::lang::{Lexeme, shortest_round_trip}; +use crate::num::Num; +use crate::ratio::BigQ; + +/// The rounding mode is the caller's policy, not the crate's. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Rounding { + HalfAwayFromZero, + HalfEven, +} + +impl Num { + /// Plain decimal text with exactly `places` fractional digits; `None` for a non-finite value. + pub fn round_to_places(&self, places: u8, mode: Rounding) -> Option { + let q = match self.exact_q() { + Some(q) => q, + None => Lexeme::read(&shortest_round_trip(self.to_f64())?)?.exact_value(), + }; + Some(render(round(&q, places, mode), places)) + } +} + +/// `q * 10^places`, rounded to an integer by `mode`. +fn round(q: &BigQ, places: u8, mode: Rounding) -> BigInt { + let scaled = q.num.abs() * BigInt::from(10).pow(places as u32); + let (whole, rest) = scaled.div_rem(&q.den); + let twice: BigInt = rest * 2; + let up = match twice.cmp(&q.den) { + std::cmp::Ordering::Greater => true, + std::cmp::Ordering::Equal => mode == Rounding::HalfAwayFromZero || whole.is_odd(), + std::cmp::Ordering::Less => false, + }; + let magnitude = whole + BigInt::from(up as u8); + match q.num.is_negative() { + true => -magnitude, + false => magnitude, + } +} + +fn render(units: BigInt, places: u8) -> String { + let digits = units.abs().to_string(); + let width = places as usize + 1; + let padded = format!("{digits:0>width$}"); + let (int, frac) = padded.split_at(padded.len() - places as usize); + let sign = if units.is_negative() { "-" } else { "" }; + match places { + 0 => format!("{sign}{int}"), + _ => format!("{sign}{int}.{frac}"), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::lang::{Typing, parse_number}; + + fn typed(s: &str) -> Num { + parse_number(s, Typing::Inferred).unwrap() + } + + #[test] + fn a_typed_value_rounds_as_it_was_typed() { + let at2 = |n: &Num| n.round_to_places(2, Rounding::HalfAwayFromZero).unwrap(); + for (text, want) in [ + ("1.005", "1.01"), + ("1.015", "1.02"), + ("0.285", "0.29"), + ("-2.675", "-2.68"), + ] { + assert_eq!(at2(&typed(text)), want, "{text}"); + assert_eq!( + at2(&Num::inexact(text.parse().unwrap())), + want, + "inexact {text}" + ); + } + assert_eq!(at2(&Num::ratio(10, 3).unwrap()), "3.33"); + assert_eq!(at2(&Num::ratio(2, 3).unwrap()), "0.67"); + assert_eq!(at2(&typed("-0.001")), "0.00"); + let even = typed("0.125") + .round_to_places(2, Rounding::HalfEven) + .unwrap(); + assert_eq!(even, "0.12"); + assert_eq!( + typed("12.5") + .round_to_places(0, Rounding::HalfEven) + .unwrap(), + "12" + ); + assert_eq!( + Num::inexact(f64::NAN).round_to_places(2, Rounding::HalfEven), + None + ); + let wide = Num::ratio(1, 3) + .unwrap() + .round_to_places(40, Rounding::HalfEven) + .unwrap(); + assert_eq!(wide, format!("0.{}", "3".repeat(40))); + } +} diff --git a/crates/numtower/src/sum.rs b/crates/numtower/src/sum.rs new file mode 100644 index 00000000..366edbb5 --- /dev/null +++ b/crates/numtower/src/sum.rs @@ -0,0 +1,396 @@ +//! Order-free, retractable sums (clause 8): ADR 009's f64 `ExactSum`, and `NumSum`, +//! which adds the exact partitions. The design: `context/code/numtower.md`, "Sums". + +use std::collections::BTreeMap; + +use num_bigint::{BigInt, BigUint}; +use num_integer::Integer; +use num_traits::{One, Zero}; + +use crate::fixed::{I256, I320}; +use crate::num::{MAX_SCALE, Num, NumRef, pow10}; +use crate::ratio::{BigQ, ratio_to_f64}; + +/// The most non-decimal denominator groups a sum reads exactly (question 16: K = 8). +pub const K: usize = 8; + +const LIMB_BITS: u32 = 32; +/// A finite f64 is `mant * 2^(shift - 1074)`, `shift <= 2045`; weights add 60 bits. +const LIMBS: usize = 70; +const CHUNK: u64 = 1 << 30; + +/// An exact sum of weighted f64s, rounded once when read (ADR 009). +#[derive(Clone, Debug)] +pub struct ExactSum { + limbs: [i64; LIMBS], + pending: u64, + nan: i128, + pos_inf: i128, + neg_inf: i128, +} + +/// Equal sums, not equal representations: carries propagated, `pending` ignored. +impl PartialEq for ExactSum { + fn eq(&self, o: &ExactSum) -> bool { + let (mut mine, mut theirs) = (self.limbs, o.limbs); + propagate(&mut mine); + propagate(&mut theirs); + let counts = |s: &ExactSum| (s.nan, s.pos_inf, s.neg_inf); + mine == theirs && counts(self) == counts(o) + } +} + +impl Eq for ExactSum {} + +impl Default for ExactSum { + fn default() -> Self { + ExactSum { + limbs: [0; LIMBS], + pending: 0, + nan: 0, + pos_inf: 0, + neg_inf: 0, + } + } +} + +impl ExactSum { + pub fn new() -> ExactSum { + ExactSum::default() + } + + /// Adds `x` with a Z-set `weight`; the weight is split into base-2^30 digits. + pub fn add(&mut self, x: f64, weight: i64) { + match (x.is_nan(), x.is_infinite(), x > 0.0) { + (true, _, _) => self.nan += weight as i128, + (false, true, true) => self.pos_inf += weight as i128, + (false, true, false) => self.neg_inf += weight as i128, + (false, false, _) => { + let magnitude = weight.unsigned_abs(); + for j in 0..3 { + let digit = (magnitude >> (30 * j)) & (CHUNK - 1); + if digit != 0 { + let signed = if weight < 0 { + -(digit as i64) + } else { + digit as i64 + }; + self.add_digit(x, signed, 30 * j); + } + } + } + } + } + + fn add_digit(&mut self, x: f64, weight: i64, extra: u32) { + if self.pending + weight.unsigned_abs() > CHUNK { + propagate(&mut self.limbs); + self.pending = 0; + } + 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 at = shift + extra; + let wide = u128::from(mant) << (at % LIMB_BITS); + let k = (at / 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 += weight.unsigned_abs(); + } + + /// Adds another accumulator's contents. + pub fn absorb(&mut self, o: &ExactSum) { + let mut theirs = o.limbs; + propagate(&mut theirs); + propagate(&mut self.limbs); + for (mine, t) in self.limbs.iter_mut().zip(theirs) { + *mine += t; + } + propagate(&mut self.limbs); + self.pending = 0; + self.nan += o.nan; + self.pos_inf += o.pos_inf; + self.neg_inf += o.neg_inf; + } + + /// The exact sum rounded to the nearest f64, ties to even. + pub fn value(&self) -> f64 { + match (self.nan != 0, self.pos_inf != 0, self.neg_inf != 0) { + (true, _, _) | (false, true, true) => f64::NAN, + (false, true, false) => f64::INFINITY, + (false, false, true) => f64::NEG_INFINITY, + (false, false, false) => ratio_to_f64(&self.total(), &(BigInt::one() << 1074)), + } + } + + /// The finite part as an integer in units of 2^-1074. + fn total(&self) -> BigInt { + let mut limbs = self.limbs; + propagate(&mut limbs); + let digits: Vec = limbs[..LIMBS - 1].iter().map(|&l| l as u32).collect(); + let low = BigInt::from(BigUint::from_slice(&digits)); + (BigInt::from(limbs[LIMBS - 1]) << (LIMB_BITS as usize * (LIMBS - 1))) + low + } +} + +/// 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 `SUM` over any `Num`s: exact when every term is exact and at most K groups remain. +#[derive(Clone, Debug, Default)] +pub struct NumSum { + decimals: I320, + groups: BTreeMap, + group_floats: ExactSum, + inexact: ExactSum, + /// Net weight per f64 bit pattern: a read is inexact while any term is present. + inexact_terms: BTreeMap, +} + +impl NumSum { + pub fn new() -> NumSum { + NumSum::default() + } + + pub fn insert(&mut self, x: &Num) { + self.add(x, 1); + } + + pub fn retract(&mut self, x: &Num) { + self.add(x, -1); + } + + /// Adds `x` with a Z-set `weight`. Exact and associative: the order cannot matter. + pub fn add(&mut self, x: &Num, weight: i64) { + match x.kind() { + NumRef::Inexact(f) => { + self.inexact.add(f, weight); + let net = self.inexact_terms.entry(f.to_bits()).or_insert(0); + *net += weight as i128; + if *net == 0 { + self.inexact_terms.remove(&f.to_bits()); + } + } + NumRef::Exact(d) => { + let unit = pow10(MAX_SCALE - d.scale()) as u128; + self.add_scaled(d.mantissa(), unit, weight); + } + NumRef::Fraction(r) => match at_scale_28(*r.denom()) { + Some(unit) => self.add_scaled(*r.numer(), unit, weight), + None => self.add_to_group(*r.numer(), *r.denom(), weight), + }, + } + } + + fn add_scaled(&mut self, n: i128, unit: u128, weight: i64) { + let magnitude = weight.unsigned_abs() as u128; + let scaled = match unit.checked_mul(magnitude) { + Some(both) => I320::from_i128(n).checked_mul_u128(both), + None => I320::from_i128(n) + .checked_mul_u128(unit) + .and_then(|t| t.checked_mul_u128(magnitude)), + }; + let term = scaled + .and_then(|t| if weight < 0 { t.checked_neg() } else { Some(t) }) + .expect("a weighted term is below 2^284"); + self.decimals = self + .decimals + .checked_add(&term) + .expect("320 bits holds 2^36 terms"); + } + + fn add_to_group(&mut self, n: i128, d: i128, weight: i64) { + let old = self.groups.get(&d).copied().unwrap_or(I256::ZERO); + let new = old + .checked_add(&I256::mul_i128(n, weight as i128)) + .expect("256 bits holds 2^64 terms"); + let value = |num: I256| ratio_to_f64(&num.to_bigint(), &BigInt::from(d)); + if !old.is_zero() { + self.group_floats.add(value(old), -1); + } + match new.is_zero() { + true => { + self.groups.remove(&d); + } + false => { + self.group_floats.add(value(new), 1); + self.groups.insert(d, new); + } + } + } + + /// The current non-decimal denominator groups. + pub fn groups(&self) -> usize { + self.groups.len() + } + + /// One read: exact over at most K groups, else the f64 accumulators (clause 8). + pub fn read(&self) -> Num { + match (!self.inexact_terms.is_empty(), self.groups.len() > K) { + (false, false) => Num::settle(&self.exact_total()), + (_, true) => { + let mut acc = self.inexact.clone(); + acc.absorb(&self.group_floats); + acc.add(self.decimal_total().to_f64(), 1); + Num::inexact(acc.value()) + } + (true, false) => { + let mut acc = self.inexact.clone(); + acc.add(self.exact_total().to_f64(), 1); + Num::inexact(acc.value()) + } + } + } + + fn decimal_total(&self) -> BigQ { + BigQ::new(self.decimals.to_bigint(), BigInt::from(pow10(MAX_SCALE))) + } + + fn exact_total(&self) -> BigQ { + let decimal = (!self.decimals.is_zero()) + .then(|| (self.decimals.to_bigint(), BigInt::from(pow10(MAX_SCALE)))); + let terms: Vec<(BigInt, BigInt)> = decimal + .into_iter() + .chain( + self.groups + .iter() + .map(|(d, n)| (n.to_bigint(), BigInt::from(*d))), + ) + .collect(); + let (num, den) = combine_over_lcm(&terms); + BigQ::new(num, den) + } +} + +/// `Σ n/d` over the lcm of the denominators, not yet reduced: the read reduces once. +pub(crate) fn combine_over_lcm(terms: &[(BigInt, BigInt)]) -> (BigInt, BigInt) { + let lcm = terms.iter().fold(BigInt::one(), |l, (_, d)| l.lcm(d)); + let num = terms + .iter() + .fold(BigInt::zero(), |s, (n, d)| s + n * (&lcm / d)); + (num, lcm) +} + +/// `10^28 / d` when `d = 2^p 5^q` with `max(p, q) <= 28`: the term joins the decimal partition. +fn at_scale_28(d: i128) -> Option { + let p = d.trailing_zeros(); + let mut rest = d >> p; + let mut q = 0; + while rest % 5 == 0 { + rest /= 5; + q += 1; + } + (rest == 1 && p.max(q) <= MAX_SCALE).then(|| (pow10(MAX_SCALE) / d) as u128) +} + +#[cfg(test)] +mod tests { + use super::*; + use num_rational::BigRational; + use proptest::prelude::*; + use proptest::test_runner::{Config, RngAlgorithm, TestRng, TestRunner}; + + type Combine = fn(&[(BigInt, BigInt)]) -> (BigInt, BigInt); + + /// Up to K groups of wide numerators whose denominators share factors. + fn groups() -> impl Strategy> { + let shared = prop::sample::select(vec![1i128, 3, 7 * 11, 1 << 20, 3i128.pow(40)]); + let group = + (any::(), any::(), shared, 1i128..(1 << 60)).prop_map(|(n, w, s, d)| { + ( + BigInt::from(n) * BigInt::from(w), + BigInt::from(s) * BigInt::from(d), + ) + }); + prop::collection::vec(group, 1..=K) + } + + fn check_grouped_read(combine: Combine, deterministic: bool) -> Result<(), String> { + let config = Config { + cases: 256, + failure_persistence: None, + ..Config::default() + }; + let mut runner = match deterministic { + true => { + TestRunner::new_with_rng(config, TestRng::deterministic_rng(RngAlgorithm::ChaCha)) + } + false => TestRunner::new(config), + }; + let result = runner.run(&groups(), |terms| { + let (num, den) = combine(&terms); + let lcm = terms.iter().fold(BigInt::one(), |l, (_, d)| l.lcm(d)); + prop_assert_eq!(&den, &lcm, "the read did not combine over the lcm"); + let oracle: BigRational = terms + .iter() + .map(|(n, d)| BigRational::new(n.clone(), d.clone())) + .sum(); + let read = BigQ::new(num, den); + let got = BigRational::new(read.num.clone(), read.den.clone()); + prop_assert_eq!(got, oracle, "the read is not the sum"); + prop_assert!(read.num.gcd(&read.den).is_one(), "the read is not reduced"); + Ok(()) + }); + result.map_err(|e| e.to_string()) + } + + #[test] + fn the_grouped_read_combines_over_the_lcm_and_reduces_once() { + check_grouped_read(combine_over_lcm, false).unwrap(); + } + + fn over_the_product(terms: &[(BigInt, BigInt)]) -> (BigInt, BigInt) { + let den = terms.iter().fold(BigInt::one(), |p, (_, d)| p * d); + let num = terms + .iter() + .fold(BigInt::zero(), |s, (n, d)| s + n * (&den / d)); + (num, den) + } + + #[test] + fn the_grouped_read_combines_over_the_lcm_and_reduces_once_catches_the_product() { + let err = check_grouped_read(over_the_product, true).expect_err("the rival passed"); + assert!( + err.contains("the read did not combine over the lcm"), + "{err}" + ); + } + + #[test] + fn a_fraction_with_a_decimal_denominator_joins_the_decimal_partition() { + let mut s = NumSum::new(); + s.insert(&Num::ratio(10i128.pow(35), 1).unwrap()); + s.insert(&Num::ratio(3, 1 << 20).unwrap()); + assert_eq!(s.groups(), 0); + s.insert(&Num::ratio(1, 1i128 << 40).unwrap()); + s.insert(&Num::ratio(1, 3).unwrap()); + assert_eq!( + s.groups(), + 2, + "2^40 has scale 40, past 28: a group of its own" + ); + s.retract(&Num::ratio(1, 3).unwrap()); + assert_eq!( + s.groups(), + 1, + "a group whose numerator returns to zero is dropped" + ); + } +} diff --git a/crates/numtower/tests/oracle/mod.rs b/crates/numtower/tests/oracle/mod.rs new file mode 100644 index 00000000..6218812b --- /dev/null +++ b/crates/numtower/tests/oracle/mod.rs @@ -0,0 +1,189 @@ +//! The `BigRational` oracle and the generators the property tests share. +//! Each property is a checker taking the implementation, so a rival can be run through it. + +#![allow(dead_code)] + +use std::cmp::Ordering; + +use num_bigint::BigInt; +use num_rational::BigRational; +use num_traits::{One, Signed, ToPrimitive, Zero}; +use numtower::{Decimal, Num, NumRef}; +use proptest::prelude::*; +use proptest::test_runner::{Config, RngAlgorithm, TestCaseError, TestRng, TestRunner}; + +pub fn q(n: &Num) -> Option { + match n.kind() { + NumRef::Exact(d) => Some(BigRational::new( + d.mantissa().into(), + BigInt::from(10).pow(d.scale()), + )), + NumRef::Fraction(r) => Some(BigRational::new((*r.numer()).into(), (*r.denom()).into())), + NumRef::Inexact(x) => BigRational::from_float(x), + } +} + +/// Whether `v` fits a stored exact arm: a decimal at scale <= 28, or an `i128` fraction. +pub fn fits(v: &BigRational) -> bool { + let limit = BigInt::one() << 96; + let decimal = (0..=28u32).any(|s| { + let scaled = v * BigRational::from_integer(BigInt::from(10).pow(s)); + scaled.is_integer() && scaled.to_integer().abs() < limit + }); + let (min, max) = (BigInt::from(i128::MIN), BigInt::from(i128::MAX)); + decimal || (&min <= v.numer() && v.numer() <= &max && v.denom() <= &max) +} + +/// Whether `f` is the round-to-nearest-even f64 of `v`: the review's interval form. +pub fn is_nearest(f: f64, v: &BigRational) -> bool { + let exact = |x: f64| BigRational::from_float(x).expect("finite"); + let two = BigRational::from_integer(BigInt::from(2)); + let half_ulp_past = |x: f64, up: bool| { + let (next, inward) = match up { + true => (x.next_up(), x.next_down()), + false => (x.next_down(), x.next_up()), + }; + match next.is_finite() { + true => (exact(x) + exact(next)) / &two, + false => exact(x) + (exact(x) - exact(inward)) / &two, + } + }; + match f.is_finite() { + false => { + let edge = half_ulp_past(f64::MAX.copysign(f), f > 0.0); + (f > 0.0 && v >= &edge) || (f < 0.0 && v <= &edge) + } + true => { + let lo = half_ulp_past(f, false); + let hi = half_ulp_past(f, true); + let even = f.to_bits() & 1 == 0; + match even { + true => &lo <= v && v <= &hi, + false => &lo < v && v < &hi, + } + } + } +} + +/// A decimal's text built digit by digit, without `rust_decimal`'s `Display`. +pub fn decimal_text(mantissa: i128, scale: u32) -> String { + let digits = format!( + "{:0>width$}", + mantissa.unsigned_abs(), + width = scale as usize + 1 + ); + let (int, frac) = digits.split_at(digits.len() - scale as usize); + let sign = if mantissa < 0 { "-" } else { "" }; + format!("{sign}{int}.{frac}0") +} + +/// The correctly rounded f64 of an exact value, independent of `numtower`'s own. +pub fn independent_f64(n: &Num) -> f64 { + match n.kind() { + NumRef::Exact(d) => decimal_text(d.mantissa(), d.scale()).parse().unwrap(), + NumRef::Fraction(_) => q(n).unwrap().to_f64().unwrap(), + NumRef::Inexact(x) => x, + } +} + +/// Same arm, same value and same scale: what a fingerprint would see. +pub fn same(a: &Num, b: &Num) -> bool { + match (a.kind(), b.kind()) { + (NumRef::Exact(x), NumRef::Exact(y)) => { + x.mantissa() == y.mantissa() && x.scale() == y.scale() + } + (NumRef::Fraction(x), NumRef::Fraction(y)) => x == y, + (NumRef::Inexact(x), NumRef::Inexact(y)) => x.to_bits() == y.to_bits(), + _ => false, + } +} + +pub fn arm(n: &Num) -> &'static str { + match n.kind() { + NumRef::Exact(_) => "exact", + NumRef::Fraction(_) => "fraction", + NumRef::Inexact(_) => "inexact", + } +} + +pub fn mantissa96() -> impl Strategy { + let max = (1i128 << 96) - 1; + prop_oneof![ + -max..=max, + -1_000_000i128..1_000_000, + -(10i128.pow(15))..10i128.pow(15), + ] +} + +pub fn decimal() -> impl Strategy { + (mantissa96(), 0u32..=28).prop_map(|(m, s)| Num::exact(Decimal::from_i128_with_scale(m, s))) +} + +fn of_width(bits: u32) -> impl Strategy { + (any::(), any::()).prop_map(move |(r, neg)| { + let v = ((r >> (128 - bits)) | (1u128 << (bits - 1))) as i128; + if neg { -v } else { v } + }) +} + +pub fn fraction() -> impl Strategy { + let widths = (1u32..=126, 1u32..=126) + .prop_flat_map(|(a, b)| (of_width(a), of_width(b))) + .prop_map(|(n, d)| Num::ratio(n, d).unwrap()); + let edge = (1i128..i128::MAX).prop_map(|d| Num::ratio(i128::MIN, d | 1).unwrap()); + prop_oneof![8 => widths, 1 => edge] +} + +pub fn finite_f64() -> impl Strategy { + use proptest::num::f64::{NEGATIVE, NORMAL, POSITIVE, SUBNORMAL, ZERO}; + prop_oneof![ + (POSITIVE | NEGATIVE | NORMAL | SUBNORMAL | ZERO), + -1e6f64..1e6, + (-1000i64..1000).prop_map(|k| k as f64 / 8.0), + ] +} + +pub fn exact() -> impl Strategy { + prop_oneof![decimal(), fraction()] +} + +pub fn any_num() -> impl Strategy { + prop_oneof![ + 3 => exact(), + 1 => finite_f64().prop_map(Num::inexact), + 1 => prop_oneof![Just(f64::NAN), Just(f64::INFINITY), Just(f64::NEG_INFINITY)].prop_map(Num::inexact), + ] +} + +/// Runs a property; `deterministic` for a rival, so its failure is not luck. +pub fn run( + strategy: S, + cases: u32, + deterministic: bool, + test: impl Fn(S::Value) -> Result<(), TestCaseError>, +) -> Result<(), String> { + let config = Config { + cases, + failure_persistence: None, + ..Config::default() + }; + let mut runner = match deterministic { + true => TestRunner::new_with_rng(config, TestRng::deterministic_rng(RngAlgorithm::ChaCha)), + false => TestRunner::new(config), + }; + runner.run(&strategy, test).map_err(|e| e.to_string()) +} + +/// A rival must fail on the property's own assertion, never on something else. +pub fn assert_caught(result: Result<(), String>, marker: &str) { + let err = result.expect_err("the rival passed: the property cannot see it"); + assert!(err.contains(marker), "failed for another reason: {err}"); +} + +pub fn oracle_cmp(a: &Num, b: &Num) -> Option { + Some(q(a)?.cmp(&q(b)?)) +} + +pub fn integer(v: &BigRational) -> bool { + v.is_integer() && !v.numer().is_zero() +} diff --git a/crates/numtower/tests/properties.rs b/crates/numtower/tests/properties.rs new file mode 100644 index 00000000..5c5a8e5e --- /dev/null +++ b/crates/numtower/tests/properties.rs @@ -0,0 +1,1367 @@ +//! ADR 005 migration step 1's property tests, each beside the rival it must catch. +//! What each rival is and why it matters: `context/test-evidence/numtower.md`. + +mod oracle; + +use std::cmp::Ordering; +use std::str::FromStr; + +use num_bigint::BigInt; +use num_rational::BigRational; +use num_traits::{ToPrimitive, Zero}; +use numtower::{Decimal, Interim, Num, NumRef, NumSum, Typing, is_number, parse_number}; +use oracle::*; +use proptest::prelude::*; + +fn typed(s: &str) -> Num { + parse_number(s, Typing::Inferred).unwrap_or_else(|| panic!("{s} is a number")) +} + +fn dec(s: &str) -> Num { + Num::exact(Decimal::from_str(s).unwrap()) +} + +// --------------------------------------------------------------------------- +// An Exact result is never a rounded one. + +type BinOp<'a> = &'a dyn Fn(&Num, &Num) -> Option; + +fn exact_op(name: &str) -> fn(&BigRational, &BigRational) -> Option { + match name { + "+" => |a, b| Some(a + b), + "-" => |a, b| Some(a - b), + "*" => |a, b| Some(a * b), + _ => |a, b| (b.numer() != &BigInt::from(0)).then(|| a / b), + } +} + +fn check_never_rounded(name: &str, op: BinOp, deterministic: bool) -> Result<(), String> { + let oracle_op = exact_op(name); + run((exact(), exact()), 256, deterministic, |(a, b)| { + if let Some(r) = op(&a, &b).filter(Num::is_exact) { + let want = oracle_op(&q(&a).unwrap(), &q(&b).unwrap()); + prop_assert!( + q(&r) == want, + "an exact result was rounded: {a} {name} {b} = {r}" + ); + } + Ok(()) + }) +} + +fn real(name: &str) -> impl Fn(&Num, &Num) -> Option + '_ { + move |a: &Num, b: &Num| match name { + "+" => Some(a + b), + "-" => Some(a - b), + "*" => Some(a * b), + _ => a.checked_div(b).ok(), + } +} + +#[test] +fn an_exact_result_is_never_a_rounded_one() { + for name in ["+", "-", "*", "/"] { + check_never_rounded(name, &real(name), false).unwrap(); + } + let third = dec("0.3333333333333333333333333333"); + let sum = &dec("1000") + &third; + assert!(sum.is_exact(), "{sum}"); + assert_eq!(q(&sum), Some(q(&dec("1000")).unwrap() + q(&third).unwrap())); +} + +/// `rust_decimal`'s own `+`, trusted as exact: it rounds past 28 digits and says nothing. +fn silently_rounding_add(a: &Num, b: &Num) -> Option { + match (a.kind(), b.kind()) { + (NumRef::Exact(x), NumRef::Exact(y)) => { + Some(x.checked_add(*y).map_or_else(|| a + b, Num::exact)) + } + _ => Some(a + b), + } +} + +#[test] +fn an_exact_result_is_never_a_rounded_one_catches_a_silently_rounding_add() { + assert_caught( + check_never_rounded("+", &silently_rounding_add, true), + "an exact result was rounded", + ); +} + +// --------------------------------------------------------------------------- +// Contagion is total. + +fn check_contagion(add: BinOp, deterministic: bool) -> Result<(), String> { + let seed = Just((dec("2"), Num::inexact(3.0))); + let cases = prop_oneof![1 => seed, 8 => (any_num(), any_num())]; + run(cases, 256, deterministic, |(a, b)| { + let results = [ + ("+", add(&a, &b)), + ("-", Some(&a - &b)), + ("*", Some(&a * &b)), + ("/", a.checked_div(&b).ok()), + ]; + for (name, r) in results { + match (r, a.is_exact() && b.is_exact()) { + (None, _) => prop_assert!( + name == "/" && b.is_zero(), + "contagion: {a} {name} {b} refused" + ), + (Some(r), false) => { + prop_assert!(!r.is_exact(), "contagion: {a} {name} {b} = {r} is exact") + } + (Some(r), true) => { + let v = exact_op(name)(&q(&a).unwrap(), &q(&b).unwrap()).unwrap(); + prop_assert_eq!( + r.is_exact(), + fits(&v), + "contagion: {} {} {} = {}", + a, + name, + b, + r + ); + prop_assert!( + r.is_exact() || is_nearest(r.to_f64(), &v), + "contagion: {a} {name} {b} = {r}" + ); + } + } + } + prop_assert_eq!( + a.checked_div(&b).is_err(), + b.is_zero(), + "contagion: div/0 for {} / {}", + a, + b + ); + Ok(()) + }) +} + +#[test] +fn contagion_is_total() { + check_contagion(&|a, b| Some(a + b), false).unwrap(); +} + +/// Keeps an integral f64 exact, as a "harmless" promotion would. +fn integral_inexact_kept_exact(a: &Num, b: &Num) -> Option { + match (a.is_exact(), b.kind()) { + (true, NumRef::Inexact(x)) if x.is_finite() && x.fract() == 0.0 && x.abs() < 1e15 => { + Some(a + &Num::exact(Decimal::from(x as i64))) + } + _ => Some(a + b), + } +} + +#[test] +fn contagion_is_total_catches_an_integral_inexact_kept_exact() { + assert_caught( + check_contagion(&integral_inexact_kept_exact, true), + "contagion", + ); +} + +/// The scale a decimal power is stored at: minimal, padded toward `min(s·k, 28)` for k > 0 only. +fn power_scale(base: &Num, k: i32, v: &BigRational) -> u32 { + let scaled = |s: u32| v * BigRational::from_integer(BigInt::from(10).pow(s)); + let fits_at = + |s: u32| scaled(s).is_integer() && scaled(s).to_integer().magnitude().bits() <= 96; + let minimal = (0..=28) + .find(|&s| scaled(s).is_integer()) + .expect("a decimal result"); + let target = match (base.kind(), k > 0) { + (NumRef::Exact(d), true) => (d.scale() * k as u32).min(28), + _ => 0, + }; + (minimal..=target.max(minimal)) + .rev() + .find(|&s| fits_at(s)) + .unwrap_or(minimal) +} + +#[test] +fn pow_is_exact_for_an_exact_base_and_integer_exponent_when_it_fits() { + let units = + prop::sample::select(vec!["1.0", "-1.0", "1.000", "-1", "0.00", "2.0"]).prop_map(dec); + let bases = prop_oneof![4 => any_num(), 1 => units]; + run((bases, -6i32..=6), 256, false, |(a, k)| { + let r = a.pow(&Num::exact(Decimal::from(k))); + match (q(&a).filter(|_| a.is_exact()), r) { + (_, Err(_)) => prop_assert!(a.is_zero() && k < 0, "{a}^{k} refused"), + (None, Ok(r)) => prop_assert!(!r.is_exact(), "{a}^{k} = {r} is exact"), + (Some(v), Ok(r)) => { + let want = v.pow(k); + prop_assert_eq!(r.is_exact(), fits(&want), "{}^{} = {}", a, k, r); + prop_assert!( + !r.is_exact() || q(&r) == Some(want.clone()), + "{a}^{k} = {r}" + ); + prop_assert!( + r.is_exact() || is_nearest(r.to_f64(), &want), + "{a}^{k} = {r}" + ); + if let NumRef::Exact(d) = r.kind() { + let scale = power_scale(&a, k, &want); + prop_assert_eq!(d.scale(), scale, "power scale: {}^{} = {}", a, k, r); + } + } + } + Ok(()) + }) + .unwrap(); +} + +#[test] +fn a_negative_power_has_its_minimal_scale() { + let at = |base: &str, k: i64| dec(base).pow(&Num::exact(Decimal::from(k))).unwrap(); + assert!(same(&at("1.0", -2), &dec("1")), "{}", at("1.0", -2)); + assert!(same(&at("-1.0", -3), &dec("-1")), "{}", at("-1.0", -3)); + assert!(same(&at("2.0", -2), &dec("0.25"))); + assert!( + same(&at("-1.0", 3), &dec("-1.000")), + "a positive power keeps s·n" + ); +} + +#[test] +fn pow_decides_its_size_before_computing_and_honours_the_ceiling() { + let third = Num::ratio(1, 3).unwrap(); + let huge = third + .pow(&Num::exact(Decimal::from(1_000_000_000))) + .unwrap(); + assert!(!huge.is_exact() && huge.to_f64() == 0.0, "{huge}"); + let past = third.pow(&Num::exact(Decimal::from(6000))).unwrap(); + assert!( + !past.is_exact(), + "3^6000 has 9,510 bits, past the 8,192-bit ceiling" + ); + let growth = dec("1.05").pow(&Num::exact(Decimal::from(28))).unwrap(); + assert!(growth.is_exact(), "1.05^28 is exact under (d)"); + assert!( + !dec("1.05") + .pow(&Num::exact(Decimal::from(29))) + .unwrap() + .is_exact() + ); + let unit = dec("-1.0") + .pow(&Num::exact(Decimal::from_i128_with_scale( + 10i128.pow(20), + 0, + ))) + .unwrap(); + assert!( + same(&unit, &dec("1.0000000000000000000000000000")), + "{unit}" + ); + assert!(dec("0").pow(&dec("-1")).is_err()); + assert!(!dec("2").pow(&dec("0.5")).unwrap().is_exact()); + assert!(same(&dec("2").pow(&dec("-2")).unwrap(), &dec("0.25"))); + assert!(same(&dec("1.50").pow(&dec("2.0")).unwrap(), &dec("2.2500"))); + assert!(same(&dec("0").pow(&dec("0")).unwrap(), &dec("1"))); +} + +// --------------------------------------------------------------------------- +// The conversion to f64 agrees with str::parse of the decimal text. + +fn check_to_f64( + conv: &dyn Fn(&Decimal) -> f64, + range: i128, + deterministic: bool, +) -> Result<(), String> { + run((-range..=range, 0u32..=28), 256, deterministic, |(m, s)| { + let d = Decimal::from_i128_with_scale(m, s); + let want: f64 = decimal_text(m, s).parse().unwrap(); + prop_assert_eq!( + conv(&d).to_bits(), + want.to_bits(), + "to_f64 disagrees with the text of {}", + d + ); + Ok(()) + }) +} + +const FULL: i128 = (1 << 96) - 1; + +#[test] +fn the_f64_conversion_agrees_with_parsing_the_decimal_text() { + check_to_f64(&|d| Num::exact(*d).to_f64(), FULL, false).unwrap(); + let review = dec("63294.607176721429540940386554"); + assert_eq!( + review.to_f64(), + "63294.607176721429540940386554".parse::().unwrap() + ); +} + +fn crate_to_f64(d: &Decimal) -> f64 { + rust_decimal::prelude::ToPrimitive::to_f64(d).unwrap() +} + +#[test] +fn the_f64_conversion_agrees_with_parsing_the_decimal_text_catches_rust_decimal_to_f64() { + assert_caught(check_to_f64(&crate_to_f64, FULL, true), "to_f64 disagrees"); + let review = Decimal::from_str("63294.607176721429540940386554").unwrap(); + assert_ne!( + crate_to_f64(&review), + "63294.607176721429540940386554".parse::().unwrap() + ); +} + +#[test] +fn rust_decimal_to_f64_agrees_on_fifteen_digit_mantissas_so_a_narrow_generator_is_vacuous() { + let narrow = 10i128.pow(15) - 1; + let result = run((-narrow..=narrow, 0u32..=15), 256, true, |(m, s)| { + let d = Decimal::from_i128_with_scale(m, s); + prop_assert_eq!( + crate_to_f64(&d).to_bits(), + decimal_text(m, s).parse::().unwrap().to_bits() + ); + Ok(()) + }); + result.unwrap(); +} + +// --------------------------------------------------------------------------- +// A fraction's conversion lies inside the round-to-nearest-even interval. + +fn wide_fraction() -> impl Strategy { + ( + 20u32..=126, + 20u32..=126, + any::(), + any::(), + any::(), + ) + .prop_map(|(a, b, x, y, neg)| { + let n = ((x >> (128 - a)) | (1 << (a - 1))) as i128; + let d = ((y >> (128 - b)) | (1 << (b - 1))) as i128; + Num::ratio(if neg { -n } else { n }, d).unwrap() + }) +} + +fn check_fraction_conversion( + conv: &dyn Fn(&Num) -> f64, + deterministic: bool, +) -> Result<(), String> { + run( + (wide_fraction(), wide_fraction()), + 256, + deterministic, + |(x, y)| { + prop_assert!( + is_nearest(conv(&x), &q(&x).unwrap()), + "outside the interval: {x}" + ); + let wide = (&Interim::from(x.clone()) * &Interim::from(y.clone())).store(); + let exact = q(&x).unwrap() * q(&y).unwrap(); + prop_assert!( + wide.is_exact() || is_nearest(wide.to_f64(), &exact), + "outside the interval: {x} * {y}" + ); + Ok(()) + }, + ) +} + +#[test] +fn a_fractions_conversion_lies_inside_the_round_to_nearest_even_interval() { + check_fraction_conversion(&Num::to_f64, false).unwrap(); +} + +fn divide_in_f64(x: &Num) -> f64 { + match x.kind() { + NumRef::Fraction(r) => *r.numer() as f64 / *r.denom() as f64, + _ => x.to_f64(), + } +} + +#[test] +fn a_fractions_conversion_lies_inside_the_round_to_nearest_even_interval_catches_an_f64_division() { + assert_caught( + check_fraction_conversion(÷_in_f64, true), + "outside the interval", + ); +} + +// --------------------------------------------------------------------------- +// Exact mixed comparison is transitive. + +/// Values crowded around one f64: equal-looking decimals, the f64 and its neighbours. +fn crowd() -> impl Strategy> { + (-1e6f64..1e6, 1i64..9).prop_map(|(x, k)| { + let shortest = typed(&format!("{x:?}")); + let nudge = Num::exact(Decimal::new(k, 25)); + vec![ + &shortest + &nudge, + Num::inexact(x), + &shortest - &nudge, + shortest.clone(), + Num::inexact(x.next_up()), + Num::inexact(x.next_down()), + shortest.checked_div(&Num::exact(Decimal::from(3))).unwrap(), + ] + }) +} + +type Cmp<'a> = &'a dyn Fn(&Num, &Num) -> Option; + +fn check_transitive(cmp: Cmp, deterministic: bool) -> Result<(), String> { + run(crowd(), 256, deterministic, |xs| { + for a in &xs { + for b in &xs { + let reversed = cmp(b, a).map(Ordering::reverse); + prop_assert_eq!( + cmp(a, b), + reversed, + "comparison: not antisymmetric: {} vs {}", + a, + b + ); + for c in &xs { + let (ab, bc, ac) = (cmp(a, b), cmp(b, c), cmp(a, c)); + let chained = ab.is_some_and(|o| o != Ordering::Greater) + && bc.is_some_and(|o| o != Ordering::Greater); + let strict = ab == Some(Ordering::Less) || bc == Some(Ordering::Less); + let want = if strict { + Ordering::Less + } else { + Ordering::Equal + }; + prop_assert!( + !chained || ac == Some(want), + "comparison: not transitive: {a}, {b}, {c}" + ); + } + } + } + Ok(()) + }) +} + +#[test] +fn exact_mixed_comparison_is_transitive() { + check_transitive(&Num::cmp_exact, false).unwrap(); + assert_ne!(dec("0.1"), Num::inexact(0.1)); +} + +#[test] +fn comparison_agrees_with_the_oracle_at_every_pairing_of_arms() { + let finite = prop_oneof![exact(), finite_f64().prop_map(Num::inexact)]; + let pools = prop_oneof![crowd(), prop::collection::vec(finite, 2..6)]; + run(pools, 256, false, |xs| { + for a in &xs { + for b in &xs { + prop_assert_eq!(a.cmp_exact(b), oracle_cmp(a, b), "{} vs {}", a, b); + } + } + Ok(()) + }) + .unwrap(); + assert_eq!(Num::inexact(f64::NAN).cmp_exact(&dec("1")), None); + assert_eq!( + dec("1").cmp_exact(&Num::inexact(f64::INFINITY)), + Some(Ordering::Less) + ); +} + +/// Revision 1's rule: the exact operand converted to f64, then compared. +fn compare_through_f64(a: &Num, b: &Num) -> Option { + match a.is_exact() == b.is_exact() { + true => a.cmp_exact(b), + false => a.to_f64().partial_cmp(&b.to_f64()), + } +} + +#[test] +fn exact_mixed_comparison_is_transitive_catches_comparing_through_f64() { + assert_caught( + check_transitive(&compare_through_f64, true), + "comparison: not transitive", + ); +} + +#[test] +fn min_and_max_prefer_the_exact_value_and_do_not_depend_on_order() { + let half = dec("0.50"); + let pairs = [ + (half.clone(), Num::inexact(0.5)), + (dec("1.5"), dec("1.50")), + (Num::inexact(0.0), Num::inexact(-0.0)), + (Num::inexact(f64::NAN), dec("2")), + ]; + for (a, b) in &pairs { + assert!(same(a.min_of(b), b.min_of(a)), "min of {a}, {b}"); + assert!(same(a.max_of(b), b.max_of(a)), "max of {a}, {b}"); + } + assert!(same(half.min_of(&Num::inexact(0.5)), &half)); + assert!(same(dec("1.5").max_of(&dec("1.50")), &dec("1.50"))); + assert!( + Num::inexact(0.0) + .min_of(&Num::inexact(-0.0)) + .to_f64() + .is_sign_negative() + ); + assert!( + !Num::inexact(-0.0) + .max_of(&Num::inexact(0.0)) + .to_f64() + .is_sign_negative() + ); + assert!(same(Num::inexact(f64::NAN).min_of(&dec("2")), &dec("2"))); +} + +#[test] +fn the_accumulator_keeps_inexact_terms_whose_weights_cancel() { + let mut s = NumSum::new(); + s.add(&Num::inexact(1.0), 1); + s.add(&Num::inexact(2.0), -1); + s.insert(&dec("0.5")); + let read = s.read(); + assert!(!read.is_exact(), "an inexact term is still present: {read}"); + assert_eq!(read.to_f64(), -0.5); + s.add(&Num::inexact(1.0), -1); + s.add(&Num::inexact(2.0), 1); + assert!(same(&s.read(), &dec("0.5")), "every inexact term retracted"); +} + +#[test] +fn a_huge_odd_exponent_keeps_the_sign_of_a_negative_base() { + let odd = Num::exact(Decimal::from(9_007_199_254_740_993i64)); + let up = dec("-2").pow(&odd).unwrap(); + assert!(!up.is_exact() && up.to_f64() == f64::NEG_INFINITY, "{up}"); + let down = dec("-2").pow(&-&odd).unwrap(); + assert!( + !down.is_exact() && down.to_f64() == 0.0 && down.to_f64().is_sign_negative(), + "{down}" + ); + let third = Num::ratio(-1, 3).unwrap().pow(&odd).unwrap(); + assert!(third.to_f64().is_sign_negative(), "{third}"); + let even = &odd + &dec("1"); + assert_eq!(dec("-2").pow(&even).unwrap().to_f64(), f64::INFINITY); +} + +#[test] +fn i128_min_is_an_exact_fraction_numerator() { + let min = Num::ratio(i128::MIN, 1).unwrap(); + assert!( + matches!(min.kind(), NumRef::Fraction(r) if *r.numer() == i128::MIN), + "{min}" + ); + let third = Num::ratio(i128::MIN, 3).unwrap(); + assert!(matches!(third.kind(), NumRef::Fraction(_)), "{third}"); + let negated = -&third; + let want = -q(&third).unwrap(); + assert!( + !negated.is_exact() && is_nearest(negated.to_f64(), &want), + "{negated}" + ); + let back = + &(&Interim::from(third.clone()) * &Interim::from(dec("-1"))) * &Interim::from(dec("-1")); + assert!( + same(&back.store(), &third), + "through the wide path and back" + ); + let halved = third.checked_div(&dec("2")).unwrap(); + assert_eq!( + q(&halved), + Some(q(&third).unwrap() / BigRational::from_integer(BigInt::from(2))) + ); + assert_eq!( + min.cmp_exact(&Num::ratio(i128::MIN + 1, 1).unwrap()), + Some(Ordering::Less) + ); + assert_eq!( + Num::from_canonical(&min.to_string()).map(|n| same(&n, &min)), + Some(true) + ); +} + +// --------------------------------------------------------------------------- +// The accumulator is independent of order. + +fn term() -> impl Strategy { + let pooled = ( + -1000i128..1000, + prop::sample::select(vec![3i128, 7, 12, 49, 1 << 40, 3 * (1 << 33)]), + ) + .prop_map(|(n, d)| Num::ratio(n, d).unwrap()); + prop_oneof![ + 4 => decimal(), + 2 => pooled, + 1 => fraction(), + 1 => finite_f64().prop_map(Num::inexact), + ] +} + +/// Z-set weights: zero, unit, small, wide and both extremes (never `i64::MIN`, so it negates). +fn weight() -> impl Strategy { + prop_oneof![ + 3 => Just(1i64), + 1 => Just(-1i64), + 1 => Just(0i64), + 2 => -3i64..=3, + 2 => (-(1i64 << 40))..(1i64 << 40), + 1 => prop_oneof![Just(i64::MAX), Just(-i64::MAX), Just((1i64 << 30) - 1)], + ] +} + +/// Many non-decimal denominators, so a list often holds more than K groups, several terms each. +const POOL: [i128; 14] = [ + 3, + 7, + 11, + 12, + 13, + 17, + 19, + 23, + 29, + 31, + 37, + 49, + 1 << 40, + 3 * (1 << 33), +]; + +fn pooled_term() -> impl Strategy { + let pooled = (-1000i128..1000, prop::sample::select(POOL.to_vec())) + .prop_map(|(n, d)| Num::ratio(n, d).unwrap()); + prop_oneof![ + 3 => decimal(), + 4 => pooled, + 1 => fraction(), + 2 => finite_f64().prop_map(Num::inexact), + ] +} + +fn weighted(terms: impl Strategy) -> impl Strategy { + (terms, weight()) +} + +fn sum_all(xs: &[Num]) -> Num { + let mut s = NumSum::new(); + xs.iter().for_each(|x| s.insert(x)); + s.read() +} + +fn sum_weighted(xs: &[(Num, i64)]) -> Num { + let mut s = NumSum::new(); + xs.iter().for_each(|(x, w)| s.add(x, *w)); + s.read() +} + +type WeightedSum<'a> = &'a dyn Fn(&[(Num, i64)]) -> Num; + +fn check_order_free(sum: WeightedSum, deterministic: bool) -> Result<(), String> { + let seed = vec![ + (dec("0.30000000000000004"), 1), + (dec("1000000000000"), 1), + (dec("-1000000000000"), 1), + ]; + let cases = prop_oneof![ + 1 => Just(seed), + 4 => prop::collection::vec(weighted(term()), 0..24), + 2 => prop::collection::vec(weighted(pooled_term()), 0..40), + ]; + let shuffled = cases.prop_flat_map(|xs| (Just(xs.clone()), Just(xs).prop_shuffle())); + run(shuffled, 256, deterministic, |(xs, ys)| { + let (a, b) = (sum(&xs), sum(&ys)); + prop_assert!(same(&a, &b), "depends on order: {a} against {b}"); + Ok(()) + }) +} + +#[test] +fn the_accumulator_is_independent_of_order() { + check_order_free(&sum_weighted, false).unwrap(); + let seed = sum_all(&[ + dec("0.30000000000000004"), + dec("1000000000000"), + dec("-1000000000000"), + ]); + assert!(same(&seed, &dec("0.30000000000000004")), "{seed}"); +} + +/// Clause 8's rival: `Decimal`'s own `+`, demoting to f64 where it overflows. +fn decimal_fold(xs: &[(Num, i64)]) -> Num { + xs.iter().fold(Num::exact(Decimal::ZERO), |acc, (x, w)| { + let x = x * &Num::exact(Decimal::from(*w)); + match (acc.kind(), x.kind()) { + (NumRef::Exact(a), NumRef::Exact(b)) => { + a.checked_add(*b).map_or_else(|| &acc + &x, Num::exact) + } + _ => &acc + &x, + } + }) +} + +#[test] +fn the_accumulator_is_independent_of_order_catches_decimal_addition() { + assert_caught(check_order_free(&decimal_fold, true), "depends on order"); +} + +#[test] +fn the_accumulator_retracts_exactly() { + let terms = || prop::collection::vec(weighted(pooled_term()), 0..20); + run((terms(), terms()), 256, false, |(xs, ys)| { + let mut s = NumSum::new(); + xs.iter().chain(&ys).for_each(|(x, w)| s.add(x, *w)); + ys.iter().rev().for_each(|(y, w)| s.add(y, -w)); + let fresh = sum_weighted(&xs); + prop_assert!( + same(&s.read(), &fresh), + "{} after retraction, {} fresh", + s.read(), + fresh + ); + Ok(()) + }) + .unwrap(); +} + +/// The oracle's partition of a weighted list, exactly as clause 8 and option (d) state it. +struct Expected { + decimal_part: BigRational, + groups: std::collections::BTreeMap, + inexact: std::collections::BTreeMap, + inexact_total: BigRational, +} + +fn ratio_of(i: i128) -> BigRational { + BigRational::from_integer(BigInt::from(i)) +} + +fn rounded(v: &BigRational) -> BigRational { + BigRational::from_float(v.to_f64().unwrap()).unwrap() +} + +impl Expected { + fn of(xs: &[(Num, i64)]) -> Expected { + let zero = || ratio_of(0); + let mut e = Expected { + decimal_part: zero(), + groups: Default::default(), + inexact: Default::default(), + inexact_total: zero(), + }; + for (x, w) in xs { + let weighted = q(x).unwrap() * ratio_of(*w as i128); + match x.kind() { + NumRef::Inexact(f) => { + *e.inexact.entry(f.to_bits()).or_insert(0) += *w as i128; + e.inexact_total += weighted; + } + NumRef::Fraction(r) if !decimal_denominator(*r.denom()) => { + *e.groups.entry(*r.denom()).or_insert_with(zero) += weighted; + } + _ => e.decimal_part += weighted, + } + } + e + } + + fn live(&self) -> impl Iterator { + self.groups.values().filter(|v| !v.is_zero()) + } +} + +#[test] +fn the_accumulator_reads_the_exact_sum_or_its_stated_roundings() { + let arms = std::cell::Cell::new([0u32; 3]); + let lists = prop_oneof![ + prop::collection::vec(weighted(term()), 0..24), + prop::collection::vec(weighted(pooled_term()), 0..40), + ]; + run(lists, 256, false, |xs| { + let got = sum_weighted(&xs); + let e = Expected::of(&xs); + let present = e.inexact.values().any(|w| *w != 0); + let exact_total = e.live().fold(e.decimal_part.clone(), |s, g| s + g); + let mut seen = arms.get(); + match (present, e.live().count() <= numtower::K) { + (false, true) => { + seen[0] += 1; + prop_assert_eq!(got.is_exact(), fits(&exact_total), "{}", got); + prop_assert!( + !got.is_exact() || q(&got) == Some(exact_total.clone()), + "{got}" + ); + prop_assert!( + got.is_exact() || is_nearest(got.to_f64(), &exact_total), + "{got}" + ); + } + (true, true) => { + seen[1] += 1; + let total = rounded(&exact_total) + &e.inexact_total; + prop_assert!(!got.is_exact() && is_nearest(got.to_f64(), &total), "{got}"); + } + (_, false) => { + seen[2] += 1; + let groups = e.live().map(rounded).fold(ratio_of(0), |s, g| s + g); + let total = groups + rounded(&e.decimal_part) + &e.inexact_total; + let near = is_nearest(got.to_f64(), &total); + prop_assert!(!got.is_exact() && near, "past K groups: {got} for {total}"); + } + } + arms.set(seen); + Ok(()) + }) + .unwrap(); + let seen = arms.get(); + assert!( + seen.iter().all(|n| *n >= 20), + "every branch must be reached often: {seen:?}" + ); +} + +fn decimal_denominator(d: i128) -> bool { + let p = d.trailing_zeros(); + let mut rest = d >> p; + let mut q5 = 0; + while rest % 5 == 0 { + rest /= 5; + q5 += 1; + } + rest == 1 && p.max(q5) <= 28 +} + +#[test] +fn more_than_k_denominators_read_inexact_even_when_the_total_fits() { + let primes = [3i128, 7, 11, 13, 17, 19, 23, 29, 31]; + let terms: Vec = primes.iter().map(|&p| Num::ratio(1, p).unwrap()).collect(); + let eight = sum_all(&terms[..8]); + assert!(matches!(eight.kind(), NumRef::Fraction(_)), "{eight}"); + let nine = sum_all(&terms); + assert!(!nine.is_exact(), "{nine}"); + let exact: BigRational = terms.iter().map(|t| q(t).unwrap()).sum(); + assert!( + fits(&exact), + "the exact total fits, so the cutoff is what demotes" + ); + let each_rounded: BigRational = terms + .iter() + .map(|t| BigRational::from_float(q(t).unwrap().to_f64().unwrap()).unwrap()) + .sum(); + assert!( + is_nearest(nine.to_f64(), &each_rounded), + "each group rounded, then the sum once" + ); +} + +#[test] +fn exact_sum_is_the_exact_weighted_total_rounded_once() { + let weights = prop_oneof![-3i64..=3, any::(), (-(1i64 << 40))..(1i64 << 40)]; + run( + prop::collection::vec((finite_f64(), weights), 0..20), + 256, + false, + |terms| { + let mut acc = numtower::ExactSum::new(); + let mut exact = BigRational::from_integer(BigInt::from(0)); + for (x, w) in &terms { + acc.add(*x, *w); + exact += BigRational::from_float(*x).unwrap() + * BigRational::from_integer(BigInt::from(*w)); + } + prop_assert!( + is_nearest(acc.value(), &exact), + "{} for {:?}", + acc.value(), + terms + ); + Ok(()) + }, + ) + .unwrap(); +} + +// --------------------------------------------------------------------------- +// Division's direct fit test. + +fn check_division(div: BinOp, deterministic: bool) -> Result<(), String> { + let seeds = prop_oneof![ + Just((dec("1"), dec("1099511627776"))), + Just((dec("1"), dec("931322574615478515625"))), + Just((dec("0.0000000000000000000000000001"), dec("0.5"))), + ]; + let cases = prop_oneof![1 => seeds, 3 => (decimal(), decimal())]; + run(cases, 256, deterministic, |(a, b)| { + if let (Some(r), Some(v)) = (div(&a, &b), exact_op("/")(&q(&a).unwrap(), &q(&b).unwrap())) { + prop_assert_eq!( + r.is_exact(), + fits(&v), + "division fit: {} / {} = {}", + a, + b, + r + ); + prop_assert!( + !r.is_exact() || q(&r) == Some(v.clone()), + "division fit: {a} / {b} = {r}" + ); + let decimal = (0..=28u32).any(|s| { + let scaled = &v * BigRational::from_integer(BigInt::from(10).pow(s)); + scaled.is_integer() && scaled.to_integer().magnitude().bits() <= 96 + }); + prop_assert_eq!( + matches!(r.kind(), NumRef::Exact(_)), + decimal, + "division fit: {} / {} = {}", + a, + b, + r + ); + } + Ok(()) + }) +} + +#[test] +fn division_decides_fit_directly() { + check_division(&|a, b| a.checked_div(b).ok(), false).unwrap(); + let frac = |n: &Num| matches!(n.kind(), NumRef::Fraction(_)); + assert!( + frac(&dec("1").checked_div(&dec("1099511627776")).unwrap()), + "1/2^40" + ); + assert!( + frac(&dec("1").checked_div(&dec("931322574615478515625")).unwrap()), + "1/5^30" + ); + let tiny = dec("0.0000000000000000000000000001") + .checked_div(&dec("0.5")) + .unwrap(); + assert!( + same(&tiny, &dec("0.0000000000000000000000000002")), + "{tiny}" + ); + let avg = &(&(&dec("10.00") + &dec("20.00")) + &dec("40.00")); + assert!(frac(&avg.checked_div(&dec("3")).unwrap())); + assert!(same( + &dec("1").checked_div(&dec("4")).unwrap(), + &dec("0.25") + )); + let third = Num::ratio(1, 3).unwrap(); + assert!( + same(&(&third * &dec("3")), &dec("1")), + "1/3 * 3 is exactly 1" + ); + assert_eq!( + Num::inexact(1.0).checked_div(&Num::inexact(-0.0)), + Err(numtower::DivisionByZero) + ); +} + +/// Multiplying back, with revision 6's scale test: sound, and incomplete. +fn multiply_back(a: &Num, b: &Num) -> Option { + match (a.kind(), b.kind()) { + (NumRef::Exact(x), NumRef::Exact(y)) if !y.is_zero() => { + let candidate = x.checked_div(*y)?; + let back = candidate.checked_mul(*y)?; + let kept = back.scale() == candidate.scale() + y.scale() && back == *x; + Some(match kept { + true => Num::exact(candidate), + false => Num::inexact(a.to_f64() / b.to_f64()), + }) + } + _ => a.checked_div(b).ok(), + } +} + +#[test] +fn division_decides_fit_directly_catches_multiplying_back() { + assert_caught(check_division(&multiply_back, true), "division fit"); +} + +// --------------------------------------------------------------------------- +// The number language: the gate and the parse agree. + +fn lexeme() -> impl Strategy { + let digits = "[0-9]{0,40}"; + let structured = ( + prop::option::of("[+-]"), + digits, + prop::option::of(("\\.", digits)), + prop::option::of(("[eE]", prop::option::of("[+-]"), "[0-9]{0,4}")), + ) + .prop_map(|(sign, int, frac, exp)| { + let frac = frac.map(|(p, f)| format!("{p}{f}")).unwrap_or_default(); + let exp = exp + .map(|(e, s, d)| format!("{e}{}{d}", s.unwrap_or_default())) + .unwrap_or_default(); + format!("{}{int}{frac}{exp}", sign.unwrap_or_default()) + }); + prop_oneof![ + 4 => structured, + 1 => "[0-9eE+\\-._/aifnIN ]{0,10}", + 1 => "[+-]?[0-9]{1,30}/[0-9]{1,30}", + 1 => prop::sample::select(vec!["inf", "-inf", "nan", "NaN", "Infinity", "-infinity", "1_000", "+.5", "1.", ".", "e5"]).prop_map(String::from), + ] +} + +fn check_language( + parse: &dyn Fn(&str, Typing) -> Option, + deterministic: bool, +) -> Result<(), String> { + run(lexeme(), 256, deterministic, |s| { + let gate = s.parse::().is_ok_and(f64::is_finite); + let parsed = parse(&s, Typing::Inferred); + prop_assert_eq!( + parsed.is_some(), + gate, + "the gate and the parse disagree on {:?}", + s + ); + prop_assert_eq!( + is_number(&s, Typing::Inferred), + gate, + "the gate and the parse disagree on {:?}", + s + ); + if let Some(n) = &parsed { + let nearest: f64 = s.parse().unwrap(); + prop_assert!( + n.to_f64() == nearest, + "the gate and the parse disagree on the value of {s:?}: {n}" + ); + let v = q(n).unwrap(); + prop_assert!( + !n.is_exact() || is_nearest(nearest, &v), + "the gate and the parse disagree on {s:?}" + ); + } + let forced = parse(&s, Typing::Forced); + prop_assert_eq!( + forced.is_some(), + is_number(&s, Typing::Forced), + "the gate and the parse disagree on forced {:?}", + s + ); + Ok(()) + }) +} + +#[test] +fn the_number_languages_gate_and_parse_agree() { + check_language(&parse_number, false).unwrap(); +} + +/// `rust_decimal`'s exact parse: it refuses `1.5e3` and accepts `1_000`. +fn from_str_exact(s: &str, typing: Typing) -> Option { + match Decimal::from_str_exact(s) { + Ok(d) => Some(Num::exact(d)), + Err(_) if typing == Typing::Forced => parse_number(s, typing), + Err(_) => None, + } +} + +#[test] +fn the_number_languages_gate_and_parse_agree_catches_from_str_exact() { + assert_caught( + check_language(&from_str_exact, true), + "the gate and the parse disagree", + ); +} + +#[test] +fn the_number_language_seed_cases() { + let inferred = |s: &str| parse_number(s, Typing::Inferred); + let forced = |s: &str| parse_number(s, Typing::Forced); + assert!(same(&typed("1.5e3"), &dec("1500"))); + assert!(same(&typed("1.50"), &dec("1.50"))); + assert_eq!(inferred("1_000").map(|n| n.to_string()), None); + let zero = typed("-0.0"); + assert!(same(&zero, &dec("0.0")) && !zero.to_f64().is_sign_negative()); + assert!(inferred("3/4").is_none()); + assert!(same(&forced("3/4").unwrap(), &dec("0.75"))); + assert!(matches!( + forced("-10/3").unwrap().kind(), + NumRef::Fraction(_) + )); + assert!(forced("3/0").is_none()); + for text in ["1e999", "-1e999", "nan", "inf", "Infinity", "NaN", "-inf"] { + assert!( + inferred(text).is_none() && forced(text).is_none(), + "{text} is text" + ); + } + let underflow = typed("1e-999"); + assert!(!underflow.is_exact() && underflow.to_f64() == 0.0); + assert_eq!(typed("1e-30"), Num::ratio(1, 10i128.pow(30)).unwrap()); + assert!(matches!(typed("1e-30").kind(), NumRef::Fraction(_))); + let long = typed("0.12345678901234567890123456789"); + let want = BigRational::new( + BigInt::from(12345678901234567890123456789i128), + BigInt::from(10).pow(29), + ); + assert!(matches!(long.kind(), NumRef::Fraction(_)) && q(&long) == Some(want)); + let big = (1i128 << 120).to_string(); + let n = typed(&big); + assert!(matches!(n.kind(), NumRef::Fraction(r) if *r.denom() == 1 && *r.numer() == 1 << 120)); + assert!(!typed("1e-40").is_exact()); + assert!(!typed("1234567890123456789012345678901234567890").is_exact()); + assert!(same(&typed("0e99999999999999999999"), &dec("0"))); +} + +// --------------------------------------------------------------------------- +// A typed value and the same value computed are equally exact. + +fn small_decimal() -> impl Strategy { + let mantissa = prop_oneof![-1000i128..1000, mantissa96()]; + (mantissa, 0u32..=28).prop_map(|(m, s)| Num::exact(Decimal::from_i128_with_scale(m, s))) +} + +type Mul<'a> = &'a dyn Fn(&Num, &Num) -> Num; + +fn check_typed_equals_computed(mul: Mul, deterministic: bool) -> Result<(), String> { + let seed = Just((dec("0.000000000000001"), dec("0.000000000000001"))); + run( + prop_oneof![1 => seed, 3 => (small_decimal(), small_decimal())], + 256, + deterministic, + |(a, b)| { + let v = q(&a).unwrap() * q(&b).unwrap(); + let ten = |k: u32| BigRational::from_integer(BigInt::from(10).pow(k)); + let k = (0..=60u32) + .find(|&k| (&v * ten(k)).is_integer()) + .expect("a decimal product"); + let text = format!("{}e-{k}", (&v * ten(k)).to_integer()); + let parsed = typed(&text); + let computed = mul(&a, &b); + prop_assert_eq!( + arm(&parsed), + arm(&computed), + "typed and computed differ: {} against {} * {}", + text, + a, + b + ); + prop_assert!( + parsed.cmp_exact(&computed) == Some(Ordering::Equal), + "typed and computed differ: {text} against {computed}" + ); + Ok(()) + }, + ) +} + +#[test] +fn a_typed_value_and_the_same_value_computed_are_equally_exact() { + check_typed_equals_computed(&|a, b| a * b, false).unwrap(); + let computed = &typed("1e-15") * &typed("1e-15"); + assert!(same(&typed("1e-30"), &computed), "{computed}"); +} + +/// Revision 6's rule: demote whenever the returned scale drops. +fn demote_on_the_hint(a: &Num, b: &Num) -> Num { + match (a.kind(), b.kind()) { + (NumRef::Exact(x), NumRef::Exact(y)) => match x.checked_mul(*y) { + Some(r) if r.scale() == x.scale() + y.scale() => Num::exact(r), + _ => Num::inexact(a.to_f64() * b.to_f64()), + }, + _ => a * b, + } +} + +#[test] +fn a_typed_value_and_the_same_value_computed_are_equally_exact_catches_demoting_on_the_hint() { + assert_caught( + check_typed_equals_computed(&demote_on_the_hint, true), + "typed and computed differ", + ); +} + +// --------------------------------------------------------------------------- +// Clause 4's stored scale, question 13 (b)'s widening, and the canonical text. + +#[test] +fn a_recomputed_result_keeps_the_stored_scale_clause_4_names() { + let product = &dec("1.2500000000000000") * &dec("4.0000000000000000"); + assert!( + same(&product, &dec("5.0000000000000000000000000000")), + "{product}" + ); + let common = &dec("12.50") + &dec("1.5"); + assert!(same(&common, &dec("14.00")), "{common}"); + let rounded = &dec("1000") + &dec("0.3333333333333333333333333333"); + assert!(matches!(rounded.kind(), NumRef::Fraction(_)), "{rounded}"); + assert!(same(&(&dec("-1.5") + &dec("1.5")), &dec("0.0"))); +} + +#[test] +fn an_intermediate_beyond_i128_comes_back_exact() { + let a = Num::ratio(10i128.pow(30), 7).unwrap(); + let wide = &Interim::from(a.clone()) * &Interim::from(a.clone()); + assert!( + wide.is_exact(), + "an intermediate is unbounded up to the ceiling" + ); + let back = wide.checked_div(&Interim::from(a.clone())).unwrap().store(); + assert!(same(&back, &a), "{back}"); + let stored = (&a * &a).to_f64(); + assert!(!(&a * &a).is_exact() && is_nearest(stored, &(q(&a).unwrap() * q(&a).unwrap()))); +} + +#[test] +fn the_canonical_text_reads_back_as_the_same_number() { + run(any_num(), 256, false, |n| { + let text = n.to_string(); + let back = Num::from_canonical(&text); + prop_assert!( + back.as_ref() + .is_some_and(|b| same(b, &n) || (n.to_f64().is_nan() && b.to_f64().is_nan())), + "{text}" + ); + Ok(()) + }) + .unwrap(); + assert_eq!(dec("12.50").to_string(), "12.50"); + assert_eq!(Num::ratio(-10, 3).unwrap().to_string(), "-10/3"); + assert_eq!(Num::inexact(1.0 / 3.0).to_string(), "#i0.3333333333333333"); + assert_eq!(Num::inexact(f64::NEG_INFINITY).to_string(), "-inf"); + assert!(Num::from_canonical("#inan").is_none()); + assert!( + Num::from_canonical("1e-999").is_none(), + "canonical exact text never reads inexact" + ); + let ten35 = Num::ratio(10i128.pow(35), 1).unwrap(); + assert_eq!(ten35.to_string(), format!("{}/1", 10i128.pow(35))); +} + +#[test] +fn exact_sums_are_equal_when_their_values_are() { + let of = |terms: &[(f64, i64)]| { + let mut s = numtower::ExactSum::new(); + terms.iter().for_each(|(x, w)| s.add(*x, *w)); + s + }; + assert_eq!( + of(&[(0.5, 1), (0.5, 1)]), + of(&[(1.0, 1)]), + "0.5 + 0.5 = 1.0" + ); + assert_eq!( + of(&[(3.25, 1), (3.25, -1)]), + numtower::ExactSum::new(), + "insert then retract" + ); + let big = (1i64 << 30) - 1; + assert_eq!( + of(&[(1.5, big), (0.25, 2)]), + of(&[(0.25, 2), (1.5, big)]), + "order" + ); + assert_ne!(of(&[(1.0, 1)]), of(&[(1.0, 2)])); + assert_ne!( + of(&[(f64::NAN, 1)]), + of(&[]), + "a counted NaN is part of the value" + ); +} + +#[test] +fn a_forced_fraction_never_reads_as_infinity() { + let huge = format!("{}/1", "9".repeat(400)); + assert_eq!( + parse_number(&huge, Typing::Forced).map(|n| n.to_string()), + None + ); + assert!( + !is_number(&huge, Typing::Forced), + "the gate and the parse agree" + ); + let negative = format!("-{huge}"); + assert!(parse_number(&negative, Typing::Forced).is_none()); + let tiny = format!("1/{}", "9".repeat(400)); + let zero = parse_number(&tiny, Typing::Forced).unwrap(); + assert!( + !zero.is_exact() && zero.to_f64() == 0.0, + "an underflow is a finite zero, as 1e-999" + ); + assert!(is_number(&tiny, Typing::Forced)); +} + +type InterimCmp<'a> = &'a dyn Fn(&Interim, &Interim) -> Option; + +fn check_interim_order(cmp: InterimCmp, deterministic: bool) -> Result<(), String> { + let eps = Interim::from(Num::ratio(1, 3 * 10i128.pow(30)).unwrap()); + run( + ( + wide_fraction(), + wide_fraction(), + wide_fraction(), + wide_fraction(), + ), + 256, + deterministic, + |(x1, x2, y1, y2)| { + let a = &Interim::from(x1.clone()) * &Interim::from(x2.clone()); + let av = q(&x1).unwrap() * q(&x2).unwrap(); + let near = av.to_f64().unwrap(); + let eps_v = + BigRational::new(BigInt::from(1), BigInt::from(3) * BigInt::from(10).pow(30)); + let others = [ + ( + &Interim::from(y1.clone()) * &Interim::from(y2.clone()), + q(&y1).unwrap() * q(&y2).unwrap(), + ), + (&a + &eps, &av + &eps_v), + (&a - &eps, &av - &eps_v), + (a.clone(), av.clone()), + ( + Interim::from(Num::inexact(near)), + BigRational::from_float(near).unwrap(), + ), + ( + Interim::from(Num::inexact(near.next_up())), + BigRational::from_float(near.next_up()).unwrap(), + ), + ]; + for (b, bv) in &others { + prop_assert_eq!( + cmp(&a, b), + Some(av.cmp(bv)), + "interim order: {:?} vs {:?}", + a, + b + ); + prop_assert_eq!( + cmp(b, &a), + Some(bv.cmp(&av)), + "interim order: {:?} vs {:?}", + b, + a + ); + } + let twin = a.clone(); + prop_assert!( + !a.is_zero() && (&a - &twin).is_zero(), + "interim zero: {a:?}" + ); + Ok(()) + }, + ) +} + +#[test] +fn an_interim_compares_exactly_without_storing() { + check_interim_order(&Interim::cmp_exact, false).unwrap(); + assert_eq!( + Interim::from(Num::inexact(f64::NAN)).cmp_exact(&Interim::from(dec("1"))), + None + ); +} + +/// Storing first: what `=` would do if `Interim` could not compare. +fn store_then_compare(a: &Interim, b: &Interim) -> Option { + a.clone().store().cmp_exact(&b.clone().store()) +} + +#[test] +fn an_interim_compares_exactly_without_storing_catches_storing_first() { + assert_caught( + check_interim_order(&store_then_compare, true), + "interim order", + ); +} diff --git a/nix/ci.nix b/nix/ci.nix index ceb4f384..f8a56de9 100644 --- a/nix/ci.nix +++ b/nix/ci.nix @@ -46,6 +46,8 @@ "improvise-wasm-client" "improvise-worker-server" "improvise-gp" + # ADR 005 clause 14: nothing depends on it yet, and improvise-core will. + "numtower" ]; pFlags = crates: lib.concatMapStringsSep " " (c: "-p ${c}") crates; diff --git a/scripts/build-doc-site.mjs b/scripts/build-doc-site.mjs index 1cf88745..cede9610 100644 --- a/scripts/build-doc-site.mjs +++ b/scripts/build-doc-site.mjs @@ -64,6 +64,7 @@ const SECTIONS = [ ["context/code/tui.md", "code/tui", "The TUI"], ["context/code/cli.md", "code/cli", "The headless CLI"], ["context/code/gp.md", "code/gp", "gp-tree and improvise-gp"], + ["context/code/numtower.md", "code/numtower", "numtower, the numeric tower"], ["context/code/kb-tool.md", "code/kb-tool", "The kb tool"], ["context/code/tooling.md", "code/tooling", "The repository's own checkers"], ], @@ -82,6 +83,7 @@ const SECTIONS = [ ["context/test-evidence/tui.md", "test-evidence/tui", "The TUI"], ["context/test-evidence/cli.md", "test-evidence/cli", "The headless CLI"], ["context/test-evidence/gp.md", "test-evidence/gp", "gp-tree and improvise-gp"], + ["context/test-evidence/numtower.md", "test-evidence/numtower", "numtower"], ["context/test-evidence/tooling.md", "test-evidence/tooling", "The repository's own gates"], ["context/test-evidence/layout-cost.md", "test-evidence/layout-cost", "Layout and keypress cost"], ["context/test-evidence/scale-bisect.md", "test-evidence/scale-bisect", "Scale bisect"],