predict: envelope (interval) arithmetic for derived contract values (DBU-628)

packages/fixed_math/sources/interval.move

@@ -0,0 +1,110 @@
+// Copyright (c) Mysten Labs, Inc.
+// SPDX-License-Identifier: Apache-2.0
+
+/// Sound interval arithmetic over non-negative 1e9-scaled fixed-point values.
+///
+/// An `Interval` carries a lower and upper bound on a derived quantity; the
+/// soundness contract is that whenever every input's true value lies inside its
+/// interval, the output's true value lies inside the output interval. Bounds are
+/// maintained with directed rounding (`math::mul`/`math::mul_up`,
+/// `math::div`/`math::div_up`), so each side absorbs its own rounding and no
+/// separate error tracking is needed. Exact protocol atoms (quantities, floors,
+/// balances) enter as zero-width intervals; consumers collapse to a scalar by
+/// reading one side (`lo` for protocol outflows, `hi` for inflows, liabilities,
+/// and reserves) at explicit, auditable call sites.
+///
+/// `Interval` deliberately has no `store` ability: envelopes are transaction-
+/// local by construction, and persistent state stores exact scalars only.
+///
+/// An operation aborts when its upper bound cannot be represented (u64 overflow
+/// on the hi side, or a subtraction whose result is definitely negative); the
+/// lower bound alone never aborts — underflow on the lo side clamps to zero,
+/// which is sound because every modeled quantity is non-negative.
+module fixed_math::interval;
+
+use fixed_math::math;
+
+const EInvalidBounds: u64 = 0;
+const EDefinitelyNegative: u64 = 1;
+
+/// Non-negative bounds with `lo <= hi`; width `hi - lo` is the carried error.
+public struct Interval has copy, drop {
+    lo: u64,
+    hi: u64,
+}
+
+/// Create an interval from validated bounds.
+public fun new(lo: u64, hi: u64): Interval {
+    assert!(lo <= hi, EInvalidBounds);
+    Interval { lo, hi }
+}
+
+/// Embed an exact value as a zero-width interval.
+public fun exact(value: u64): Interval {
+    Interval { lo: value, hi: value }
+}
+
+/// Lower bound: the collapse side for protocol outflows.
+public fun lo(self: &Interval): u64 {
+    self.lo
+}
+
+/// Upper bound: the collapse side for protocol inflows, liabilities, reserves.
+public fun hi(self: &Interval): u64 {
+    self.hi
+}
+
+/// Carried error: `hi - lo`.
+public fun width(self: &Interval): u64 {
+    self.hi - self.lo
+}
+
+/// Widen both sides by an absolute error bound; the lower side clamps at zero.
+/// Used to fold a certified evaluation-error constant into a computed value.
+public fun widen(self: &Interval, err: u64): Interval {
+    let lo = if (self.lo > err) self.lo - err else 0;
+    Interval { lo, hi: self.hi + err }
+}
+
+/// Sum: widths add.
+public fun add(self: &Interval, other: &Interval): Interval {
+    Interval { lo: self.lo + other.lo, hi: self.hi + other.hi }
+}
+
+/// Difference: `[lo - other.hi, hi - other.lo]`; widths add. The lower side
+/// clamps at zero (the true value may still be positive when the bounds cross);
+/// aborts only when the difference is definitely negative — for a quantity known
+/// non-negative that means broken inputs, not dust.
+public fun sub(self: &Interval, other: &Interval): Interval {
+    assert!(self.hi >= other.lo, EDefinitelyNegative);
+    let lo = if (self.lo > other.hi) self.lo - other.hi else 0;
+    Interval { lo, hi: self.hi - other.lo }
+}
+
+/// Scaled product: lo side rounds down, hi side rounds up (+1 ulp width).
+public fun mul(self: &Interval, other: &Interval): Interval {
+    Interval {
+        lo: math::mul(self.lo, other.lo),
+        hi: math::mul_up(self.hi, other.hi),
+    }
+}
+
+/// Scaled quotient: divides the lo side by the divisor's hi and the hi side by
+/// the divisor's lo. A divisor whose lower bound is zero might be zero and
+/// aborts (`math::EInputZero`).
+public fun div(self: &Interval, other: &Interval): Interval {
+    Interval {
+        lo: math::div(self.lo, other.hi),
+        hi: math::div_up(self.hi, other.lo),
+    }
+}
+
+/// Pointwise minimum: sound because min is monotone in both arguments.
+public fun min(self: &Interval, other: &Interval): Interval {
+    Interval { lo: self.lo.min(other.lo), hi: self.hi.min(other.hi) }
+}
+
+/// Pointwise maximum: sound because max is monotone in both arguments.
+public fun max(self: &Interval, other: &Interval): Interval {
+    Interval { lo: self.lo.max(other.lo), hi: self.hi.max(other.hi) }
+}

packages/fixed_math/sources/math.move

@@ -68,19 +68,40 @@ const INV_13_U128: u128 = 76_923_077;
 /// Returns the fixed-point scale: `500_000_000` is 50% and `1_000_000_000` is 100%.
 public macro fun float_scaling(): u64 { 1_000_000_000 }
 
+/// `|x|` beyond which `normal_cdf`/`normal_pdf` saturate to their exact tail
+/// limits. Callers that clamp a `normal_cdf` argument share this bound so the
+/// clamp cannot drift from the saturation threshold.
+public macro fun normal_cdf_saturation(): u64 { 8 * float_scaling!() }
+
 // === Public Functions ===
 
 /// Multiply two 1e9-scaled fixed-point values, rounding down.
 public fun mul(x: u64, y: u64): u64 {
     (((x as u128) * (y as u128)) / F) as u64
 }
 
+/// Multiply two 1e9-scaled fixed-point values, rounding up. `x*y + F - 1` cannot
+/// overflow u128 for any u64 operands, so no intermediate guard is needed.
+/// Directed counterpart of `mul`; interval arithmetic requires both sides.
+public fun mul_up(x: u64, y: u64): u64 {
+    (((x as u128) * (y as u128) + (F - 1)) / F) as u64
+}
+
 /// Divides two 1e9-scaled fixed-point values, rounding down.
 /// Aborts when `y` is zero or the quotient does not fit in `u64`.
 public fun div(x: u64, y: u64): u64 {
     (((x as u128) * F) / (y as u128)) as u64
 }
 
+/// Divides two 1e9-scaled fixed-point values, rounding up.
+/// Aborts when `y` is zero or the quotient does not fit in `u64`. `x*F + y - 1`
+/// cannot overflow u128 for any u64 operands, so no intermediate guard is needed.
+/// Directed counterpart of `div`; interval arithmetic requires both sides.
+public fun div_up(x: u64, y: u64): u64 {
+    assert!(y > 0, EInputZero);
+    (((x as u128) * F + (y as u128) - 1) / (y as u128)) as u64
+}
+
 /// Multiplies two raw integers, divides by a raw denominator, and rounds down.
 /// Aborts when the denominator is zero or the quotient does not fit in `u64`.
 public fun mul_div_down(x: u64, y: u64, denominator: u64): u64 {
@@ -153,7 +174,7 @@ public fun exp(x: &i64::I64): u64 {
 public fun normal_cdf(x: &i64::I64): u64 {
     let x_mag = x.magnitude();
     let x_negative = x.is_negative();
-    if (x_mag > 8 * float_scaling!()) {
+    if (x_mag > normal_cdf_saturation!()) {
         return if (x_negative) { 0 } else { float_scaling!() }
     };
     (normal_cdf_u128((x_mag as u128), x_negative) as u64)
@@ -163,7 +184,7 @@ public fun normal_cdf(x: &i64::I64): u64 {
 /// Returns a 1e9-scaled density with absolute error at most 50 raw units; tails beyond `|8|` round to zero.
 public fun normal_pdf(x: &i64::I64): u64 {
     let x_mag = x.magnitude();
-    if (x_mag > 8 * float_scaling!()) return 0;
+    if (x_mag > normal_cdf_saturation!()) return 0;
 
     let x_sq_half = (((x_mag as u128) * (x_mag as u128)) / (2 * F)) as u64;
     let exponent = i64::from_parts(x_sq_half, true);
@@ -179,6 +200,13 @@ public fun sqrt(x: u64, precision: u64): u64 {
     (sqrt_u128(scaled) / multiplier) as u64
 }
 
+/// Integer square root (floor) of a `u128`. For a value scaled by `1e18` the
+/// result is its `1e9`-scaled square root; the pricing variance path relies on
+/// this to keep `sqrt(w)` precise where `w` is only a few raw units at `1e9`.
+public fun isqrt(x: u128): u128 {
+    sqrt_u128(x)
+}
+
 /// 10^n for small non-negative n. Capped at 18 because 10^19 overflows u64.
 public fun pow10(n: u64): u64 {
     assert!(n <= 18, EPow10ExponentTooLarge);

packages/fixed_math/tests/interval/interval_tests.move

@@ -0,0 +1,143 @@
+// Copyright (c) Mysten Labs, Inc.
+// SPDX-License-Identifier: Apache-2.0
+
+#[test_only]
+module fixed_math::interval_tests;
+
+use fixed_math::{interval, math::{Self, float_scaling as float}};
+use std::unit_test::assert_eq;
+
+// Expected values are hand-derived (comments show the work) — never computed
+// through the module under test (unit-tests rule 1). The algebra's soundness
+// laws (truth-containment under directed rounding, width addition, abort
+// semantics) were validated against an exact-rational reference implementation
+// by property fuzzing before this module was written; the cases here pin the
+// Move implementation at exact hand-checkable points.
+
+#[test]
+fun exact_has_zero_width() {
+    let x = interval::exact(5);
+    assert_eq!(x.lo(), 5);
+    assert_eq!(x.hi(), 5);
+    assert_eq!(x.width(), 0);
+}
+
+#[test]
+fun new_keeps_bounds() {
+    let x = interval::new(2, 3);
+    assert_eq!(x.lo(), 2);
+    assert_eq!(x.hi(), 3);
+    assert_eq!(x.width(), 1);
+}
+
+#[test, expected_failure(abort_code = interval::EInvalidBounds)]
+fun new_inverted_bounds_aborts() {
+    interval::new(3, 2);
+    abort 999
+}
+
+#[test]
+fun add_sums_bounds_and_widths() {
+    // [2,3] + [10,20] = [12,23]; width 1 + 10 = 11
+    let s = interval::new(2, 3).add(&interval::new(10, 20));
+    assert_eq!(s.lo(), 12);
+    assert_eq!(s.hi(), 23);
+    assert_eq!(s.width(), 11);
+}
+
+#[test]
+fun sub_crosses_bounds_and_widths_add() {
+    // [10,20] - [2,3] = [10-3, 20-2] = [7,18]; width 10 + 1 = 11
+    let d = interval::new(10, 20).sub(&interval::new(2, 3));
+    assert_eq!(d.lo(), 7);
+    assert_eq!(d.hi(), 18);
+    assert_eq!(d.width(), 11);
+}
+
+#[test]
+fun sub_of_self_clamps_lo_and_keeps_width() {
+    // X - X where X = [2,3]: true difference is 0, but the algebra cannot see
+    // correlation: [2-3 -> clamp 0, 3-2] = [0,1] (the dependency-problem witness)
+    let x = interval::new(2, 3);
+    let d = x.sub(&x);
+    assert_eq!(d.lo(), 0);
+    assert_eq!(d.hi(), 1);
+}
+
+#[test, expected_failure(abort_code = interval::EDefinitelyNegative)]
+fun sub_definitely_negative_aborts() {
+    // hi(3) < other.lo(10): every point of the difference is negative
+    interval::new(2, 3).sub(&interval::new(10, 20));
+    abort 999
+}
+
+#[test]
+fun mul_directs_rounding_outward() {
+    // [2,3] * [5,7] raw: lo = floor(2*5/1e9) = 0, hi = ceil(3*7/1e9) = 1
+    let p = interval::new(2, 3).mul(&interval::new(5, 7));
+    assert_eq!(p.lo(), 0);
+    assert_eq!(p.hi(), 1);
+}
+
+#[test]
+fun mul_exact_operands_stay_tight() {
+    // [1e9,2e9] * [3e9,4e9] = [3e9, 8e9] exactly (1.0*3.0=3.0, 2.0*4.0=8.0)
+    let p = interval::new(float!(), 2 * float!()).mul(&interval::new(3 * float!(), 4 * float!()));
+    assert_eq!(p.lo(), 3 * float!());
+    assert_eq!(p.hi(), 8 * float!());
+}
+
+#[test]
+fun mul_of_one_raw_unit_carries_one_ulp() {
+    // exact(1) * exact(1): true product 1e-9 raw; lo floors to 0, hi ceils to 1
+    let p = interval::exact(1).mul(&interval::exact(1));
+    assert_eq!(p.lo(), 0);
+    assert_eq!(p.hi(), 1);
+}
+
+#[test]
+fun div_crosses_divisor_bounds() {
+    // [6e9,8e9] / [2e9,4e9]: lo = floor(6/4 * 1e9) = 1.5e9, hi = ceil(8/2 * 1e9) = 4e9
+    let q = interval::new(6 * float!(), 8 * float!()).div(&interval::new(2 * float!(), 4 * float!()));
+    assert_eq!(q.lo(), 1_500_000_000);
+    assert_eq!(q.hi(), 4 * float!());
+}
+
+#[test, expected_failure(abort_code = math::EInputZero)]
+fun div_by_possibly_zero_divisor_aborts() {
+    // divisor lo == 0: the divisor might be zero, hi side must abort
+    interval::new(1, 2).div(&interval::new(0, float!()));
+    abort 999
+}
+
+#[test]
+fun min_is_pointwise() {
+    // min([2,10],[5,7]) = [2,7]
+    let m = interval::new(2, 10).min(&interval::new(5, 7));
+    assert_eq!(m.lo(), 2);
+    assert_eq!(m.hi(), 7);
+}
+
+#[test]
+fun max_is_pointwise() {
+    // max([2,10],[5,7]) = [5,10]
+    let m = interval::new(2, 10).max(&interval::new(5, 7));
+    assert_eq!(m.lo(), 5);
+    assert_eq!(m.hi(), 10);
+}
+
+#[test]
+fun widen_extends_both_sides() {
+    // [5,10] widened by 3 = [2,13]
+    let w = interval::new(5, 10).widen(3);
+    assert_eq!(w.lo(), 2);
+    assert_eq!(w.hi(), 13);
+}
+
+#[test]
+fun widen_clamps_lo_at_zero() {
+    // [1,10] widened by 5: lo clamps to 0, hi = 15
+    let w = interval::new(1, 10).widen(5);
+    assert_eq!(w.lo(), 0);
+    assert_eq!(w.hi(), 15);
+}

packages/fixed_math/tests/math/math_tests.move

@@ -108,6 +108,58 @@ fun div_rounds_down_to_integer_unit() {
     assert_eq!(math::div(1, 3), 333_333_333);
 }
 
+#[test]
+fun mul_up_ceils_one_raw_unit_of_dust() {
+    // ceil(1 * 1 / 1e9) = 1; mul floors the same product to 0
+    assert_eq!(math::mul_up(1, 1), 1);
+}
+
+#[test]
+fun mul_up_exact_product_does_not_round() {
+    // 1.5 * 2.25 = 3.375 exactly at 1e9 scale
+    assert_eq!(math::mul_up(float!() + float!() / 2, 2 * float!() + float!() / 4), 3_375_000_000);
+}
+
+#[test]
+fun mul_up_rounds_half_up_to_next_unit() {
+    // 7 * 0.5 = 3.5 raw units: ceil = 4, floor = 3
+    assert_eq!(math::mul_up(7, float!() / 2), 4);
+}
+
+#[test]
+fun mul_up_of_zero_is_zero() {
+    assert_eq!(math::mul_up(0, float!()), 0);
+}
+
+#[test]
+fun div_up_ceils_half_unit_quotient() {
+    // ceil(1 * 1e9 / 2e9) = ceil(0.5) = 1; div floors it to 0
+    assert_eq!(math::div_up(1, 2 * float!()), 1);
+}
+
+#[test]
+fun div_up_exact_division_does_not_round() {
+    // 5 / 2 = 2.5e9 exactly
+    assert_eq!(math::div_up(5 * float!(), 2 * float!()), 2_500_000_000);
+}
+
+#[test]
+fun div_up_rounds_repeating_quotient_up() {
+    // ceil(1 * 1e9 / 3) = 333_333_334; div floors to 333_333_333
+    assert_eq!(math::div_up(1, 3), 333_333_334);
+}
+
+#[test]
+fun div_up_of_zero_is_zero() {
+    assert_eq!(math::div_up(0, float!()), 0);
+}
+
+#[test, expected_failure(abort_code = math::EInputZero)]
+fun div_up_zero_divisor_aborts() {
+    math::div_up(float!(), 0);
+    abort 999
+}
+
 #[test]
 fun mul_div_down_floors_raw_integer_ratio() {
     // floor(10 * 10 / 6) = 16