the-q 0.2.1

//! Error-propagation (Lipschitz) lemmas (V7), the layer under interval
//! arithmetic and the accumulation bounds.
//!
//! Stated division-free through `frac_diff_le`. `add`/`sub` have constant 1;
//! `mul` is bounded by `|a|·|b - b'| + |b'|·|a - a'|`; `recip` costs
//! `md²/mn²` for a denominator bounded below by `mn/md`; `div` composes the
//! two. For `mul` and `div` the bare identities (`lemma_mul_lipschitz`,
//! `lemma_div_lipschitz`) are kept beside the bounds because callers use both.

use verus_builtin_macros::verus;

#[allow(unused_imports)]
use verus_builtin::*;
#[allow(unused_imports)]
use vstd::prelude::*;

#[allow(unused_imports)]
use crate::model::*;
#[allow(unused_imports)]
use crate::types::Rat;

verus! {

/// `|n1/d1 - n2/d2| <= en/ed`, written without division.
///
/// Multiplying through by `d1·d2·ed` (all positive) turns the statement into
/// `|n1·d2 - n2·d1| · ed <= en · d1 · d2`.
pub open spec fn frac_diff_le(n1: int, d1: int, n2: int, d2: int, en: int, ed: int) -> bool {
    abs_int(n1 * d2 - n2 * d1) * ed <= en * (d1 * d2)
}

/// `|u ± v| <= |u| + |v|`.
pub proof fn lemma_abs_triangle(u: int, v: int)
    ensures
        abs_int(u + v) <= abs_int(u) + abs_int(v),
        abs_int(u - v) <= abs_int(u) + abs_int(v),
{
}

/// `|u · v| == |u| · |v|` with both signs unknown.
pub proof fn lemma_abs_prod(u: int, v: int)
    ensures
        abs_int(u * v) == abs_int(u) * abs_int(v),
{
    if v > 0 {
        lemma_abs_mul_pos(u, v);
        assert(abs_int(v) == v);
    } else if v == 0 {
        assert(u * v == 0) by (nonlinear_arith)
            requires
                v == 0,
        ;
        assert(abs_int(u) * 0 == 0) by (nonlinear_arith);
    } else {
        // Fold the sign of `v` out and reuse the positive case. Every rewrite
        // is explicit. The solver moves from `|u·(-v)|` to `|u·v|` one step at
        // a time, because each bridge is an equality *under* an `abs_int`
        // rather than a bare arithmetic step.
        lemma_abs_mul_pos(u, -v);
        assert(abs_int(u * (-v)) == abs_int(u) * (-v));
        assert(u * (-v) == -(u * v)) by (nonlinear_arith);
        assert(abs_int(u * (-v)) == abs_int(-(u * v)));
        assert(abs_int(-(u * v)) == abs_int(u * v));
        assert(abs_int(v) == -v);
        assert(abs_int(u) * (-v) == abs_int(u) * abs_int(v));
    }
}

/// Scaling a `frac_diff_le` by a positive constant. `en/ed == (en·c)/(ed·c)`.
///
/// This lemma applies whenever two bounds carry different error denominators.
/// It brings them onto a common denominator before they combine.
pub proof fn lemma_frac_diff_scale(
    n1: int,
    d1: int,
    n2: int,
    d2: int,
    en: int,
    ed: int,
    c: int,
)
    requires
        c > 0,
        frac_diff_le(n1, d1, n2, d2, en, ed),
    ensures
        frac_diff_le(n1, d1, n2, d2, en * c, ed * c),
{
    let a = abs_int(n1 * d2 - n2 * d1);
    lemma_mul_le_mono(c, a * ed, en * (d1 * d2));
    assert(c * (a * ed) == a * (ed * c)) by (nonlinear_arith);
    assert(c * (en * (d1 * d2)) == (en * c) * (d1 * d2)) by (nonlinear_arith);
}

/// Multiplication by a non-negative constant preserves `<=`.
///
/// This is a one-line fact. The bounds below apply it a dozen times, and the
/// multiplier is a different five-factor product each time. A named lemma keeps
/// each application a rewrite rather than another nonlinear goal.
pub proof fn lemma_mul_le_mono(u: int, v: int, w: int)
    requires
        u >= 0,
        v <= w,
    ensures
        u * v <= u * w,
{
    assert(u * v <= u * w) by (nonlinear_arith)
        requires
            u >= 0,
            v <= w,
    ;
}

/// Two `<=` between non-negatives multiply.
pub proof fn lemma_mul_le_mono2(a: int, b: int, c: int, d: int)
    requires
        0 <= a <= c,
        0 <= b <= d,
    ensures
        a * b <= c * d,
{
    lemma_mul_le_mono(a, b, d);
    lemma_mul_le_mono(d, a, c);
    assert(a * d == d * a) by (nonlinear_arith);
    assert(d * c == c * d) by (nonlinear_arith);
}

/// **Addition is 1-Lipschitz in each argument.**
///
/// If `|a - a'| <= e1` and `|b - b'| <= e2` then `|(a+b) - (a'+b')| <= e1 + e2`.
pub proof fn lemma_add_lipschitz(
    an: int,
    ad: int,
    bn: int,
    bd: int,
    a2n: int,
    a2d: int,
    b2n: int,
    b2d: int,
    e1: int,
    e2: int,
    ed: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
        ed > 0,
        e1 >= 0,
        e2 >= 0,
        frac_diff_le(an, ad, a2n, a2d, e1, ed),
        frac_diff_le(bn, bd, b2n, b2d, e2, ed),
    ensures
        frac_diff_le(
            an * bd + bn * ad,
            ad * bd,
            a2n * b2d + b2n * a2d,
            a2d * b2d,
            e1 + e2,
            ed,
        ),
{
    // (a+b) - (a'+b') == (a - a') + (b - b'). Also |x + y| <= |x| + |y|.
    lemma_triangle(an, ad, a2n, a2d, bn, bd, b2n, b2d, e1, e2, ed);
}

/// The triangle inequality, in the cross-multiplied form the addition lemma
/// needs.
#[verifier::rlimit(20)]
pub proof fn lemma_triangle(
    an: int,
    ad: int,
    a2n: int,
    a2d: int,
    bn: int,
    bd: int,
    b2n: int,
    b2d: int,
    e1: int,
    e2: int,
    ed: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
        ed > 0,
        e1 >= 0,
        e2 >= 0,
        abs_int(an * a2d - a2n * ad) * ed <= e1 * (ad * a2d),
        abs_int(bn * b2d - b2n * bd) * ed <= e2 * (bd * b2d),
    ensures
        abs_int((an * bd + bn * ad) * (a2d * b2d) - (a2n * b2d + b2n * a2d) * (ad * bd)) * ed
            <= (e1 + e2) * ((ad * bd) * (a2d * b2d)),
{
    // The difference factors as
    //   (a - a')·(bd·b2d) + (b - b')·(ad·a2d)
    // with positive weights. Apply |x+y| <= |x|+|y| and multiply the two
    // hypotheses by those weights.
    //
    // The factorisation goes to the solver in four small ring steps. Handed
    // over whole it exhausts the resource limit. Eight variables and degree
    // four exceed what the nonlinear tactic handles in one goal.
    assert((an * bd + bn * ad) * (a2d * b2d) == (an * bd) * (a2d * b2d) + (bn * ad) * (a2d * b2d))
        by (nonlinear_arith);
    assert((an * bd) * (a2d * b2d) == (an * a2d) * (bd * b2d)) by (nonlinear_arith);
    assert((bn * ad) * (a2d * b2d) == (bn * b2d) * (ad * a2d)) by (nonlinear_arith);
    assert((a2n * b2d + b2n * a2d) * (ad * bd) == (a2n * b2d) * (ad * bd) + (b2n * a2d) * (ad
        * bd)) by (nonlinear_arith);
    assert((a2n * b2d) * (ad * bd) == (a2n * ad) * (bd * b2d)) by (nonlinear_arith);
    assert((b2n * a2d) * (ad * bd) == (b2n * bd) * (ad * a2d)) by (nonlinear_arith);
    assert((an * a2d - a2n * ad) * (bd * b2d) == (an * a2d) * (bd * b2d) - (a2n * ad) * (bd * b2d))
        by (nonlinear_arith);
    assert((bn * b2d - b2n * bd) * (ad * a2d) == (bn * b2d) * (ad * a2d) - (b2n * bd) * (ad * a2d))
        by (nonlinear_arith);
    assert((an * bd + bn * ad) * (a2d * b2d) - (a2n * b2d + b2n * a2d) * (ad * bd) == (an * a2d
        - a2n * ad) * (bd * b2d) + (bn * b2d - b2n * bd) * (ad * a2d));
    assert(abs_int(
        (an * a2d - a2n * ad) * (bd * b2d) + (bn * b2d - b2n * bd) * (ad * a2d),
    ) <= abs_int(an * a2d - a2n * ad) * (bd * b2d) + abs_int(bn * b2d - b2n * bd) * (ad * a2d))
        by (nonlinear_arith)
        requires
            bd > 0,
            b2d > 0,
            ad > 0,
            a2d > 0,
    ;
    assert(abs_int(an * a2d - a2n * ad) * (bd * b2d) * ed <= e1 * (ad * a2d) * (bd * b2d))
        by (nonlinear_arith)
        requires
            bd > 0,
            b2d > 0,
            abs_int(an * a2d - a2n * ad) * ed <= e1 * (ad * a2d),
    ;
    assert(abs_int(bn * b2d - b2n * bd) * (ad * a2d) * ed <= e2 * (bd * b2d) * (ad * a2d))
        by (nonlinear_arith)
        requires
            ad > 0,
            a2d > 0,
            abs_int(bn * b2d - b2n * bd) * ed <= e2 * (bd * b2d),
    ;
    // Assemble the result. Scale the triangle inequality by `ed`, then chain
    // the two weighted hypotheses. Each step names its operands, thus the
    // blocks stay small enough to discharge.
    let dd = (an * bd + bn * ad) * (a2d * b2d) - (a2n * b2d + b2n * a2d) * (ad * bd);
    let t = abs_int(an * a2d - a2n * ad) * (bd * b2d) + abs_int(bn * b2d - b2n * bd) * (ad * a2d);
    assert(abs_int(dd) <= t);
    assert(abs_int(dd) * ed <= t * ed) by (nonlinear_arith)
        requires
            ed > 0,
            abs_int(dd) <= t,
    ;
    assert(t * ed == abs_int(an * a2d - a2n * ad) * (bd * b2d) * ed + abs_int(
        bn * b2d - b2n * bd,
    ) * (ad * a2d) * ed) by (nonlinear_arith)
        requires
            t == abs_int(an * a2d - a2n * ad) * (bd * b2d) + abs_int(bn * b2d - b2n * bd) * (ad
                * a2d),
    ;
    assert(e1 * (ad * a2d) * (bd * b2d) + e2 * (bd * b2d) * (ad * a2d) == (e1 + e2) * ((ad * bd) * (
    a2d * b2d))) by (nonlinear_arith);
}

/// **Multiplication on a bounded domain.**
///
/// `|a·b - a'·b'| <= |a|·|b - b'| + |b'|·|a - a'|`. Every opinion component
/// lies on `[0, 1]`. There both coefficients are at most `1`, thus the errors
/// add.
#[verifier::rlimit(20)]
pub proof fn lemma_mul_lipschitz(
    an: int,
    ad: int,
    bn: int,
    bd: int,
    a2n: int,
    a2d: int,
    b2n: int,
    b2d: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
    ensures
        (an * bn) * (a2d * b2d) - (a2n * b2n) * (ad * bd) == (an * a2d) * (bn * b2d - b2n * bd) + (
        b2n * bd) * (an * a2d - a2n * ad),
{
    // Distribute both products, then match term by term. The two middle terms
    // are the same monomial and cancel.
    assert((an * a2d) * (bn * b2d - b2n * bd) == (an * a2d) * (bn * b2d) - (an * a2d) * (b2n * bd))
        by (nonlinear_arith);
    assert((b2n * bd) * (an * a2d - a2n * ad) == (b2n * bd) * (an * a2d) - (b2n * bd) * (a2n * ad))
        by (nonlinear_arith);
    assert((an * a2d) * (b2n * bd) == (b2n * bd) * (an * a2d)) by (nonlinear_arith);
    assert((an * bn) * (a2d * b2d) == (an * a2d) * (bn * b2d)) by (nonlinear_arith);
    assert((a2n * b2n) * (ad * bd) == (b2n * bd) * (a2n * ad)) by (nonlinear_arith);
}

/// Product bound: `|a - a'| <= e1/ed`, `|b - b'| <= e2/ed`, `|a| <= ca`,
/// `|b'| <= cb` give `|a·b - a'·b'| <= (ca·e2 + cb·e1)/ed`. The magnitude
/// bounds sit on the diagonal [`lemma_mul_lipschitz`] produces (`a` unprimed,
/// `b` primed). [`lemma_mul_lipschitz_unit`] is the `ca = cb = 1` case.
#[verifier::rlimit(40)]
pub proof fn lemma_mul_lipschitz_bound(
    an: int,
    ad: int,
    bn: int,
    bd: int,
    a2n: int,
    a2d: int,
    b2n: int,
    b2d: int,
    e1: int,
    e2: int,
    ed: int,
    ca: int,
    cb: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
        ed > 0,
        e1 >= 0,
        e2 >= 0,
        ca >= 0,
        cb >= 0,
        abs_int(an) <= ca * ad,
        abs_int(b2n) <= cb * b2d,
        frac_diff_le(an, ad, a2n, a2d, e1, ed),
        frac_diff_le(bn, bd, b2n, b2d, e2, ed),
    ensures
        frac_diff_le(an * bn, ad * bd, a2n * b2n, a2d * b2d, ca * e2 + cb * e1, ed),
{
    // Name every absolute value first: inline, `abs_int` unfolds to an
    // if-then-else under every product and the solver exceeds its limit.
    let x = an * a2d - a2n * ad;
    let y = bn * b2d - b2n * bd;
    let p = (an * bn) * (a2d * b2d) - (a2n * b2n) * (ad * bd);
    let ax = abs_int(x);
    let ay = abs_int(y);
    let ap = abs_int(p);
    let pa = abs_int(an);
    let pb = abs_int(b2n);

    // Step 1: the identity, unchanged.
    lemma_mul_lipschitz(an, ad, bn, bd, a2n, a2d, b2n, b2d);
    assert(p == (an * a2d) * y + (b2n * bd) * x);

    // Step 2: |p| <= (pa·ay)·a2d + (pb·ax)·bd. Each product re-associates so
    // that the positive denominator is the outer factor. That is the shape
    // `lemma_abs_mul_pos` needs. The remaining unknown-sign pair goes through
    // `lemma_abs_prod`.
    lemma_abs_triangle((an * a2d) * y, (b2n * bd) * x);
    assert((an * a2d) * y == (an * y) * a2d) by (nonlinear_arith);
    assert((b2n * bd) * x == (b2n * x) * bd) by (nonlinear_arith);
    lemma_abs_mul_pos(an * y, a2d);
    lemma_abs_mul_pos(b2n * x, bd);
    lemma_abs_prod(an, y);
    lemma_abs_prod(b2n, x);
    assert(ap <= (pa * ay) * a2d + (pb * ax) * bd);

    // Step 3: scale by `ed` and regroup so that each hypothesis appears
    // verbatim as a factor. `ay·ed` and `ax·ed` are exactly the two
    // `frac_diff_le`s. The regrouping uses three small identities instead of
    // one wide one.
    lemma_mul_le_mono(ed, ap, (pa * ay) * a2d + (pb * ax) * bd);
    assert(ap * ed == ed * ap) by (nonlinear_arith);
    assert(ed * ((pa * ay) * a2d + (pb * ax) * bd) == ed * ((pa * ay) * a2d) + ed * ((pb * ax) * bd))
        by (nonlinear_arith);
    assert(ed * ((pa * ay) * a2d) == (pa * a2d) * (ay * ed)) by (nonlinear_arith);
    assert(ed * ((pb * ax) * bd) == (pb * bd) * (ax * ed)) by (nonlinear_arith);

    // Step 4: discharge the two hypotheses.
    assert(pa * a2d >= 0) by (nonlinear_arith)
        requires
            pa >= 0,
            a2d > 0,
    ;
    assert(pb * bd >= 0) by (nonlinear_arith)
        requires
            pb >= 0,
            bd > 0,
    ;
    lemma_mul_le_mono(pa * a2d, ay * ed, e2 * (bd * b2d));
    lemma_mul_le_mono(pb * bd, ax * ed, e1 * (ad * a2d));

    // Step 5: replace the magnitudes by their ceilings. Products appear in full:
    // a `let` binding is invisible inside `by (nonlinear_arith)`.
    assert(a2d * (e2 * (bd * b2d)) >= 0) by (nonlinear_arith)
        requires
            a2d > 0,
            e2 >= 0,
            bd > 0,
            b2d > 0,
    ;
    assert(bd * (e1 * (ad * a2d)) >= 0) by (nonlinear_arith)
        requires
            bd > 0,
            e1 >= 0,
            ad > 0,
            a2d > 0,
    ;
    lemma_mul_le_mono(a2d * (e2 * (bd * b2d)), pa, ca * ad);
    lemma_mul_le_mono(bd * (e1 * (ad * a2d)), pb, cb * b2d);
    assert((a2d * (e2 * (bd * b2d))) * pa == (pa * a2d) * (e2 * (bd * b2d))) by (nonlinear_arith);
    assert((bd * (e1 * (ad * a2d))) * pb == (pb * bd) * (e1 * (ad * a2d))) by (nonlinear_arith);
    assert((a2d * (e2 * (bd * b2d))) * (ca * ad) == (ca * e2) * ((ad * bd) * (a2d * b2d)))
        by (nonlinear_arith);
    assert((bd * (e1 * (ad * a2d))) * (cb * b2d) == (cb * e1) * ((ad * bd) * (a2d * b2d)))
        by (nonlinear_arith);

    // Step 6: add the two halves.
    assert((ca * e2) * ((ad * bd) * (a2d * b2d)) + (cb * e1) * ((ad * bd) * (a2d * b2d)) == (ca * e2
        + cb * e1) * ((ad * bd) * (a2d * b2d))) by (nonlinear_arith);
}

/// The product bound with `|a|, |b'| <= 1`: errors add.
pub proof fn lemma_mul_lipschitz_unit(
    an: int,
    ad: int,
    bn: int,
    bd: int,
    a2n: int,
    a2d: int,
    b2n: int,
    b2d: int,
    e1: int,
    e2: int,
    ed: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
        ed > 0,
        e1 >= 0,
        e2 >= 0,
        abs_int(an) <= ad,
        abs_int(b2n) <= b2d,
        frac_diff_le(an, ad, a2n, a2d, e1, ed),
        frac_diff_le(bn, bd, b2n, b2d, e2, ed),
    ensures
        frac_diff_le(an * bn, ad * bd, a2n * b2n, a2d * b2d, e1 + e2, ed),
{
    assert(1 * ad == ad);
    assert(1 * b2d == b2d);
    lemma_mul_lipschitz_bound(an, ad, bn, bd, a2n, a2d, b2n, b2d, e1, e2, ed, 1, 1);
    assert(1 * e2 + 1 * e1 == e1 + e2);
}

/// The algebraic core of the quotient bound:
/// `a/b - a'/b' == (b'·(a - a') + a·(b' - b)) / (b·b')`.
#[verifier::rlimit(20)]
pub proof fn lemma_div_lipschitz(
    an: int,
    ad: int,
    bn: int,
    bd: int,
    a2n: int,
    a2d: int,
    b2n: int,
    b2d: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
        bn != 0,
        b2n != 0,
    ensures
        // With A/B == (an·bd)/(ad·bn) and A'/B' == (a2n·b2d)/(a2d·b2n), the
        // cross-multiplied difference splits into the numerator perturbation
        // X == an·a2d - a2n·ad and the denominator perturbation
        // Y == bn·b2d - b2n·bd, weighted by bd·b2n and a2n·ad respectively.
        (an * bd) * (a2d * b2n) - (a2n * b2d) * (ad * bn) == (bd * b2n) * (an * a2d - a2n * ad) - (
        a2n * ad) * (bn * b2d - b2n * bd),
{
    // This has the same shape as the multiplication case. Distribute, and the
    // two `a2n·ad·b2n·bd` terms cancel.
    assert((bd * b2n) * (an * a2d - a2n * ad) == (bd * b2n) * (an * a2d) - (bd * b2n) * (a2n * ad))
        by (nonlinear_arith);
    assert((a2n * ad) * (bn * b2d - b2n * bd) == (a2n * ad) * (bn * b2d) - (a2n * ad) * (b2n * bd))
        by (nonlinear_arith);
    assert((bd * b2n) * (a2n * ad) == (a2n * ad) * (b2n * bd)) by (nonlinear_arith);
    assert((an * bd) * (a2d * b2n) == (bd * b2n) * (an * a2d)) by (nonlinear_arith);
    assert((a2n * b2d) * (ad * bn) == (a2n * ad) * (bn * b2d)) by (nonlinear_arith);
}

/// Reciprocal bound: for `b, b' >= mn/md > 0` (given as `mn·bd <= md·bn`) and
/// `|b - b'| <= e2/ed`, `|1/b - 1/b'| <= (e2·md²)/(ed·mn²)`.
pub proof fn lemma_recip_lipschitz_bound(
    bn: int,
    bd: int,
    b2n: int,
    b2d: int,
    e2: int,
    ed: int,
    mn: int,
    md: int,
)
    requires
        bd > 0,
        b2d > 0,
        ed > 0,
        bn > 0,
        b2n > 0,
        mn > 0,
        md > 0,
        e2 >= 0,
        mn * bd <= md * bn,
        mn * b2d <= md * b2n,
        frac_diff_le(bn, bd, b2n, b2d, e2, ed),
    ensures
        frac_diff_le(bd, bn, b2d, b2n, e2 * (md * md), ed * (mn * mn)),
{
    let y = bn * b2d - b2n * bd;
    let ay = abs_int(y);

    // Flipping both fractions negates the cross-difference, which `abs_int`
    // does not see.
    assert(bd * b2n - b2d * bn == -y) by (nonlinear_arith)
        requires
            y == bn * b2d - b2n * bd,
    ;
    assert(abs_int(-y) == ay);

    // Scale the hypothesis `ay·ed <= e2·(bd·b2d)` by `mn^2`. Place the two
    // factors of `mn` next to the two denominators they displace.
    assert(mn * mn >= 0) by (nonlinear_arith)
        requires
            mn > 0,
    ;
    lemma_mul_le_mono(mn * mn, ay * ed, e2 * (bd * b2d));
    assert((mn * mn) * (ay * ed) == ay * (ed * (mn * mn))) by (nonlinear_arith);
    assert((mn * mn) * (e2 * (bd * b2d)) == e2 * ((mn * bd) * (mn * b2d))) by (nonlinear_arith);

    // Spend both copies of the lower bound at once. This is the `md^2`.
    assert(mn * bd >= 0) by (nonlinear_arith)
        requires
            mn > 0,
            bd > 0,
    ;
    assert(mn * b2d >= 0) by (nonlinear_arith)
        requires
            mn > 0,
            b2d > 0,
    ;
    lemma_mul_le_mono2(mn * bd, mn * b2d, md * bn, md * b2n);
    lemma_mul_le_mono(e2, (mn * bd) * (mn * b2d), (md * bn) * (md * b2n));
    assert(e2 * ((md * bn) * (md * b2n)) == (e2 * (md * md)) * (bn * b2n)) by (nonlinear_arith);
}

/// Quotient bound: with `b, b' >= mn/md > 0` (division-free as `mn·bd <= md·bn`)
/// and `|a| <= ca`,
/// `|a/b - a'/b'| <= (ca·e2·md² + md·e1·mn²) / (ed·mn²)`.
/// Proved as [`lemma_recip_lipschitz_bound`] fed into
/// [`lemma_mul_lipschitz_bound`]. Divisors must be strictly positive.
pub proof fn lemma_div_lipschitz_bound(
    an: int,
    ad: int,
    bn: int,
    bd: int,
    a2n: int,
    a2d: int,
    b2n: int,
    b2d: int,
    e1: int,
    e2: int,
    ed: int,
    ca: int,
    mn: int,
    md: int,
)
    requires
        ad > 0,
        bd > 0,
        a2d > 0,
        b2d > 0,
        ed > 0,
        bn > 0,
        b2n > 0,
        mn > 0,
        md > 0,
        e1 >= 0,
        e2 >= 0,
        ca >= 0,
        mn * bd <= md * bn,
        mn * b2d <= md * b2n,
        abs_int(an) <= ca * ad,
        frac_diff_le(an, ad, a2n, a2d, e1, ed),
        frac_diff_le(bn, bd, b2n, b2d, e2, ed),
    ensures
        frac_diff_le(
            an * bd,
            ad * bn,
            a2n * b2d,
            a2d * b2n,
            ca * (e2 * (md * md)) + md * (e1 * (mn * mn)),
            ed * (mn * mn),
        ),
{
    // Put both errors over the common denominator `ed·mn^2`. The reciprocal
    // bound already arrives there, thus the dividend gets rescaled.
    assert(mn * mn > 0) by (nonlinear_arith)
        requires
            mn > 0,
    ;
    assert(ed * (mn * mn) > 0) by (nonlinear_arith)
        requires
            ed > 0,
            mn > 0,
    ;
    lemma_frac_diff_scale(an, ad, a2n, a2d, e1, ed, mn * mn);
    lemma_recip_lipschitz_bound(bn, bd, b2n, b2d, e2, ed, mn, md);

    // The product bound bounds `|1/b'|`, the primed side, not `|1/b|`.
    // `mn·b2d <= md·b2n` with `mn >= 1` gives `b2d <= md·b2n`, thus `md`
    // serves.
    assert(b2d <= mn * b2d) by (nonlinear_arith)
        requires
            mn >= 1,
            b2d > 0,
    ;
    assert(abs_int(b2d) == b2d);
    assert(abs_int(b2d) <= md * b2n);

    // `a/b == a · (1/b)`. The reciprocal is the second argument, thus its
    // denominator is `bn` and its numerator is `bd`.
    lemma_mul_lipschitz_bound(
        an,
        ad,
        bd,
        bn,
        a2n,
        a2d,
        b2d,
        b2n,
        e1 * (mn * mn),
        e2 * (md * md),
        ed * (mn * mn),
        ca,
        md,
    );
}

/// Triangle inequality on fractions: `|x - y| <= e1/E` and `|y - z| <= e2/E`
/// give `|x - z| <= (e1 + e2)/E`.
pub proof fn lemma_frac_triangle(
    xn: int,
    xd: int,
    yn: int,
    yd: int,
    zn: int,
    zd: int,
    e1: int,
    e2: int,
    ee: int,
)
    requires
        xd > 0,
        yd > 0,
        zd > 0,
        ee > 0,
        abs_int(xn * yd - yn * xd) * ee <= e1 * (xd * yd),
        abs_int(yn * zd - zn * yd) * ee <= e2 * (yd * zd),
    ensures
        abs_int(xn * zd - zn * xd) * ee <= (e1 + e2) * (xd * zd),
{
    assert((xn * zd - zn * xd) * yd == (xn * yd - yn * xd) * zd + (yn * zd - zn * yd) * xd)
        by (nonlinear_arith);
    assert(abs_int((xn * zd - zn * xd) * yd) == abs_int(xn * zd - zn * xd) * yd)
        by (nonlinear_arith)
        requires
            yd > 0,
    ;
    assert(abs_int((xn * yd - yn * xd) * zd + (yn * zd - zn * yd) * xd) <= abs_int(
        xn * yd - yn * xd,
    ) * zd + abs_int(yn * zd - zn * yd) * xd) by (nonlinear_arith)
        requires
            zd > 0,
            xd > 0,
    ;
    // Scale each hypothesis by the third denominator.
    assert((abs_int(xn * yd - yn * xd) * ee) * zd <= (e1 * (xd * yd)) * zd) by (nonlinear_arith)
        requires
            zd > 0,
            abs_int(xn * yd - yn * xd) * ee <= e1 * (xd * yd),
    ;
    assert((abs_int(yn * zd - zn * yd) * ee) * xd <= (e2 * (yd * zd)) * xd) by (nonlinear_arith)
        requires
            xd > 0,
            abs_int(yn * zd - zn * yd) * ee <= e2 * (yd * zd),
    ;
    assert(((e1 + e2) * (xd * zd)) * yd == (e1 * (xd * zd)) * yd + (e2 * (xd * zd)) * yd)
        by (nonlinear_arith);
    assert((e1 * (xd * yd)) * zd == (e1 * (xd * zd)) * yd) by (nonlinear_arith);
    assert((e2 * (yd * zd)) * xd == (e2 * (xd * zd)) * yd) by (nonlinear_arith);
    // The triangle inequality is on `|X|·yd`. The hypotheses are on `|A|·ee`
    // and `|B|·ee`. Scaling by `ee` joins them, and that step is nonlinear.
    assert(abs_int(xn * zd - zn * xd) * yd <= abs_int(xn * yd - yn * xd) * zd + abs_int(
        yn * zd - zn * yd,
    ) * xd);
    assert((abs_int(xn * zd - zn * xd) * yd) * ee <= (abs_int(xn * yd - yn * xd) * zd + abs_int(
        yn * zd - zn * yd,
    ) * xd) * ee) by (nonlinear_arith)
        requires
            ee > 0,
            abs_int(xn * zd - zn * xd) * yd <= abs_int(xn * yd - yn * xd) * zd + abs_int(
                yn * zd - zn * yd,
            ) * xd,
    ;
    assert((abs_int(xn * zd - zn * xd) * yd) * ee == (abs_int(xn * zd - zn * xd) * ee) * yd)
        by (nonlinear_arith);
    assert((abs_int(xn * yd - yn * xd) * zd + abs_int(yn * zd - zn * yd) * xd) * ee == (abs_int(
        xn * yd - yn * xd,
    ) * ee) * zd + (abs_int(yn * zd - zn * yd) * ee) * xd) by (nonlinear_arith);
    assert((abs_int(xn * zd - zn * xd) * ee) * yd <= ((e1 + e2) * (xd * zd)) * yd);
    assert(abs_int(xn * zd - zn * xd) * ee <= (e1 + e2) * (xd * zd)) by (nonlinear_arith)
        requires
            yd > 0,
            (abs_int(xn * zd - zn * xd) * ee) * yd <= ((e1 + e2) * (xd * zd)) * yd,
    ;
}

/// V8 induction step: one rounded `add` on an accumulator within `k` units is
/// within `k + 1`. This step's R3 unit (made absolute by `max(1, |v|) <= m`)
/// plus the carried error, which addition passes through unchanged.
pub proof fn lemma_abs_error_step(prev: Rat, pn: int, pd: int, next: Rat, r: Rat, k: nat, m: int)
    requires
        prev.wf(),
        next.wf(),
        r.wf(),
        pd > 0,
        m >= 1,
        within_abs_error(prev, pn, pd, k, m),
        within_error_bound(r, crate::q::add_n(prev, next), crate::q::prod_d(prev, next)),
        max_int(
            crate::q::prod_d(prev, next),
            abs_int(crate::q::add_n(prev, next)),
        ) <= m * crate::q::prod_d(prev, next),
    ensures
        within_abs_error(
            r,
            pn * next.d() + next.n() * pd,
            pd * next.d(),
            (k + 1) as nat,
            m,
        ),
{
    let ad = prev.d();
    let an = prev.n();
    let bd = next.d();
    let bn = next.n();
    let en = crate::q::add_n(prev, next);
    let ed = crate::q::prod_d(prev, next);
    let tn = pn * bd + bn * pd;
    let td = pd * bd;
    let e = pow2(precision_b());
    lemma_pow2_pos(precision_b());
    assert(ad > 0 && bd > 0 && ed == ad * bd && ed > 0) by (nonlinear_arith)
        requires
            ad > 0,
            bd > 0,
            ed == ad * bd,
    ;
    assert(pd * bd > 0) by (nonlinear_arith)
        requires
            pd > 0,
            bd > 0,
    ;
    assert(td > 0);
    // (a) this step's rounding error, in absolute form.
    assert(abs_int(r.n() * ed - en * r.d()) * e <= m * (r.d() * ed)) by (nonlinear_arith)
        requires
            abs_int(r.n() * ed - en * r.d()) * e <= r.d() * max_int(ed, abs_int(en)),
            max_int(ed, abs_int(en)) <= m * ed,
            r.d() > 0,
    ;
    // (b) the carried error. The `next` term cancels. The difference between
    // the step's exact value and the target is exactly bd^2 times the
    // accumulator's own error. Expanding en·td - tn·ed shows this. The two
    // `bn` terms are the same monomial and cancel, which leaves bd^2 times the
    // accumulator's own error.
    assert((an * bd + bn * ad) * (pd * bd) == (an * bd) * (pd * bd) + (bn * ad) * (pd * bd))
        by (nonlinear_arith);
    assert((pn * bd + bn * pd) * (ad * bd) == (pn * bd) * (ad * bd) + (bn * pd) * (ad * bd))
        by (nonlinear_arith);
    assert((bn * ad) * (pd * bd) == (bn * pd) * (ad * bd)) by (nonlinear_arith);
    assert((an * bd) * (pd * bd) == (bd * bd) * (an * pd)) by (nonlinear_arith);
    assert((pn * bd) * (ad * bd) == (bd * bd) * (pn * ad)) by (nonlinear_arith);
    assert((bd * bd) * (an * pd) - (bd * bd) * (pn * ad) == (bd * bd) * (an * pd - pn * ad))
        by (nonlinear_arith);
    assert(en * td - tn * ed == (bd * bd) * (an * pd - pn * ad));
    assert(ed * td == (bd * bd) * (ad * pd)) by (nonlinear_arith)
        requires
            ed == ad * bd,
            td == pd * bd,
    ;
    assert(abs_int((bd * bd) * (an * pd - pn * ad)) == (bd * bd) * abs_int(an * pd - pn * ad))
        by (nonlinear_arith)
        requires
            bd > 0,
    ;
    assert(abs_int(en * td - tn * ed) * e <= ((k as int) * m) * (ed * td)) by (nonlinear_arith)
        requires
            bd > 0,
            abs_int(an * pd - pn * ad) * e <= (k as int) * m * (ad * pd),
            abs_int(en * td - tn * ed) == (bd * bd) * abs_int(an * pd - pn * ad),
            ed * td == (bd * bd) * (ad * pd),
    ;
    lemma_frac_triangle(r.n(), r.d(), en, ed, tn, td, m, (k as int) * m, e);
    assert(m + (k as int) * m == ((k + 1) as int) * m) by (nonlinear_arith);
}

/// The relative-error unit this crate uses throughout, as its reciprocal
/// `2^61`. This form keeps the bounds division-free.
pub open spec fn unit_error() -> int {
    pow2(precision_b())
}

/// A single-step bound implies the `k`-step bound for `k >= 1`.
pub proof fn lemma_bound_1_implies_k(r: Rat, n: int, d: int, k: nat)
    requires
        d > 0,
        r.wf(),
        k >= 1,
        within_error_bound(r, n, d),
    ensures
        within_error_bound_k(r, n, d, k),
{
    assert(r.d() * crate::model::max_int(d, abs_int(n)) <= (k as int) * r.d()
        * crate::model::max_int(d, abs_int(n))) by (nonlinear_arith)
        requires
            k >= 1,
            r.d() > 0,
            d > 0,
    ;
}

} // verus!