diff --git a/V1-ROADMAP.md b/V1-ROADMAP.md index 3ad4139c..692f0414 100644 --- a/V1-ROADMAP.md +++ b/V1-ROADMAP.md @@ -1,6 +1,6 @@ # dashu Roadmap — v0.5.x and beyond -Last updated: 2026-07-11 +Last updated: 2026-07-20 Feature work deferred out of the **v0.5.0** release. The v0.5.x items are all **additive** (no breaking changes) and safe to ship as point releases on top of 0.5.0; the post-v1 items @@ -33,11 +33,22 @@ are longer-term goals. File:line references are anchors from the v0.5.0 tree and ### Correctness -- **Guaranteed-correct rounding (Ziv retry loop).** ✅ *Partially delivered.* `exp`, `exp_m1`, - `ln`, `ln_1p` are now guaranteed-correctly rounded via a Ziv retry loop in `dashu-float` - (`Context::ziv`, driven by the `ErrorBounds` preimage). Remaining: trig, hyperbolic, `powf`, - `hypot`, and inheriting the loop across `dashu-cmplx`'s complex transcendentals (which currently - route through the now-Ziv-backed real primitives but aren't themselves Ziv-wrapped). +- **Guaranteed-correct rounding (Ziv retry loop).** ✅ *Delivered for all real and complex + transcendentals.* `exp`, `exp_m1`, `ln`, `ln_1p`, the trigonometric family (`sin`, `cos`, + `sin_cos`, `tan`, `asin`, `acos`, `atan`, `atan2`), the hyperbolic family (`sinh`, `cosh`, + `sinh_cosh`, `tanh`, `asinh`, `acosh`, `atanh`), `hypot`, `powi`, and `powf` are + guaranteed-correctly rounded via a Ziv retry loop (`Context::ziv`/`ziv_pair`, driven by the + `ErrorBounds` preimage). Series transcendentals (trig, atan) use near-correct `_compute` cores the + wrapper certifies; composition transcendentals treat the Ziv-correct primitives as black boxes. + `powf`'s `exp(y·ln x)` amplification is handled by a work-precision radius + (`result.ulp()·(|y·ln x|+1)·(B+8)`) that shrinks as `B^{-guard}`, so the containment test + converges. `powi` (repeated squaring) carries a radius that scales with the squaring-compounding + error (`2^bit_len(n)·ulp`) and drops to zero when the chain is exact — required under the directed + rounding modes (`FBig`/`CBig` default to `mode::Zero`), where an exactly-representable result sits + on a one-sided rounding boundary. `dashu-cmplx`'s complex transcendentals + (`exp`/`log`/`powf`/`powi`/`sin`/`cos`/`tan`/`sin_cos`/`asin`/`acos`/`atan`/`sqrt`) are + correctly rounded via their own complex Ziv driver (certifies both parts); `asin`/`acos` use the + factored `1-z²=(1-z)(1+z)` to stay accurate right up to `z = ±1`. - **Signed-zero preservation in `CBig` zero short-circuits.** `sin_cos` and `sqr` take a fast path on exactly-zero input that returns `+0` components, so several Annex-G / IEEE signed-zero cases are not preserved (all numerically equal to `+0`, hence deferred): @@ -58,8 +69,8 @@ Consolidated from the original `CBig` design. All additive. complex-valued functions themselves are deferred.) - **More transcendentals** — complex `fma` (fused multiply-add — hard to round correctly), `rootofunity`, complex `agm`, and `exp2`/`exp10`/`log2`/`log10`. -- **Vector ops** — `dot`/mean helpers and a correctly-rounded (exact-accumulating) `Sum` for - `CBig`. (Fold-based `Sum`/`Product` for `CBig` already exist; `Sum` is not yet correctly-rounded.) +- **Vector ops** — `dot`/mean helpers. (`Sum` for `CBig` is already exact-accumulating, matching + `FBig: Sum`; `Product` for `CBig` remains a fold, matching `FBig: Product`.) - **Independent re/im rounding** — a `CRound` trait giving MPC `mpc_rnd_t` parity (0.5 uses a single `R` for both parts). - **Third-party integration** — `CBig` `serde`/`rkyv`/`zeroize`, and diff --git a/complex/CHANGELOG.md b/complex/CHANGELOG.md index 4afe6d11..60769c08 100644 --- a/complex/CHANGELOG.md +++ b/complex/CHANGELOG.md @@ -7,6 +7,49 @@ parts for each of the real/imaginary components (built on `Repr::new_const`). The `cbig!` literal macro now works in `const` position for coefficients that fit in a `DoubleWord`; larger coefficients fall back to the runtime heap path. +- **Correct rounding for the complex transcendentals** via a Ziv retry loop (`complex/src/ziv.rs`). + `exp`, `log`, `powf`, `sin`/`cos`/`tan`/`sin_cos`, `asin`/`acos`/`atan`, and `sqrt` now certify + *both* the real and imaginary parts — rounding each to the target precision and retrying with more + guard digits while either part's error interval straddles a rounding boundary — matching + `dashu-float`'s real transcendentals. Each transcendental reports a provable per-part error radius + (`result.ulp() × C`, plus an amplification term where the composition magnifies error: `log`'s + `ln|z|` near `|z| = 1`, `powf`'s result magnitude). `tan` uses the cancellation-free double-angle + form `(sin 2x + i·sinh 2y)/(cos 2x + cosh 2y)` (the naive `sin z/cos z` cancels in the real part for + large `|Im z|`), so it is accurate for all finite `|Im z|`. `abs` delegates directly to + `dashu-float`'s already-correctly-rounded `hypot`, dropping a double-rounding re-round. The + well-conditioned regime is guaranteed-correctly rounded; `asin`/`acos` compute the inner `1-z²` in + the factored form `(1-z)(1+z)` (Sterbenz-exact near the singularities `z = ±1`, where the direct + `1 - z²` would cancel against the `sqr` rounding error), so they stay accurate right up to `z = ±1`. + (`atan` was already sound near its singularities `z = ±i`: its `1 ± iz` is built from an exact + rotation and an exact-significand add/sub, so there is no cancellation to factor out.) + +### Change +- **(breaking, bound)** The complex transcendentals (`exp`, `ln`, `powf`, `powi`, `sin`/`cos`/`tan`/ + `sin_cos`, `asin`/`acos`/`atan`) and their `CachedCBig` forwarders now require `R: ErrorBounds`, + inherited from `dashu-float`'s Ziv-backed real `exp`/`ln`/`sinh_cosh` primitives. `powi` is now + correctly rounded via its own Ziv loop (repeated squaring, with a per-part radius that scales with + the squaring-compounding error and drops to zero when the chain is exact — required under the + directed rounding modes, `CBig`'s default). The field arithmetic + (`add`/`sub`/`mul`/`div`/`sqr`/`inv`) remains `R: Round`. + +### Fix +- **`Sum` for `CBig` is now correctly rounded.** It exact-accumulates the real and imaginary + parts independently via `FBig`'s exact-accumulating `Sum` (lossless `Repr` accumulation, a single + final rounding per axis at the `Context::max` target), instead of the previous componentwise fold + whose per-step rounding could lose low-order terms (e.g. `1e20 + 1 - 1e20` folded to `0` at low + precision; it now correctly yields `1`). `Product` stays a fold, matching `FBig: Product`. +- **Unlimited-precision handling**, centralized in the Ziv driver: it asserts a limited context up + front, so every transcendental panics on precision 0 as documented (the exact special-value + shortcuts — `exp(0)=1`, `log(0)=-∞`, `powf(z,0)=1`, `sqrt(±0)`, `sin/cos(0)`, etc. — still bypass + it). The prior per-function `guard()`/`assert_limited()` scaffolding is removed (`guard()` itself + is gone, now unused). +- **Arithmetic at unlimited precision is correct.** `mul`/`sqr`/`norm` use `Context::work_context`, + which is the exact `self.float()` (precision 0) at unlimited, so they are exact there. `div`/`inv` + use the same path and panic at unlimited via `dashu-float`'s `div` (a quotient isn't exactly + representable in general). `abs` asserts limited before delegating to the float `hypot`. +- **`no_std` build of the test-only Ziv retry counter.** The `LAST_ZIV_RETRIES` `thread_local!` + requires `std`, so the crate failed to compile under `--no-default-features`. It's now gated behind + `feature = "std"` (with its uses and the counter-reading tests), so the Ziv driver is `no_std`-clean. ## 0.5.0 (Initial release) diff --git a/complex/src/cbig_cached_ops.rs b/complex/src/cbig_cached_ops.rs index 1d913bde..397749e4 100644 --- a/complex/src/cbig_cached_ops.rs +++ b/complex/src/cbig_cached_ops.rs @@ -11,7 +11,7 @@ use core::iter::{Product, Sum}; use core::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign}; use dashu_base::Inverse; -use dashu_float::round::Round; +use dashu_float::round::{ErrorBounds, Round}; use dashu_float::FBig; use dashu_int::{IBig, Word}; @@ -372,29 +372,35 @@ macro_rules! delegate_to_cbig { } impl CachedCBig { - // cache-threading transcendentals - forward_cached!(ln => log); - forward_cached!(exp => exp); - forward_cached!(sin => sin); - forward_cached!(cos => cos); - forward_cached!(tan => tan); - forward_cached!(asin => asin); - forward_cached!(acos => acos); - forward_cached!(atan => atan); - // non-cache delegations delegate_to_cbig!(sqr); - delegate_to_cbig!(sqrt); delegate_to_cbig!(conj); delegate_to_cbig!(proj); delegate_to_cbig!(mul_i(negative: bool)); - delegate_to_cbig!(powi(exp: IBig)); /// The squared modulus `re² + im²` (a real [`FBig`]) (see [`CBig::norm`]). #[inline] pub fn norm(&self) -> FBig { self.cbig.norm() } +} + +// Transcendentals route through the Ziv-backed (or Ziv-dependent) real/complex methods, which +// require `R: ErrorBounds`. +impl CachedCBig { + // cache-threading transcendentals + delegate_to_cbig!(powi(exp: IBig)); + forward_cached!(ln => log); + forward_cached!(exp => exp); + forward_cached!(sin => sin); + forward_cached!(cos => cos); + forward_cached!(tan => tan); + forward_cached!(asin => asin); + forward_cached!(acos => acos); + forward_cached!(atan => atan); + + // `sqrt` and `abs` route through `Context::hypot`, which is now Ziv-backed (`R: ErrorBounds`). + delegate_to_cbig!(sqrt); /// The modulus `|z|` (a real [`FBig`]) (see [`CBig::abs`]). #[inline] diff --git a/complex/src/div.rs b/complex/src/div.rs index 933eccdd..53bf7b77 100644 --- a/complex/src/div.rs +++ b/complex/src/div.rs @@ -21,7 +21,7 @@ impl Context { if z.is_zero() { return Ok(riemann(*self)); // 1/0 = ∞ } - let gctx = self.guard(DIV_GUARD); + let gctx = self.work_context(DIV_GUARD); let p = self.precision(); let (x, y) = (z.re(), z.im()); // n = x² + y² @@ -41,7 +41,7 @@ impl Context { if let Some(special) = div_special(z, w) { return special; } - let gctx = self.guard(DIV_GUARD); + let gctx = self.work_context(DIV_GUARD); let p = self.precision(); let (x, y) = (z.re(), z.im()); let (u, v) = (w.re(), w.im()); @@ -109,7 +109,7 @@ impl Context { } return Ok(riemann(*self)); // z/±0 (z≠0) = ∞ } - let gctx = self.guard(DIV_GUARD); + let gctx = self.work_context(DIV_GUARD); let p = self.precision(); let re = gctx.div(z.re(), s.repr())?.value().with_precision(p); let im = gctx.div(z.im(), s.repr())?.value().with_precision(p); @@ -268,4 +268,14 @@ mod tests { assert!(q_neg.is_infinite()); assert!(q_neg == q_pos); // both are +∞ + i·0 } + + // `div` at unlimited precision panics via the float layer (`work_context` → `self.float()` → + // float `div`, which rejects unlimited: a quotient isn't exactly representable in general). + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_div_unlimited_panics() { + let one = C::ONE; + let i = C::I; + let _ = &one / &i; // 1/i at unlimited + } } diff --git a/complex/src/exp.rs b/complex/src/exp.rs index b7aacb02..45ceb97c 100644 --- a/complex/src/exp.rs +++ b/complex/src/exp.rs @@ -10,9 +10,9 @@ use crate::cbig::CBig; use crate::repr::{combine_parts, exact, reborrow_cache, riemann, CfpResult, Context}; use dashu_base::Approximation::*; -use dashu_base::{BitTest, Sign}; -use dashu_float::round::Round; -use dashu_float::{ConstCache, FBig, FpError}; +use dashu_base::{Abs, BitTest, Sign}; +use dashu_float::round::ErrorBounds; +use dashu_float::{ConstCache, Context as FloatCtxt, FBig, FpError}; use dashu_int::{IBig, Word}; /// Guard digits (base-B) for `exp`. Composes a real `exp`, a `sin_cos`, and two products. @@ -22,9 +22,78 @@ const EXP_GUARD: usize = 14; /// cancellation-prone path, so a larger guard than the bare arithmetic ops. const POWF_GUARD: usize = 22; -impl Context { - /// Complex exponential under this context (context layer). Reuses `dashu-float`'s `exp` and - /// `sin_cos`; the cache is threaded into both (the convenience layer passes `None`). +impl Context { + /// Raise a complex number to an integer power under this context (context layer), correctly + /// rounded via a Ziv loop over the binary-exponentiation (repeated-squaring) chain. No cache. + /// + /// `powi(z, 0) = 1`; a negative exponent computes `(1/z)^|n|` directly, so the + /// sign-dependent overflow/underflow propagates from the closure with `?`. Repeated squaring + /// compounds the relative error (it roughly doubles per step), so after `bit_len(n)` squarings + /// the per-part error is bounded by about `2^nlen · ulp`, which the radius reflects; complex + /// `sqr`/`mul` are near-correct (a few ulp, not 0.5), so the bound carries an extra margin. + pub fn powi(&self, z: &CBig, exp: IBig) -> CfpResult { + let (sign, n) = exp.into_parts(); + if n.is_zero() { + return Ok(Exact(CBig::ONE)); + } + let negative = sign == Sign::Negative; + if n.is_one() { + // |n| == 1: z (positive) or 1/z (negative), a single op. + return if negative { + self.inv(z) + } else { + Ok(Exact(z.clone())) + }; + } + + let p = self.precision(); + let nlen = n.bit_len(); + // Initial guard scales with `nlen` (the squaring-compounding loss) plus a margin for the + // near-correct complex `sqr`/`mul`; sized so the first attempt certifies a non-tie result. + let initial_guard = nlen + 6; + let [re, im] = self.ziv(initial_guard, |guard| { + let pw = p + guard; + let gctx = Context::new(pw); + // start from z (positive exponent, always exact) or its working-precision reciprocal + // (negative exponent, exact only when 1/z is exactly representable). + let (start, mut exact) = if negative { + let inv = gctx.inv(z)?; + let ex = matches!(inv, Exact(_)); + (inv.value(), ex) + } else { + (z.clone(), true) + }; + // left-to-right binary exponentiation, tracking whether every step rounded Exact. + let mut acc = start.clone(); + for i in (0..nlen - 1).rev() { + let s = gctx.sqr(&acc)?; + exact = exact && matches!(s, Exact(_)); + acc = s.value(); + if n.bit(i) { + let m = gctx.mul(&acc, &start)?; + exact = exact && matches!(m, Exact(_)); + acc = m.value(); + } + } + let (re, im) = acc.into_parts(); + // re-root to the working precision (parts may be exact constants). + let re = re.with_precision(pw).value(); + let im = im.with_precision(pw).value(); + let shift = (nlen + 3) as isize; + // When the whole chain is exact the result is the mathematically exact zⁿ: report a zero + // radius. This is required under the directed rounding modes (the `CBig` default), where + // an exactly-representable result lies on a one-sided rounding boundary that no nonzero + // radius can fit inside. (See `dashu-float`'s `powi` for the same reasoning.) + let re_r = if exact { FBig::ZERO } else { re.ulp() << shift }; + let im_r = if exact { FBig::ZERO } else { im.ulp() << shift }; + Ok([(re.clone(), re_r), (im.clone(), im_r)]) + })?; + Ok(combine_parts(re, im)) + } + + /// Complex exponential under this context (context layer). Computes `e^x·(cos y + i·sin y)` + /// from `dashu-float`'s (correctly-rounded) `exp` and `sin_cos`, wrapped in a Ziv loop that + /// certifies both parts; the cache is threaded into both (the convenience layer passes `None`). /// /// Special values: `exp(0) = 1`; `exp(+inf + i·finite) = +∞` (Riemann point); /// `exp(-inf + i·finite) = 0`; an infinite imaginary part makes the trig undefined @@ -48,52 +117,32 @@ impl Context { }; } - let gctx = self.guard(EXP_GUARD); + // `e^x·(cos y + i·sin y)`. The float `exp`/`sin_cos` are correctly-rounded at the working + // precision, so each contributes ~0 to the composition radius; only the two products round, + // at a few working-ULPs. The Ziv driver asserts a limited context (the special-value + // shortcuts above are exact and need no precision); overflow from a large real part + // propagates from the closure with `?`. let p = self.precision(); - let ex = gctx.exp(z.re(), reborrow_cache(&mut cache))?.value(); - let (sin_y, cos_y) = gctx.sin_cos(z.im(), reborrow_cache(&mut cache)); - let cos_y = cos_y?.value(); - let sin_y = sin_y?.value(); - let re = gctx.mul(ex.repr(), cos_y.repr())?.value().with_precision(p); - let im = gctx.mul(ex.repr(), sin_y.repr())?.value().with_precision(p); + let [re, im] = self.ziv(EXP_GUARD, |guard| { + let gctx = FloatCtxt::::new(p + guard); + let ex = gctx.exp(z.re(), reborrow_cache(&mut cache))?.value(); + let (sin_y, cos_y) = gctx.sin_cos(z.im(), reborrow_cache(&mut cache)); + let cos_y = cos_y?.value(); + let sin_y = sin_y?.value(); + let re = gctx.mul(ex.repr(), cos_y.repr())?.value(); + let im = gctx.mul(ex.repr(), sin_y.repr())?.value(); + Ok([(re.clone(), re.ulp() * 6), (im.clone(), im.ulp() * 6)]) + })?; Ok(combine_parts(re, im)) } - /// Raise a complex number to an integer power under this context (context layer), via repeated - /// squaring (branch-cut-free, cheaper than `exp(n·log z)`). No cache. - /// - /// `powi(z, 0) = 1`; a negative exponent computes `powi(z, |n|)` then inverts. - pub fn powi(&self, z: &CBig, exp: IBig) -> CfpResult { - let (sign, n) = exp.into_parts(); - if n.is_zero() { - return Ok(Exact(CBig::ONE)); - } - let negative = sign == Sign::Negative; - let bitlen = n.bit_len(); - // left-to-right binary exponentiation, starting from the leading set bit - let mut acc = z.clone(); - for i in (0..bitlen - 1).rev() { - acc = self.sqr(&acc)?.value(); - if n.bit(i) { - acc = self.mul(&acc, z)?.value(); - } - } - // The intermediate rounding flags are folded away (the value is near-correctly rounded); - // for a negative exponent the final `inv` carries its own flags. - if negative { - self.inv(&acc) - } else { - Ok(Exact(acc)) - } - } - /// Raise `base` to a complex power under this context (context layer): `exp(w·log base)` on the - /// principal branch, evaluated at `p + POWF_GUARD` and re-rounded. `powf(0, 0) = 1` (matching - /// `FBig::powf`). + /// principal branch, correctly rounded via a Ziv loop. `powf(0, 0) = 1` (matching `FBig::powf`). /// - /// Unlike `exp`, this drives whole-[`CBig`] operations (`log`/`mul`/`exp`), so it builds a - /// complex working [`Context`] at guard precision directly rather than the float - /// `Context::guard` (which yields a `FloatCtxt` for per-part math). + /// The result's error is amplified by the exponent magnitude `‖w·log base‖`: the outer `exp` + /// multiplies the error in `w·log base` by the result magnitude, so the per-part radius carries + /// a data-dependent `‖w·log base‖` factor (mirroring `FBig::powf`). Overflow (a large exponent) + /// propagates from the closure with `?`. pub fn powf( &self, base: &CBig, @@ -103,35 +152,49 @@ impl Context { if w.is_zero() { return Ok(Exact(CBig::ONE)); // powf(z, 0) = 1, incl. powf(0, 0) } - let gctx = Context::new(self.precision() + POWF_GUARD); - let log_z = gctx.log(base, reborrow_cache(&mut cache))?.value(); - let wlogz = gctx.mul(w, &log_z)?.value(); - let hi = gctx.exp(&wlogz, reborrow_cache(&mut cache))?.value(); let p = self.precision(); - let (re, im) = hi.into_parts(); - Ok(combine_parts(re.with_precision(p), im.with_precision(p))) + let [re, im] = self.ziv(POWF_GUARD, |guard| { + let pw = p + guard; + let gctx = Context::new(pw); + let log_z = gctx.log(base, reborrow_cache(&mut cache))?.value(); + let wlogz = gctx.mul(w, &log_z)?.value(); + let hi = gctx.exp(&wlogz, reborrow_cache(&mut cache))?.value(); + // The outer `exp` amplifies the error in `w·log base` by the result magnitude, so the + // radius scales with `‖w·log base‖`. `re()`/`im()` are raw `Repr`s — wrap to take the + // absolute value; the L1 norm `|re|+|im|` upper-bounds the magnitude. + let fctx = gctx.float(); + let l1 = FBig::::from_repr(wlogz.re().clone(), fctx).abs() + + FBig::::from_repr(wlogz.im().clone(), fctx).abs(); + let amp = (l1 + FBig::::ONE) * 16i32; + let (re, im) = hi.into_parts(); + // re-root to the working precision (`exp`/`log` may return exact constants). + let re = re.with_precision(pw).value(); + let im = im.with_precision(pw).value(); + Ok([(re.clone(), re.ulp() * &), (im.clone(), im.ulp() * &)]) + })?; + Ok(combine_parts(re, im)) } } -impl CBig { - /// Complex exponential `e^z` (convenience layer). +impl CBig { + /// Integer power (convenience layer). /// /// # Panics /// - /// Panics if the precision is unlimited or on an indeterminate special value. + /// Panics on an indeterminate / out-of-domain result (e.g. `0⁻¹`). #[inline] - pub fn exp(&self) -> Self { - self.context().unwrap_cfp(self.context().exp(self, None)) + pub fn powi(&self, exp: IBig) -> Self { + self.context().unwrap_cfp(self.context().powi(self, exp)) } - /// Integer power (convenience layer). + /// Complex exponential `e^z` (convenience layer). /// /// # Panics /// - /// Panics on an indeterminate / out-of-domain result (e.g. `0⁻¹`). + /// Panics if the precision is unlimited or on an indeterminate special value. #[inline] - pub fn powi(&self, exp: IBig) -> Self { - self.context().unwrap_cfp(self.context().powi(self, exp)) + pub fn exp(&self) -> Self { + self.context().unwrap_cfp(self.context().exp(self, None)) } /// Complex power `self^w` (convenience layer). @@ -165,8 +228,10 @@ mod tests { #[test] fn exp_one_is_e() { - // exp(1+0i) = e ≈ 2.71828…; check 2 < e < 3 via the real part - let e = C::ONE.exp(); + // exp(1+0i) = e ≈ 2.71828…; check 2 < e < 3 via the real part. Use a *limited*-precision + // input — `exp` rejects unlimited precision (it would otherwise silently compute at the + // fixed `EXP_GUARD`). + let e = c(1, 0).exp(); let (re, _im) = e.into_parts(); assert!(re > F::from(2)); assert!(re < F::from(3)); @@ -195,6 +260,19 @@ mod tests { assert!(r.im().is_pos_zero()); } + #[test] + fn exp_huge_real_overflows() { + // exp of a huge real part overflows the isize exponent range. The error propagates from the + // Ziv closure via `?` (no hoisted probe), and the convenience layer saturates it to +∞. + let huge = F::from_parts(IBig::from(1) << 100, 0) + .with_precision(53) + .value(); + let z = CBig::from_parts(huge, F::from(0).with_precision(53).value()); + let e = z.exp(); + assert!(e.re().is_infinite()); + assert_eq!(e.re().sign(), Sign::Positive); + } + #[test] fn powi_zero_is_one() { assert!(c(3, 4).powi(0.into()) == C::ONE); @@ -233,4 +311,18 @@ mod tests { let z = c(2, 1); assert!(z.powf(&C::ONE) == z); } + + // exp/powf on an unlimited-precision CBig must panic, not silently compute at the fixed guard + // precision (`C::ONE` / `C::ZERO` are unlimited-precision constants). + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_exp_unlimited_precision_panics() { + let _ = C::ONE.exp(); + } + + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_powf_unlimited_precision_panics() { + let _ = C::ONE.powf(&C::ONE); + } } diff --git a/complex/src/iter.rs b/complex/src/iter.rs index acd26150..b0376ab4 100644 --- a/complex/src/iter.rs +++ b/complex/src/iter.rs @@ -1,24 +1,38 @@ //! Implementation of core::iter traits for [`CBig`]. //! -//! `Sum`/`Product` fold with the binary `+`/`*` operators (componentwise perfectly-rounded real -//! ops on each part). The impls are concrete (`Sum`/`Sum<&CBig>`, `Product`/`Product<&CBig>`) -//! rather than generic over `Add`/`Mul`, matching the narrowed iter surface used for `FBig`. A -//! correctly-rounded (exact-accumulating) `Sum` for `CBig` is a possible future refinement. +//! `Sum` is exact-accumulating: since complex addition is componentwise, the real and imaginary +//! parts are summed independently via `FBig`'s exact-accumulating `Sum` (each axis accumulates its +//! `Repr`s losslessly and rounds once), so the result is the correctly-rounded complex sum. `Product` +//! folds with the binary `*` operator (componentwise near-correct). The impls are concrete +//! (`Sum`/`Sum<&CBig>`, `Product`/`Product<&CBig>`) rather than generic over `Add`/`Mul`, matching +//! the narrowed iter surface used for `FBig`. use crate::cbig::CBig; +use alloc::vec::Vec; use core::iter::{Product, Sum}; use dashu_float::round::Round; +use dashu_float::FBig; use dashu_int::Word; +/// Exact-accumulating complex sum: split into real/imaginary [`FBig`] streams and sum each via +/// `FBig`'s `Sum` (which accumulates `Repr`s exactly and rounds once at the target context). An +/// empty stream sums to `CBig::ZERO` (both axes are `FBig::ZERO`). +fn precise_sum( + parts: impl Iterator, FBig)>, +) -> CBig { + let (re_parts, im_parts): (Vec>, Vec>) = parts.unzip(); + CBig::from_parts(re_parts.into_iter().sum(), im_parts.into_iter().sum()) +} + impl Sum for CBig { fn sum>(iter: I) -> Self { - iter.fold(CBig::ZERO, |acc, x| acc + x) + precise_sum(iter.map(CBig::into_parts)) } } impl<'a, R: Round, const B: Word> Sum<&'a CBig> for CBig { fn sum>(iter: I) -> Self { - iter.fold(CBig::ZERO, |acc, x| acc + x) + precise_sum(iter.map(|c| c.clone().into_parts())) } } @@ -39,6 +53,7 @@ mod tests { use crate::cbig::CBig; use dashu_float::round::mode::HalfAway; use dashu_float::FBig; + use dashu_int::IBig; type C = CBig; type F = FBig; @@ -77,4 +92,31 @@ mod tests { let p: C = core::iter::empty::().product(); assert_eq!(p, C::ONE); } + + #[test] + fn sum_exact_cancellation() { + // z + (-z) = +0 on both axes (canonical +0, not -0, under HalfAway). + let z = C::from_parts(F::from(3), F::from(4)); + let s: C = [z, -C::from_parts(F::from(3), F::from(4))] + .into_iter() + .sum(); + assert_eq!(s, C::ZERO); + } + + #[test] + fn sum_exact_accumulates_wide_magnitudes() { + // At precision 4, (1e20 + 1 - 1e20) folds to 0 (the +1 is lost when added to 1e20), but + // exact accumulation preserves it: the sum is 1 on each axis. This distinguishes the + // exact-accumulating `Sum` from a naive fold. + let mk = |sig: i32, exp: isize| { + F::from_parts(IBig::from(sig), exp) + .with_precision(4) + .value() + }; + let big = C::from_parts(mk(1, 20), mk(1, 20)); + let one = C::from_parts(mk(1, 0), mk(1, 0)); + let neg_big = C::from_parts(mk(-1, 20), mk(-1, 20)); + let s: C = [big, one, neg_big].into_iter().sum(); + assert_eq!(s, C::from_parts(mk(1, 0), mk(1, 0))); + } } diff --git a/complex/src/lib.rs b/complex/src/lib.rs index f7c47a60..829d0415 100644 --- a/complex/src/lib.rs +++ b/complex/src/lib.rs @@ -85,6 +85,7 @@ mod parse; mod repr; mod root; mod third_party; +mod ziv; // All the public items from third_party will be exposed #[allow(unused_imports)] diff --git a/complex/src/log.rs b/complex/src/log.rs index d79fb9d5..dfec708b 100644 --- a/complex/src/log.rs +++ b/complex/src/log.rs @@ -2,14 +2,14 @@ use crate::cbig::CBig; use crate::repr::{combine_parts, exact, reborrow_cache, riemann, CfpResult, Context}; -use dashu_float::round::Round; -use dashu_float::{ConstCache, FBig, Repr}; -use dashu_int::Word; +use dashu_float::round::ErrorBounds; +use dashu_float::{ConstCache, Context as FloatCtxt, FBig, Repr}; +use dashu_int::{IBig, Word}; /// Guard digits (base-B) for `log`. Composes `hypot` (for `|z|`), `ln`, and `atan2`. const LOG_GUARD: usize = 14; -impl Context { +impl Context { /// Complex natural logarithm under this context (context layer). `log z = ln|z| + i·arg(z)`, /// with the imaginary part in `]−π, π]`. The cache threads into `ln` and `atan2`. /// @@ -31,22 +31,40 @@ impl Context { return Ok(riemann(*self)); // log(∞) = +∞ (Riemann point) } - let gctx = self.guard(LOG_GUARD); + // `ln|z| + i·arg(z)`. The float `hypot`/`ln`/`atan2` are correctly-rounded at the working + // precision. The imaginary part (`atan2` of the exact parts) carries only `atan2`'s own + // rounding; the real part `ln|z|` additionally propagates `hypot`'s relative error through + // `ln`, which dominates near `|z| = 1` (where `ln|z| → 0`) — so its radius carries an extra + // absolute `B^{1-pw}` term. The Ziv driver asserts a limited context (the `log(0)`/`log(∞)` + // shortcuts above are exact and need no precision). let p = self.precision(); - // ln|z| - let r = gctx.hypot(z.re(), z.im())?.value(); - let ln_r = gctx.ln(r.repr(), reborrow_cache(&mut cache))?.value(); - // arg(z) = atan2(im, re) - let arg = gctx - .atan2(z.im(), z.re(), reborrow_cache(&mut cache))? - .value(); - let re = ln_r.with_precision(p); - let im = arg.with_precision(p); + let [re, im] = self.ziv(LOG_GUARD, |guard| { + let pw = p + guard; + let gctx = FloatCtxt::::new(pw); + // ln|z| + let r = gctx.hypot(z.re(), z.im())?.value(); + let ln_r = gctx.ln(r.repr(), reborrow_cache(&mut cache))?.value(); + // arg(z) = atan2(im, re) + let arg = gctx + .atan2(z.im(), z.re(), reborrow_cache(&mut cache))? + .value(); + // The float transcendentals return unlimited-precision exact constants for exact cases + // (e.g. `ln 1 = 0`, `atan2(0,1) = 0`); re-root to the working precision so `.ulp()` + // (which rejects unlimited) is well-defined. + let ln_r = ln_r.with_precision(pw).value(); + let arg = arg.with_precision(pw).value(); + // `B^{1-pw}` upper-bounds `hypot`'s propagated error `ulp(r)/|r| ≤ B^{1-pw}`, which + // dominates `ulp(ln_r)` when `|ln_r| < 1` (`|z| ≈ 1`). + let propagated = FBig::::from_parts(IBig::from(1), 1 - pw as isize); + let re_rad = ln_r.ulp() * 4 + propagated * 4; + let im_rad = arg.ulp() * 4; + Ok([(ln_r, re_rad), (arg, im_rad)]) + })?; Ok(combine_parts(re, im)) } } -impl CBig { +impl CBig { /// Complex natural logarithm (principal branch; convenience layer). /// /// # Panics @@ -82,7 +100,9 @@ mod tests { #[test] fn ln_one_is_zero() { - assert!(C::ONE.ln() == C::ZERO); + // Use a *limited*-precision input — `log` rejects unlimited precision (it would otherwise + // silently compute at the fixed `LOG_GUARD`). + assert!(c(1, 0).ln() == C::ZERO); } #[test] @@ -102,4 +122,12 @@ mod tests { assert!(l.re().is_infinite()); assert_eq!(l.re().sign(), Sign::Negative); } + + // log on an unlimited-precision CBig must panic, not silently compute at LOG_GUARD digits + // (the `log(0)` / `log(∞)` shortcuts above are exact, so they don't hit this). + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_log_unlimited_precision_panics() { + let _ = C::ONE.ln(); + } } diff --git a/complex/src/math/trig.rs b/complex/src/math/trig.rs index 9f3fabfb..b7d7b983 100644 --- a/complex/src/math/trig.rs +++ b/complex/src/math/trig.rs @@ -6,17 +6,18 @@ use crate::cbig::CBig; use crate::repr::{combine_parts, reborrow_cache, CfpResult, Context}; -use dashu_float::round::Round; -use dashu_float::{ConstCache, FBig, FpError, Repr}; +use dashu_float::round::ErrorBounds; +use dashu_float::{ConstCache, Context as FloatCtxt, FBig, FpError, Repr}; use dashu_int::{IBig, Word}; /// Guard digits (base-B) for the forward trig. Composes real `sin_cos` + `sinh_cosh` + two /// products; the cancellation near the trig zeros is absorbed by the re-round. const TRIG_GUARD: usize = 16; -impl Context { - /// Simultaneously compute `sin z` and `cos z` (context layer). Returns `(sin, cos)` each as a - /// [`CfpResult`]. An infinite input maps to [`FpError::Indeterminate`] (the C99 NaN cases). +impl Context { + /// Simultaneously compute `sin z` and `cos z` (context layer), correctly rounded via a shared + /// Ziv loop. Returns `(sin, cos)` each as a [`CfpResult`]. An infinite input maps to + /// [`FpError::Indeterminate`] (the C99 NaN cases). pub fn sin_cos( &self, z: &CBig, @@ -37,48 +38,35 @@ impl Context { return (zero, one); } - let gctx = self.guard(TRIG_GUARD); + // `sin z = sinx·coshy + i·cosx·sinhy`, `cos z = cosx·coshy − i·sinx·sinhy`. The four products + // share one evaluation of the real `sin_cos`/`sinh_cosh` (each correctly-rounded at the + // working precision, contributing ~0); only the products round, a few working-ULPs each. A + // single 4-part Ziv loop certifies all of `sin` and `cos` together. let p = self.precision(); - let (sinx, cosx) = gctx.sin_cos(z.re(), reborrow_cache(&mut cache)); - let sinx = match sinx { - Ok(v) => v.value(), - Err(e) => return (Err(e), Err(FpError::Indeterminate)), - }; - let cosx = match cosx { - Ok(v) => v.value(), - Err(e) => return (Err(FpError::Indeterminate), Err(e)), - }; - let (sinhy_res, coshy_res) = gctx.sinh_cosh(z.im(), reborrow_cache(&mut cache)); - let sinhy = match sinhy_res { - Ok(v) => v.value(), - Err(e) => return (Err(e), Err(FpError::Indeterminate)), - }; - let coshy = match coshy_res { - Ok(v) => v.value(), - Err(e) => return (Err(FpError::Indeterminate), Err(e)), - }; - - // sin z = (sinx·coshy) + i·(cosx·sinhy); cos z = (cosx·coshy) − i·(sinx·sinhy). - // `sin_cos` returns a tuple, so the products are matched explicitly (no `?`). - let prod = |a: &FBig, b: &FBig| -> Result<_, FpError> { - Ok(gctx.mul(a.repr(), b.repr())?.value().with_precision(p)) - }; - let sin_re = match prod(&sinx, &coshy) { - Ok(v) => v, - Err(e) => return (Err(e), Err(FpError::Indeterminate)), - }; - let sin_im = match prod(&cosx, &sinhy) { - Ok(v) => v, - Err(e) => return (Err(e), Err(FpError::Indeterminate)), - }; - let cos_re = match prod(&cosx, &coshy) { - Ok(v) => v, - Err(e) => return (Err(FpError::Indeterminate), Err(e)), - }; - let neg_sinx = -sinx; - let cos_im = match prod(&neg_sinx, &sinhy) { - Ok(v) => v, - Err(e) => return (Err(FpError::Indeterminate), Err(e)), + let parts = self.ziv(TRIG_GUARD, |guard| { + let gctx = FloatCtxt::::new(p + guard); + let (sinx, cosx) = gctx.sin_cos(z.re(), reborrow_cache(&mut cache)); + let sinx = sinx?.value(); + let cosx = cosx?.value(); + let (sinhy, coshy) = gctx.sinh_cosh(z.im(), reborrow_cache(&mut cache)); + let sinhy = sinhy?.value(); + let coshy = coshy?.value(); + let sin_re = gctx.mul(sinx.repr(), coshy.repr())?.value(); + let sin_im = gctx.mul(cosx.repr(), sinhy.repr())?.value(); + let cos_re = gctx.mul(cosx.repr(), coshy.repr())?.value(); + let neg_sinx = -sinx; // cos z's imaginary part is −sinx·sinhy + let cos_im = gctx.mul(neg_sinx.repr(), sinhy.repr())?.value(); + Ok([ + (sin_re.clone(), sin_re.ulp() * 8), + (sin_im.clone(), sin_im.ulp() * 8), + (cos_re.clone(), cos_re.ulp() * 8), + (cos_im.clone(), cos_im.ulp() * 8), + ]) + }); + let [sin_re, sin_im, cos_re, cos_im] = match parts { + Ok(arr) => arr, + // an overflow (e.g. `cosh` of a huge imaginary part) fails both sin and cos together. + Err(e) => return (Err(e), Err(e)), }; (Ok(combine_parts(sin_re, sin_im)), Ok(combine_parts(cos_re, cos_im))) } @@ -103,21 +91,56 @@ impl Context { self.sin_cos(z, cache).1 } - /// Complex tangent `sin z / cos z` (context layer). + /// Complex tangent (context layer), correctly rounded via a Ziv loop, using the cancellation-free + /// double-angle identity + /// + /// `tan(x+iy) = (sin 2x + i·sinh 2y) / (cos 2x + cosh 2y)`. + /// + /// The denominator `cos 2x + cosh 2y` is a sum of a bounded term (`cos 2x ∈ [−1, 1]`) and a + /// term `≥ 1` (`cosh 2y`), so it never catastrophically cancels — unlike `sin z / cos z`, whose + /// `sin·conj(cos)` real part cancels from `~cosh²y` down to `O(1)` for large `|Im z|`. The result + /// is accurate for all finite `|Im z|`; the only small-denominator points are the real-axis poles + /// (`y = 0, x = π/2 + kπ`), where the large value is genuine, not an artifact. pub fn tan( &self, z: &CBig, - cache: Option<&mut ConstCache>, + mut cache: Option<&mut ConstCache>, ) -> CfpResult { - let (sin_z, cos_z) = self.sin_cos(z, cache); - let sin_z = sin_z?; - let cos_z = cos_z?; - self.div(&sin_z.value(), &cos_z.value()) + let p = self.precision(); + let [re, im] = self.ziv(TRIG_GUARD, |guard| { + let pw = p + guard; + let gctx = FloatCtxt::::new(pw); + // 2x, 2y (exact doublings — same significand, exponent +1). + let x2 = gctx.add(z.re(), z.re())?.value(); + let y2 = gctx.add(z.im(), z.im())?.value(); + let (sin2x, cos2x) = gctx.sin_cos(x2.repr(), reborrow_cache(&mut cache)); + let sin2x = sin2x?.value(); + let cos2x = cos2x?.value(); + let (sinh2y, cosh2y) = gctx.sinh_cosh(y2.repr(), reborrow_cache(&mut cache)); + let sinh2y = sinh2y?.value(); + let cosh2y = cosh2y?.value(); + // D = cos 2x + cosh 2y (a benign sum: a bounded term plus one ≥ 1). + let denom = gctx.add(cos2x.repr(), cosh2y.repr())?.value(); + let re = gctx.div(sin2x.repr(), denom.repr())?.value(); + let im = gctx.div(sinh2y.repr(), denom.repr())?.value(); + // re-root to the working precision (`sin_cos`/`sinh_cosh`/`div` may return exact + // constants for exact cases such as `tan(0) = 0`). + let re = re.with_precision(pw).value(); + let im = im.with_precision(pw).value(); + Ok([(re.clone(), re.ulp() * 8), (im.clone(), im.ulp() * 8)]) + })?; + Ok(combine_parts(re, im)) } - /// Inverse sine `asin z = -i·log(iz + sqrt(1-z²))` (context layer, Kahan form). The argument of - /// the inner `log` always has positive real part, so the branch cut comes entirely from the - /// `sqrt`; an infinite input maps to [`FpError::Indeterminate`]. + /// Inverse sine `asin z = -i·log(iz + sqrt(1-z²))` (context layer, Kahan form), correctly + /// rounded via a Ziv loop. The argument of the inner `log` always has positive real part, so the + /// branch cut comes entirely from the `sqrt`; an infinite input maps to + /// [`FpError::Indeterminate`]. The `1-z²` under the `sqrt` is computed in the factored form + /// `(1-z)(1+z)`, which is Sterbenz-exact near `z = ±1` (where the direct `1-z²` would + /// catastrophically cancel against the `sqr` rounding error), so the well-conditioned regime + /// extends right up to the singularities. A generous constant radius covers the + /// square/subtract/sqrt/add/log composition; the Ziv retries absorb the `sqrt` amplification as + /// `1-z² → 0`. pub fn asin( &self, z: &CBig, @@ -126,21 +149,35 @@ impl Context { if z.is_infinite() { return Err(FpError::Indeterminate); } - let gctx = Context::new(self.precision() + ITRIG_GUARD); let p = self.precision(); - let one = CBig::ONE; - let z2 = gctx.sqr(z)?.value(); - let one_m_z2 = gctx.sub(&one, &z2)?.value(); - let sqrt_term = gctx.sqrt(&one_m_z2)?.value(); - let iz = z.mul_i(false); // exact rotation - let w = gctx.add(&iz, &sqrt_term)?.value(); - let log_w = gctx.log(&w, reborrow_cache(&mut cache))?.value(); - let asin_z = log_w.mul_i(true); // -i·log(w) - let (re, im) = asin_z.into_parts(); - Ok(combine_parts(re.with_precision(p), im.with_precision(p))) + let [re, im] = self.ziv(ITRIG_GUARD, |guard| { + let pw = p + guard; + let gctx = Context::new(pw); + let one = CBig::ONE; + // Factor `1-z² = (1-z)(1+z)`. Near `z = ±1` the direct `1 - z²` subtracts a value + // dominated by the `sqr` rounding error from 1 (catastrophic cancellation); the factored + // form is Sterbenz-exact there (`1-z` is computed exactly, since the subtraction's + // significand difference is exact), so the radius stays sound right up to the singularity. + let one_m_z = gctx.sub(&one, z)?.value(); + let one_p_z = gctx.add(&one, z)?.value(); + let one_m_z2 = gctx.mul(&one_m_z, &one_p_z)?.value(); + let sqrt_term = gctx.sqrt(&one_m_z2)?.value(); + let iz = z.mul_i(false); // exact rotation + let w = gctx.add(&iz, &sqrt_term)?.value(); + let log_w = gctx.log(&w, reborrow_cache(&mut cache))?.value(); + let asin_z = log_w.mul_i(true); // -i·log(w) + let (re, im) = asin_z.into_parts(); + // re-root to the working precision (`log` may return an exact constant for exact cases). + let re = re.with_precision(pw).value(); + let im = im.with_precision(pw).value(); + Ok([(re.clone(), re.ulp() * 20), (im.clone(), im.ulp() * 20)]) + })?; + Ok(combine_parts(re, im)) } - /// Inverse cosine `acos z = -i·log(z + i·sqrt(1-z²))` (context layer, Kahan form). + /// Inverse cosine `acos z = -i·log(z + i·sqrt(1-z²))` (context layer, Kahan form), correctly + /// rounded via a Ziv loop. Same composition and singularity structure as `asin` (including the + /// factored `1-z² = (1-z)(1+z)` near `z = ±1`). pub fn acos( &self, z: &CBig, @@ -149,21 +186,32 @@ impl Context { if z.is_infinite() { return Err(FpError::Indeterminate); } - let gctx = Context::new(self.precision() + ITRIG_GUARD); let p = self.precision(); - let one = CBig::ONE; - let z2 = gctx.sqr(z)?.value(); - let one_m_z2 = gctx.sub(&one, &z2)?.value(); - let sqrt_term = gctx.sqrt(&one_m_z2)?.value(); - let i_sqrt = sqrt_term.mul_i(false); // i·sqrt(1-z²) - let w = gctx.add(z, &i_sqrt)?.value(); - let log_w = gctx.log(&w, reborrow_cache(&mut cache))?.value(); - let acos_z = log_w.mul_i(true); // -i·log(w) - let (re, im) = acos_z.into_parts(); - Ok(combine_parts(re.with_precision(p), im.with_precision(p))) + let [re, im] = self.ziv(ITRIG_GUARD, |guard| { + let pw = p + guard; + let gctx = Context::new(pw); + let one = CBig::ONE; + // Factored `1-z² = (1-z)(1+z)` — Sterbenz-exact near `z = ±1` (see `asin`). + let one_m_z = gctx.sub(&one, z)?.value(); + let one_p_z = gctx.add(&one, z)?.value(); + let one_m_z2 = gctx.mul(&one_m_z, &one_p_z)?.value(); + let sqrt_term = gctx.sqrt(&one_m_z2)?.value(); + let i_sqrt = sqrt_term.mul_i(false); // i·sqrt(1-z²) + let w = gctx.add(z, &i_sqrt)?.value(); + let log_w = gctx.log(&w, reborrow_cache(&mut cache))?.value(); + let acos_z = log_w.mul_i(true); // -i·log(w) + let (re, im) = acos_z.into_parts(); + let re = re.with_precision(pw).value(); + let im = im.with_precision(pw).value(); + Ok([(re.clone(), re.ulp() * 20), (im.clone(), im.ulp() * 20)]) + })?; + Ok(combine_parts(re, im)) } - /// Inverse tangent `atan z = (i/2)·(log(1-iz) - log(1+iz))` (context layer). + /// Inverse tangent `atan z = (i/2)·(log(1-iz) - log(1+iz))` (context layer), correctly rounded + /// via a Ziv loop. The two logs nearly cancel for small `z`; near `z = ±i` one of `1∓iz` + /// vanishes and its log diverges. The Ziv retries absorb the cancellation in the + /// well-conditioned regime. pub fn atan( &self, z: &CBig, @@ -174,27 +222,33 @@ impl Context { // 1±iz terms become infinite and the log diverges — report Indeterminate for now. return Err(FpError::Indeterminate); } - let gctx = Context::new(self.precision() + ITRIG_GUARD); let p = self.precision(); - let one = CBig::ONE; - let iz = z.mul_i(false); - let a = gctx.sub(&one, &iz)?.value(); // 1 - iz - let b = gctx.add(&one, &iz)?.value(); // 1 + iz - let log_a = gctx.log(&a, reborrow_cache(&mut cache))?.value(); - let log_b = gctx.log(&b, reborrow_cache(&mut cache))?.value(); - let diff = gctx.sub(&log_a, &log_b)?.value(); - let i_half_diff = diff.mul_i(false); // i·diff, then /2 below - let two: CBig = IBig::from(2).into(); - let atan_z = gctx.div(&i_half_diff, &two)?.value(); - let (re, im) = atan_z.into_parts(); - Ok(combine_parts(re.with_precision(p), im.with_precision(p))) + let [re, im] = self.ziv(ITRIG_GUARD, |guard| { + let pw = p + guard; + let gctx = Context::new(pw); + let one = CBig::ONE; + let iz = z.mul_i(false); + let a = gctx.sub(&one, &iz)?.value(); // 1 - iz + let b = gctx.add(&one, &iz)?.value(); // 1 + iz + let log_a = gctx.log(&a, reborrow_cache(&mut cache))?.value(); + let log_b = gctx.log(&b, reborrow_cache(&mut cache))?.value(); + let diff = gctx.sub(&log_a, &log_b)?.value(); + let i_half_diff = diff.mul_i(false); // i·diff, then /2 below + let two: CBig = IBig::from(2).into(); + let atan_z = gctx.div(&i_half_diff, &two)?.value(); + let (re, im) = atan_z.into_parts(); + let re = re.with_precision(pw).value(); + let im = im.with_precision(pw).value(); + Ok([(re.clone(), re.ulp() * 20), (im.clone(), im.ulp() * 20)]) + })?; + Ok(combine_parts(re, im)) } } /// Guard digits (base-B) for the inverse trig (squares, a sqrt, logs, and a divide). const ITRIG_GUARD: usize = 18; -impl CBig { +impl CBig { /// Complex sine (convenience layer). Panics on an indeterminate special value. #[inline] pub fn sin(&self) -> Self { @@ -279,24 +333,38 @@ mod tests { assert!(im.abs_cmp(&F::from_parts(1.into(), -12)).is_le()); } + #[test] + fn tan_large_imaginary_is_near_i() { + use dashu_base::{Abs, AbsOrd}; + // tan(x + i·100) ≈ i: real part → 0, imaginary → tanh(100) ≈ 1. The cancellation-free + // double-angle form computes this accurately; the naive `sin/cos` division would cancel the + // real part to noise for such a large `|Im z|` (the motivating case for the new formula). + let (re, im) = c(1, 100).tan().into_parts(); + let tol = F::from_parts(1.into(), -40); + assert!(re.abs().abs_cmp(&tol).is_le()); + assert!((im - F::from(1)).abs().abs_cmp(&tol).is_le()); + } + #[test] fn sin_i_is_i_sinh_one() { - // sin(i) = i·sinh(1) = i·1.1752… ; purely imaginary - let s = C::I.sin(); + // sin(i) = i·sinh(1) = i·1.1752… ; purely imaginary. Use a *limited*-precision input — + // `sin` rejects unlimited precision (it would otherwise silently compute at `TRIG_GUARD`). + let s = c(0, 1).sin(); assert!(s.re().significand().is_zero()); assert!(!s.im().significand().is_zero()); } #[test] fn asin_zero_is_zero() { - assert!(C::ZERO.asin() == C::ZERO); + // limited-precision input (asin rejects unlimited precision) + assert!(c(0, 0).asin() == C::ZERO); } #[test] fn asin_one_is_half_pi() { use dashu_base::{Abs, AbsOrd}; - // asin(1) = π/2 - let (re, im) = C::ONE.asin().into_parts(); + // asin(1) = π/2 (limited-precision input — asin rejects unlimited precision) + let (re, im) = c(1, 0).asin().into_parts(); let half_pi = F::from_parts(15707963267948966i64.into(), -16) .with_precision(60) .value(); @@ -310,7 +378,8 @@ mod tests { #[test] fn acos_zero_is_half_pi() { use dashu_base::{Abs, AbsOrd}; - let (re, _im) = C::ZERO.acos().into_parts(); + // limited-precision input (acos rejects unlimited precision) + let (re, _im) = c(0, 0).acos().into_parts(); let half_pi = F::from_parts(15707963267948966i64.into(), -16) .with_precision(60) .value(); @@ -323,8 +392,8 @@ mod tests { #[test] fn atan_one_is_quarter_pi() { use dashu_base::{Abs, AbsOrd}; - // atan(1) = π/4 - let (re, _im) = C::ONE.atan().into_parts(); + // atan(1) = π/4 (limited-precision input — atan rejects unlimited precision) + let (re, _im) = c(1, 0).atan().into_parts(); let quarter_pi = F::from_parts(7853981633974483i64.into(), -16) .with_precision(60) .value(); @@ -341,4 +410,19 @@ mod tests { let r = z.sin().asin(); assert!(r == z); } + + // The trig functions reject unlimited precision. `sin`/`cos`/`tan` do so via `guard`; the + // inverse trig (`asin`/`acos`/`atan`) build their work context directly and assert explicitly + // (like `powf`). The zero shortcuts (`C::ZERO.sin()` etc.) bypass the check, as they're exact. + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_sin_unlimited_panics() { + let _ = C::I.sin(); + } + + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_asin_unlimited_panics() { + let _ = C::ONE.asin(); + } } diff --git a/complex/src/misc.rs b/complex/src/misc.rs index 85531016..56dadf73 100644 --- a/complex/src/misc.rs +++ b/complex/src/misc.rs @@ -4,7 +4,7 @@ use crate::cbig::CBig; use crate::repr::{exact, CfpResult, Context}; use core::ops::Neg; use dashu_base::Sign; -use dashu_float::round::Round; +use dashu_float::round::{ErrorBounds, Round}; use dashu_float::{FBig, FpResult, Repr}; use dashu_int::Word; @@ -12,10 +12,6 @@ use dashu_int::Word; /// small fixed guard comfortably settles the accumulated rounding of two squarings and an add. const NORM_GUARD: usize = 8; -/// Guard digits (base-B) for `abs`. The inner `hypot` already carries its own guard; this extra -/// margin absorbs the final re-round to the CBig precision. -const ABS_GUARD: usize = 8; - impl CBig { /// The complex conjugate `x - iy`. Exact (sign flip of the imaginary part, including `-0`/`-inf`). #[inline] @@ -43,25 +39,6 @@ impl CBig { pub fn norm(&self) -> FBig { self.context.float().unwrap_fp(self.context.norm(self)) } - - /// The modulus `|z| = sqrt(re² + im²)` (a real [`FBig`]). A thin composition over - /// [`dashu_float::Context::hypot`] (the overflow-safe scaled sum-of-squares), evaluated at guard - /// precision and re-rounded. Near-correctly rounded. - /// - /// # Panics - /// - /// Panics if the precision is unlimited. - #[inline] - pub fn abs(&self) -> FBig { - self.context.float().unwrap_fp(self.context.abs(self)) - } - - /// The argument (phase) `atan2(im, re) ∈ ]-π, π]`. The branch cut lies on `]−∞, 0]`; signed zero - /// and infinities are handled per the C99 Annex G `atan2` table (reused from `dashu-float`). - #[inline] - pub fn arg(&self) -> FBig { - self.context.float().unwrap_fp(self.context.arg(self, None)) - } } impl Neg for CBig { @@ -139,12 +116,32 @@ impl Context { self.float(), ))); } - let gctx = self.guard(NORM_GUARD); + // `work_context` (not `guard`): at unlimited precision `sqr`/`add` are exact, so `norm` + // is too — it shouldn't panic or round there. + let gctx = self.work_context(NORM_GUARD); let re2 = gctx.unwrap_fp(gctx.sqr(z.re())); let im2 = gctx.unwrap_fp(gctx.sqr(z.im())); let n = gctx.unwrap_fp(gctx.add(re2.repr(), im2.repr())); Ok(n.with_precision(self.precision())) } +} + +impl Context { + /// The modulus `|z| = hypot(re, im)` (context layer), correctly rounded. Returns `+∞` for an + /// infinite input. `abs` is a single real operation, so it delegates directly to + /// [`dashu_float::Context::hypot`] (itself Ziv-correctly-rounded at the target precision) — + /// computing at `p + guard` and re-rounding would be a *double* rounding that can break correct + /// rounding. + /// + /// # Panics + /// + /// Panics if the precision is unlimited. + pub fn abs(&self, z: &CBig) -> FpResult> { + // Assert limited up front to keep the "transcendentals reject unlimited" contract — the + // float `hypot` would otherwise short-circuit unlimited to an exact value. + self.assert_limited(); + self.float().hypot(z.re(), z.im()) + } /// The argument `atan2(im, re)` (context layer). Delegates to `dashu-float`'s Annex-G `atan2`; /// the cache threads into it (the convenience layer passes `None`). @@ -155,17 +152,26 @@ impl Context { ) -> FpResult> { self.float().atan2(z.im(), z.re(), cache) } +} - /// The modulus `|z| = hypot(re, im)` (context layer). Near-correctly rounded; returns `+∞` for - /// an infinite input. Thin composition over [`dashu_float::Context::hypot`]. +impl CBig { + /// The modulus `|z| = sqrt(re² + im²)` (a real [`FBig`]). A thin composition over + /// [`dashu_float::Context::hypot`] (the overflow-safe scaled sum-of-squares, Ziv-correctly + /// rounded), evaluated at guard precision and re-rounded. /// /// # Panics /// /// Panics if the precision is unlimited. - pub fn abs(&self, z: &CBig) -> FpResult> { - let gctx = self.guard(ABS_GUARD); - let h = gctx.hypot(z.re(), z.im())?; - Ok(h.value().with_precision(self.precision())) + #[inline] + pub fn abs(&self) -> FBig { + self.context.float().unwrap_fp(self.context.abs(self)) + } + + /// The argument (phase) `atan2(im, re) ∈ ]-π, π]`. The branch cut lies on `]−∞, 0]`; signed zero + /// and infinities are handled per the C99 Annex G `atan2` table (reused from `dashu-float`). + #[inline] + pub fn arg(&self) -> FBig { + self.context.float().unwrap_fp(self.context.arg(self, None)) } } @@ -238,4 +244,11 @@ mod tests { assert!(a > F::ZERO); assert!(a < F::ONE); } + + // abs on an unlimited-precision CBig must panic, not silently compute |z| at ABS_GUARD digits. + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_abs_unlimited_precision_panics() { + let _ = C::ONE.abs(); + } } diff --git a/complex/src/mul.rs b/complex/src/mul.rs index c85b2a1a..6703627f 100644 --- a/complex/src/mul.rs +++ b/complex/src/mul.rs @@ -21,7 +21,7 @@ impl Context { if z.is_zero() { return Ok(exact(FBig::ZERO, FBig::ZERO)); } - let gctx = self.guard(MUL_GUARD); + let gctx = self.work_context(MUL_GUARD); let p = self.precision(); let (x, y) = (z.re(), z.im()); // real part: x² - y² @@ -43,7 +43,7 @@ impl Context { } return Ok(riemann(Context::max(z.context(), w.context()))); // ∞·finite = Riemann infinity } - let gctx = self.guard(MUL_GUARD); + let gctx = self.work_context(MUL_GUARD); let p = self.precision(); let (x, y) = (z.re(), z.im()); let (u, v) = (w.re(), w.im()); @@ -66,7 +66,7 @@ impl Context { } return Ok(riemann(*self)); } - let gctx = self.guard(MUL_GUARD); + let gctx = self.work_context(MUL_GUARD); let p = self.precision(); let re = gctx.mul(z.re(), s.repr())?.value().with_precision(p); let im = gctx.mul(z.im(), s.repr())?.value().with_precision(p); @@ -191,4 +191,12 @@ mod tests { let p2 = &s * &z; assert_eq!(p2.re().significand(), &6.into()); } + + // At unlimited precision `mul`/`sqr` are exact (`work_context` uses `self.float()`), so this + // must NOT panic (unlike the transcendentals). i·i = −1 exactly. + #[test] + fn mul_at_unlimited_is_exact() { + let i = C::I; + assert_eq!(&i * &i, C::NEG_ONE); + } } diff --git a/complex/src/repr.rs b/complex/src/repr.rs index 2b25d707..4b8d0aa7 100644 --- a/complex/src/repr.rs +++ b/complex/src/repr.rs @@ -64,17 +64,34 @@ impl Context { self.0.precision() } + /// Reject unlimited precision (0): [`Self::guard`] calls this (the guard-digit recipe is a + /// limited-precision technique), and `powf` calls it directly (it builds its work context + /// bypassing `guard`). Without it, `self.guard(G)` would silently make a finite `0 + G` context + /// and the transcendental would compute at ~`G` digits — the `dashu-float` layer only sees `G` + /// and can't catch it. Private counterpart of `dashu_float`'s internal check. + pub(crate) fn assert_limited(&self) { + if self.precision() == 0 { + panic!("precision cannot be 0 (unlimited) for this operation!") + } + } + /// The inner float context used to drive the real-part math (copied, since it is `Copy`). #[inline] pub(crate) const fn float(&self) -> FloatCtxt { self.0 } - /// Build a transient float working context at `p + g` guard digits — the guard-digit recipe - /// (§6.1 of the design doc) evaluates each component at extra precision and re-rounds to `p`. + /// The work context for an arithmetic op: a finite `p + g` guard-digit context at limited + /// precision, or the exact `self.float()` (precision 0) when the context is unlimited. So + /// `mul`/`sqr`/`norm` are exact at unlimited; `div`/`inv` still panic there via the float + /// layer's own check (a quotient isn't exactly representable in general). #[inline] - pub(crate) fn guard(&self, g: usize) -> FloatCtxt { - FloatCtxt::new(self.precision() + g) + pub(crate) fn work_context(&self, g: usize) -> FloatCtxt { + if self.precision() != 0 { + FloatCtxt::new(self.precision() + g) + } else { + self.float() + } } /// Unwrap a [`CfpResult`], returning the [`CBig`] value directly. diff --git a/complex/src/root.rs b/complex/src/root.rs index dcb2e041..486a6e98 100644 --- a/complex/src/root.rs +++ b/complex/src/root.rs @@ -3,8 +3,8 @@ use crate::cbig::CBig; use crate::repr::{combine_parts, exact, CfpResult, Context}; use dashu_base::Sign; -use dashu_float::round::Round; -use dashu_float::{FBig, Repr}; +use dashu_float::round::{ErrorBounds, Round}; +use dashu_float::{Context as FloatCtxt, FBig, Repr}; use dashu_int::Word; /// Guard digits (base-B) for `sqrt`. Composes `hypot` + two real `sqrt`s + adds; a modest fixed @@ -19,7 +19,7 @@ fn signed_inf(sign: Sign) -> Repr { } } -impl Context { +impl Context { /// Principal square root of a complex number (context layer). /// /// The result has non-negative real part; when the real part is zero the imaginary part is @@ -30,45 +30,50 @@ impl Context { return special; } - let gctx = self.guard(SQRT_GUARD); + // Principal sqrt via the cancellation-free form: for x ≥ 0, `a = sqrt((r+x)/2)`, + // `b = y/(2a)`; for x < 0, `b = sign(y)·sqrt((r-x)/2)`, `a = y/(2b)` — this avoids the + // near-cancellation in `r−x` when `|y| ≪ |x|`. The float `hypot`/`sqrt` are correctly-rounded + // at the working precision, and the adds/divs/mul each round at a few working-ULPs, so a + // small constant radius certifies both parts. The Ziv driver asserts a limited context (the + // special-value shortcut above is exact). let p = self.precision(); - let two = FBig::from_repr(Repr::new(2.into(), 0), gctx); - let x = z.re(); - let y = z.im(); - - // r = |z| (overflow-safe). Use the cancellation-free form: for x ≥ 0 compute `a` from - // `(r+x)/2` (large) and `b = y/(2a)`; for x < 0 compute `b` from `(r-x)/2` (large) and - // `a = y/(2b)`. This avoids subtracting nearly-equal magnitudes when |y| ≪ |x|. - let r = gctx.hypot(x, y)?.value(); - let (a, b) = if x.sign() != Sign::Negative { - // x ≥ 0 - let rpx = gctx.add(r.repr(), x)?.value(); - let half_rpx = gctx.div(rpx.repr(), two.repr())?.value(); - let a = gctx.sqrt(half_rpx.repr())?.value(); - let two_a = gctx.mul(two.repr(), a.repr())?.value(); - let b = gctx.div(y, two_a.repr())?.value(); - (a, b) - } else { - // x < 0: b carries the sign of y - let rmx = gctx.sub(r.repr(), x)?.value(); // r − x = r + |x| - let half_rmx = gctx.div(rmx.repr(), two.repr())?.value(); - let b_mag = gctx.sqrt(half_rmx.repr())?.value(); - let b = if y.sign() == Sign::Negative { - -b_mag + let [re, im] = self.ziv(SQRT_GUARD, |guard| { + let gctx = FloatCtxt::::new(p + guard); + let two = FBig::from_repr(Repr::new(2.into(), 0), gctx); + let x = z.re(); + let y = z.im(); + let r = gctx.hypot(x, y)?.value(); + let (a, b) = if x.sign() != Sign::Negative { + // x ≥ 0 + let rpx = gctx.add(r.repr(), x)?.value(); + let half_rpx = gctx.div(rpx.repr(), two.repr())?.value(); + let a = gctx.sqrt(half_rpx.repr())?.value(); + let two_a = gctx.mul(two.repr(), a.repr())?.value(); + let b = gctx.div(y, two_a.repr())?.value(); + (a, b) } else { - b_mag + // x < 0: b carries the sign of y + let rmx = gctx.sub(r.repr(), x)?.value(); // r − x = r + |x| + let half_rmx = gctx.div(rmx.repr(), two.repr())?.value(); + let b_mag = gctx.sqrt(half_rmx.repr())?.value(); + let b = if y.sign() == Sign::Negative { + -b_mag + } else { + b_mag + }; + let two_b = gctx.mul(two.repr(), b.repr())?.value(); + let a = gctx.div(y, two_b.repr())?.value(); + (a, b) }; - let two_b = gctx.mul(two.repr(), b.repr())?.value(); - let a = gctx.div(y, two_b.repr())?.value(); - (a, b) - }; - let re = a.with_precision(p); - let im = b.with_precision(p); + let a_rad = a.ulp() * 10; + let b_rad = b.ulp() * 10; + Ok([(a, a_rad), (b, b_rad)]) + })?; Ok(combine_parts(re, im)) } } -impl CBig { +impl CBig { /// Principal square root (convenience layer). /// /// # Panics @@ -199,4 +204,12 @@ mod tests { assert!(s.re().is_infinite()); assert!(s.im().is_pos_zero()); } + + // `sqrt` at unlimited precision panics via `guard` (the special-value shortcut above only + // catches zero/infinity, so a finite nonzero input reaches the guard context). + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn complex_sqrt_unlimited_panics() { + let _ = C::I.sqrt(); + } } diff --git a/complex/src/ziv.rs b/complex/src/ziv.rs new file mode 100644 index 00000000..f3ee217e --- /dev/null +++ b/complex/src/ziv.rs @@ -0,0 +1,262 @@ +//! The complex Ziv retry loop — guaranteed-correct rounding for complex transcendentals. +//! +//! Each complex transcendental (`exp`, `log`, …) approximates its result parts at a working +//! precision `p + guard`, together with a provable absolute error radius per part, and this driver +//! rounds each part to the target precision and checks that the whole per-part error interval +//! `[ã − E, ã + E]` lies inside the candidate's [`ErrorBounds`] preimage — the set of reals that +//! round to it. If every part fits, the rounded value is *guaranteed* correct; otherwise the driver +//! retries with more guard digits, sharing the recomputed working evaluation across all parts +//! (retrying while *any* part straddles a rounding boundary). The loop provably terminates (a true +//! tie is resolved deterministically by the mode), with a large sanity cap as an unreachable +//! backstop. +//! +//! The driver is generic in the part count `N`: most functions compute one complex result (`N = 2` +//! parts), while `sin_cos` computes two (`N = 4`, sharing one evaluation of the real `sin_cos` and +//! `sinh_cosh`). This mirrors `dashu-float`'s Ziv loop, with two complex-specific choices: +//! * The approximation closure is **fallible** — a complex composition can overflow mid-evaluation +//! (e.g. `exp` of a large real part), and propagating that [`FpError`] with `?` avoids hoisting +//! a per-function overflow probe out of the loop. Such errors are terminal (precision-independent) +//! and never trigger a retry. +//! * The per-part containment test runs on [`FBig`]s at unlimited precision +//! ([`FloatCtxt::new(0)`], where `+`/`−` are lossless), expressing `dashu-float`'s raw-`Repr` +//! interval arithmetic through the public `FBig` API so this crate need not reach into float's +//! internal arithmetic. The two are exactly equivalent. +//! +//! As in the float layer, each correctly-rounded float sub-operation contributes ~0 to the +//! composition radius (it is certified by its own Ziv loop); only the arithmetic composition steps +//! contribute, bounded by a small constant × working-ULP (plus any input-error amplification, which +//! each caller folds into its radius — e.g. `log`'s `ln|z|` near `|z| = 1`). + +use core::cmp::Ordering; + +use dashu_base::Approximation; +use dashu_float::round::{ErrorBounds, Rounding}; +use dashu_float::{Context as FloatCtxt, FBig, FpError}; +use dashu_int::Word; + +use crate::repr::Context; + +/// Maximum number of Ziv retries before falling back to the best-effort rounded value. +/// +/// A sanity backstop only — the loop converges as soon as the working precision is large enough +/// that no part's error interval straddles a rounding boundary, which happens in one attempt for +/// essentially all inputs (the guard-digit heuristic is sized for that) and in a handful of +/// attempts only for inputs pathologically close to a tie. The cap exists so a bug in an error +/// radius can never produce an infinite loop. +const MAX_ZIV_RETRIES: usize = 32; + +// A test-only retry counter (extra attempts beyond the first), mirroring `dashu-float`'s counter +// so tests can assert the loop converges on the first attempt for typical inputs — i.e. `0` means +// first-attempt success. +#[cfg(all(test, feature = "std"))] +thread_local! { + pub(crate) static LAST_ZIV_RETRIES: core::cell::Cell = const { core::cell::Cell::new(0) }; +} + +impl Context { + /// Correctly round a complex transcendental's `N` result parts to this context's precision using + /// a Ziv retry loop that certifies **every** part. + /// + /// `approx(guard)` evaluates the function at working precision `self.precision() + guard` and + /// returns `Ok([(v₀, e₀), …])` — each part's value and its provable absolute error radius, all + /// as [`FBig`]s at the working context — or an [`FpError`] (overflow / domain), which propagates + /// immediately without retry (such errors are terminal, not precision problems). The driver + /// rounds each part to the target precision and retries while *any* part's error interval + /// straddles a rounding boundary, sharing the recomputed evaluation. Returns each part's rounded + /// value (with its inexactness flag) for the caller to assemble. + /// + /// Rejects an unlimited context (the recipe is a limited-precision technique); each caller's + /// exact special-value shortcuts run before this driver and are unaffected. + pub(crate) fn ziv( + &self, + initial_guard: usize, + mut approx: impl FnMut(usize) -> Result<[(FBig, FBig); N], FpError>, + ) -> Result<[Approximation, Rounding>; N], FpError> { + self.assert_limited(); + let p = self.precision(); + let mut guard = initial_guard; + let mut last = None; + #[cfg(all(test, feature = "std"))] + LAST_ZIV_RETRIES.with(|c| c.set(0)); + for _ in 0..MAX_ZIV_RETRIES { + let parts = approx(guard)?; + // Round each part to the target precision. `with_precision` consumes the value, but the + // containment test still needs the working-precision original, so round clones. + let candidates = parts.clone().map(|(v, _)| v.with_precision(p)); + if parts + .iter() + .zip(candidates.iter()) + .all(|((v, e), c)| Self::contained::(v, e, c.value_ref())) + { + return Ok(candidates); + } + last = Some(candidates); + + // Grow the guard aggressively so a near-tie resolves in a couple of retries, while the + // first attempt (with the heuristic guard) handles the common case (matches float). + let step = core::cmp::max(guard, p / 2).max(1); + guard += step; + #[cfg(all(test, feature = "std"))] + LAST_ZIV_RETRIES.with(|c| c.set(c.get() + 1)); + } + + // Unreachable in practice: return the best-effort parts from the last attempt, matching the + // pre-Ziv near-correct behavior rather than looping forever. + Ok(last.expect("MAX_ZIV_RETRIES is non-zero")) + } + + /// Per-part containment test: is the approximation's error interval `[value ± radius]` + /// entirely inside the rounding preimage of `target` (every real in `[target − lb, target + rb]` + /// rounds to `target` under `R`)? + /// + /// The interval arithmetic runs on [`FBig`]s at unlimited precision ([`FloatCtxt::new(0)`]), + /// where addition is lossless — no rounding can drop a guard digit and mis-decide the call (a + /// wrong call here yields a wrong ULP). The sums are compared rather than the differences + /// (algebraically identical for exact arithmetic, reading as one shared inequality per endpoint): + /// `value − radius ≥ target − lb ⟺ value + lb ≥ target + radius` + /// `value + radius ≤ target + rb ⟺ target + rb ≥ value + radius` + fn contained( + value: &FBig, + radius: &FBig, + target: &FBig, + ) -> bool { + let (lb, rb, incl_l, incl_r) = R::error_bounds::(target); + let x = FloatCtxt::::new(0); + let left = x + .add(value.repr(), lb.repr()) + .unwrap() + .value() + .cmp(&x.add(target.repr(), radius.repr()).unwrap().value()); + let right = x + .add(target.repr(), rb.repr()) + .unwrap() + .value() + .cmp(&x.add(value.repr(), radius.repr()).unwrap().value()); + let left_ok = if incl_l { + left != Ordering::Less + } else { + left == Ordering::Greater + }; + let right_ok = if incl_r { + right != Ordering::Less + } else { + right == Ordering::Greater + }; + left_ok && right_ok + } +} + +#[cfg(all(test, feature = "std"))] +mod tests { + use super::*; + use crate::cbig::CBig; + use dashu_float::round::mode; + + type F = FBig; + + // Build a modest limited-precision complex input from f64 parts (away from poles/singularities). + fn z(re: f64, im: f64) -> CBig { + let mk = |v: f64| F::try_from(v).unwrap().with_precision(53).value(); + CBig::from_parts(mk(re), mk(im)) + } + + // The guard constants should let every migrated transcendental converge in ≤1 retry for typical + // inputs (Ziv certifies correctness; the guard only sizes the first-attempt hit rate). This + // catches gross guard mis-sizing, not the occasional near-tie retry. + #[test] + fn ziv_few_retries_for_typical_inputs() { + let ctx: Context = Context::new(53); + const MAX_RETRIES: usize = 1; + let cases = [ + z(0.5, 0.3), + z(1.0, -0.5), + z(2.0, 0.25), + z(0.25, 0.75), + z(1.5, 0.4), + ]; + for c in cases { + let assert_few = |name: &str, retries: usize| { + assert!( + retries <= MAX_RETRIES, + "{name} took {retries} retries (expected <= {MAX_RETRIES})" + ); + }; + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.exp(&c, None)); + assert_few("exp", LAST_ZIV_RETRIES.with(|cell| cell.get())); + + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.log(&c, None)); + assert_few("log", LAST_ZIV_RETRIES.with(|cell| cell.get())); + + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.sin_cos(&c, None)); + assert_few("sin_cos", LAST_ZIV_RETRIES.with(|cell| cell.get())); + + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.sqrt(&c)); + assert_few("sqrt", LAST_ZIV_RETRIES.with(|cell| cell.get())); + + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.asin(&c, None)); + assert_few("asin", LAST_ZIV_RETRIES.with(|cell| cell.get())); + + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.atan(&c, None)); + assert_few("atan", LAST_ZIV_RETRIES.with(|cell| cell.get())); + + LAST_ZIV_RETRIES.with(|cell| cell.set(usize::MAX)); + drop(ctx.powf(&c, &z(0.7, 0.2), None)); + assert_few("powf", LAST_ZIV_RETRIES.with(|cell| cell.get())); + } + } + + // An exact approximation (both radii 0) is accepted on the first attempt (0 retries). + #[test] + fn ziv_accepts_exact_first_attempt() { + let ctx: Context = Context::new(10); + LAST_ZIV_RETRIES.with(|c| c.set(usize::MAX)); + let r = ctx.ziv(4, |_| Ok([(F::ONE, F::ZERO), (F::from(2u8), F::ZERO)])); + assert!(r.is_ok()); + assert_eq!(LAST_ZIV_RETRIES.with(|c| c.get()), 0); + } + + // An approximation whose second part's error interval straddles a boundary must retry; once the + // shrinking radius makes the interval unambiguous, the driver accepts. + #[test] + fn ziv_retries_until_contained() { + let ctx: Context = Context::new(4); + let r = ctx.ziv(2, |guard| { + // value 1+1i, imaginary radius 2^(-guard): large on the first attempt, tiny later. + Ok([(F::ONE, F::ZERO), (F::ONE, F::ONE >> guard as isize)]) + }); + drop(r.unwrap()); + assert!(LAST_ZIV_RETRIES.with(|c| c.get()) >= 1); + } + + // The driver certifies all N parts together: here the 4-part form (as `sin_cos` uses) retries + // while any one part straddles a boundary. + #[test] + fn ziv_retries_until_all_four_contained() { + let ctx: Context = Context::new(4); + let r = ctx.ziv(2, |guard| { + let rad = F::ONE >> guard as isize; + Ok([ + (F::ONE, F::ZERO), + (F::from(2u8), F::ZERO), + (F::from(3u8), rad), // this part forces the retry + (F::from(4u8), F::ZERO), + ]) + }); + drop(r.unwrap()); + assert!(LAST_ZIV_RETRIES.with(|c| c.get()) >= 1); + } + + // The driver rejects unlimited precision (its callers' exact shortcuts run first). + #[test] + #[should_panic(expected = "precision cannot be 0")] + fn ziv_rejects_unlimited() { + let ctx: Context = Context::new(0); + drop(ctx.ziv(4, |_| Ok([(F::ONE, F::ZERO), (F::ONE, F::ZERO)]))); + } +} diff --git a/complex/tests/transcendental_prop.rs b/complex/tests/transcendental_prop.rs index 3030466f..7b8bd61b 100644 --- a/complex/tests/transcendental_prop.rs +++ b/complex/tests/transcendental_prop.rs @@ -32,6 +32,30 @@ fn small_strategy() -> impl Strategy { }) } +/// Modest real part (away from the poles `π/2 + kπ`) with a moderate imaginary part up to ~40 — +/// large enough that the naive `sin z / cos z` division would catastrophically cancel in the real +/// part, small enough that `cosh(2·im)` stays cheap. For the cancellation-free `tan` oracle. +fn tan_strategy() -> impl Strategy { + ((1i64..(1i64 << 20), -4isize..4isize), 1i64..40i64).prop_map(|((re_sig, re_exp), im_num)| { + let re = F::from_parts(re_sig.into(), re_exp) + .with_precision(P) + .value(); + let im = F::from_parts(im_num.into(), 0).with_precision(P).value(); + CBig::from_parts(re, im) + }) +} + +/// Real part in roughly `(0.9, 1.1)` (straddling the `asin`/`acos` singularity `z = 1`) with a tiny +/// imaginary part — exercises the factored `1-z² = (1-z)(1+z)` form right at `z = ±1`, where the +/// direct `1 - z²` would catastrophically cancel against the `sqr` rounding error. +fn near_one_strategy() -> impl Strategy { + (90i64..110, -2i64..2).prop_map(|(re_num, im_num)| { + let re = F::from_parts(re_num.into(), -2).with_precision(P).value(); + let im = F::from_parts(im_num.into(), -8).with_precision(P).value(); + CBig::from_parts(re, im) + }) +} + fn within_ulps(a: &F, b: &F, k: u32) -> bool { if a == b { return true; @@ -64,7 +88,19 @@ proptest! { let hi = Context::new(2 * P); let rp = lo.sqrt(&z).unwrap().value(); let r2 = reround_hi(hi.sqrt(&z).unwrap().value()); - prop_assert!(within_ulps_cbig(&rp, &r2, 2)); + // sqrt is correctly rounded via Ziv, so 1 ULP. + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); + } + + #[test] + fn tan_self_oracle(z in tan_strategy()) { + // tan uses the cancellation-free double-angle form, so it stays correct for moderate-large + // |Im z| (where the naive sin/cos division would cancel). 1 ULP. + let lo = Context::new(P); + let hi = Context::new(2 * P); + let rp = lo.tan(&z, None).unwrap().value(); + let r2 = reround_hi(hi.tan(&z, None).unwrap().value()); + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); } #[test] @@ -73,7 +109,9 @@ proptest! { let hi = Context::new(2 * P); let rp = lo.exp(&z, None).unwrap().value(); let r2 = reround_hi(hi.exp(&z, None).unwrap().value()); - prop_assert!(within_ulps_cbig(&rp, &r2, 2)); + // exp is correctly rounded via Ziv, so the low- and (re-rounded) high-precision results + // agree to within 1 ULP (they match exactly except at rare near-tie double-rounding cases). + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); } #[test] @@ -82,7 +120,9 @@ proptest! { let hi = Context::new(2 * P); let rp = lo.log(&z, None).unwrap().value(); let r2 = reround_hi(hi.log(&z, None).unwrap().value()); - prop_assert!(within_ulps_cbig(&rp, &r2, 2)); + // log is correctly rounded via Ziv (its real part even carries an amplification term near + // |z|=1), so the results agree to within 1 ULP. + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); } #[test] @@ -104,8 +144,20 @@ proptest! { let cp = lo.cos(&z, None).unwrap().value(); let s2 = reround_hi(hi.sin(&z, None).unwrap().value()); let c2 = reround_hi(hi.cos(&z, None).unwrap().value()); - prop_assert!(within_ulps_cbig(&sp, &s2, 2)); - prop_assert!(within_ulps_cbig(&cp, &c2, 2)); + // sin/cos are correctly rounded via a shared Ziv loop, so 1 ULP. + prop_assert!(within_ulps_cbig(&sp, &s2, 1)); + prop_assert!(within_ulps_cbig(&cp, &c2, 1)); + } + + #[test] + fn abs_self_oracle(z in cbig_strategy()) { + // abs delegates directly to the (Ziv-correctly-rounded) real hypot, so abs@p agrees with + // abs@2p re-rounded to p exactly (0 ULP) save for rare near-tie double-rounding (≤1 ULP). + let lo = Context::new(P); + let hi = Context::new(2 * P); + let lp = lo.abs(&z).unwrap().value(); + let hp = hi.abs(&z).unwrap().value().with_precision(P).value(); + prop_assert!(within_ulps(&lp, &hp, 1)); } #[test] @@ -129,7 +181,8 @@ proptest! { let hi = Context::new(2 * P); let rp = lo.asin(&z, None).unwrap().value(); let r2 = reround_hi(hi.asin(&z, None).unwrap().value()); - prop_assert!(within_ulps_cbig(&rp, &r2, 4)); + // asin is correctly rounded via Ziv (well-conditioned regime), 1 ULP. + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); } #[test] @@ -138,6 +191,43 @@ proptest! { let hi = Context::new(2 * P); let rp = lo.atan(&z, None).unwrap().value(); let r2 = reround_hi(hi.atan(&z, None).unwrap().value()); - prop_assert!(within_ulps_cbig(&rp, &r2, 4)); + // atan is correctly rounded via Ziv (well-conditioned regime), 1 ULP. + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); + } + + #[test] + fn powf_self_oracle((base, w) in (small_strategy(), small_strategy())) { + // modest base/exponent keep `w·log base` bounded (no overflow); base has positive real part + // so `log` stays off its branch cut. powf is correctly rounded via Ziv, 1 ULP. + let lo = Context::new(P); + let hi = Context::new(2 * P); + let rp = lo.powf(&base, &w, None).unwrap().value(); + let r2 = reround_hi(hi.powf(&base, &w, None).unwrap().value()); + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); + } + + #[test] + fn powi_self_oracle(z in small_strategy(), n in -8i32..=8) { + prop_assume!(n != 0); + let lo = Context::new(P); + let hi = Context::new(2 * P); + let rp = lo.powi(&z, n.into()).unwrap().value(); + let r2 = reround_hi(hi.powi(&z, n.into()).unwrap().value()); + // powi is correctly rounded via Ziv (positive and negative integer exponents), 1 ULP. + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); + } + + #[test] + fn asin_acos_near_one_self_oracle(z in near_one_strategy()) { + // z straddles the singularity at 1 (and the symmetric point -1 by conjugation), where the + // factored `1-z²=(1-z)(1+z)` keeps `asin`/`acos` accurate; the self-oracle checks 1 ULP. + let lo = Context::new(P); + let hi = Context::new(2 * P); + let rp = lo.asin(&z, None).unwrap().value(); + let r2 = reround_hi(hi.asin(&z, None).unwrap().value()); + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); + let rp = lo.acos(&z, None).unwrap().value(); + let r2 = reround_hi(hi.acos(&z, None).unwrap().value()); + prop_assert!(within_ulps_cbig(&rp, &r2, 1)); } } diff --git a/float/CHANGELOG.md b/float/CHANGELOG.md index 29b41346..b2fe3902 100644 --- a/float/CHANGELOG.md +++ b/float/CHANGELOG.md @@ -6,9 +6,104 @@ - `Repr::new_const`: a `const`-evaluable, normalized `Repr` constructor from a `DoubleWord` significand (the `const` counterpart of `Repr::new`). `FBig::from_parts_const` now delegates to it, and the complex literal macro uses it. -- `FBig::log2` / `Context::log2` / `CachedFBig::log2`: base-2 logarithm, correctly rounded via - `ln(x)/ln(2)` evaluated at an elevated working precision. Previously only the f32-precision - `log2_bounds` magnitude estimate was available, so directed `log2` was wrong by many ULPs. +- `Repr::is_int` is now `const` (it only inspects the exponent and the infinity sentinel). +- **Exact `Add`/`Sub`/`Mul` operators for `Repr`** (in a new `repr_ops` module). A `Repr` carries no + precision limit, so add/sub/mul on it are lossless — these are the shared primitives the crate uses + for exact intermediates (the Ziv containment test, the correctly-rounded `Sum`, and the `FBig` + multiply path). `Mul` saturates exponent overflow/underflow to the signed infinity/zero sentinels, + so the operator is infallible; the internal `make_mul_repr`/`unwrap_mul_repr!` helpers are removed + in favor of the operator. No behavior change to `FBig` arithmetic. +- **Guaranteed-correct rounding for `exp`, `exp_m1`, `ln`, `ln_1p`** via a Ziv retry loop. A generic + `Context::ziv` driver rounds the working-precision approximation to the target precision and + verifies, against the `ErrorBounds` rounding preimage, that the approximation's provable error + interval lies entirely inside one rounding bin; if not, it retries with more guard digits. The + loop provably terminates and preserves the `Exact`/`Inexact` flag. Series evaluation is factored + into near-correct `ln_compute`/`exp_compute` cores (`R: Round`) that the Ziv wrapper certifies. +- **Tightened `exp` guard digits** (now a performance knob, since Ziv — not the guard count — + guarantees correctness): the `Bⁿ`-powering guard is halved from `2n` to `n` and the series guard + drops its conservative `+ 2`. +- **Guaranteed-correct rounding for most remaining transcendentals** via the Ziv loop: the + trigonometric family (`sin`, `cos`, `sin_cos`, `tan`, `asin`, `acos`, `atan`, `atan2`), the + hyperbolic family (`sinh`, `cosh`, `sinh_cosh`, `tanh`, `asinh`, `acosh`, `atanh`), and `hypot`. + The trig series (`sin`/`cos`/`atan`) + are factored into near-correct `_compute` cores like `exp_compute`; the composition-based functions + treat the now-Ziv-correct `exp`/`ln`/`atan` as black boxes and count only their arithmetic. The + trig argument reduction folds a `|k|·ulp(π/2)` reduction-error term into the radius so the + containment test stays sound for huge `|x|`. `ziv_pair` certifies both halves of `sin_cos`/ + `sinh_cosh`. +- **Guaranteed-correct rounding for `powf` (non-integer exponent)** via the Ziv loop. `x^y = + exp(y·ln x)`, and `exp` amplifies the rounding of `y·ln x` by the result magnitude — so the + radius is `result.ulp()·(|y·ln x|+1)·(B+8)` taken at the *working* precision: it shrinks as + `B^{-guard}`, so the containment test converges (a radius computed at unlimited precision would + be constant across retries and never settle for a value near a rounding boundary). An + integer-valued exponent delegates to `powi` (binary exponentiation), which also admits a + negative base — its sign fixed by the exponent's parity — so `powf(-x, n)` is in domain for + integer `n`. +- **Guaranteed-correct rounding for `powi` (integer exponent)** via the Ziv loop. Binary + exponentiation (repeated squaring) compounds the relative error — it roughly doubles per + squaring — so after `bit_len(n)` squarings the error is bounded by about `2^nlen · ulp`, which + the Ziv radius reflects (`ulp_w << (nlen + 1)`). A negative exponent computes `(1/base)^|n|` + directly, so the sign-dependent overflow/underflow falls out naturally. When the squaring chain + rounds `Exact` (the result is exactly representable, e.g. an integer power that fits), the true + error is zero and the radius is reported as zero — this is required under the *directed* rounding + modes (`Zero`/`Down`/`Up`/`Away`, the `FBig` default), where an exactly-representable result lies + on a one-sided rounding boundary that no nonzero radius can fit inside. +- **`FBig::log2` / `Context::log2` / `CachedFBig::log2`**: base-2 logarithm, correctly rounded via + the Ziv loop. `log2(x) = ln(x)/ln(2)`, and rounding the two operands separately before the single + division is only near-correct — under a directed mode, enlarging the positive denominator shrinks + the quotient, so the directed result is unbounded — so each `ln_compute` reports its provable + radius and the two radii are carried through the division as an outward-rounded interval `[lo, hi]` + the Ziv driver certifies. An exact power of two short-circuits to the integer `log2(x)`; this is + required under directed modes, where the exactly-representable integer sits on a one-sided rounding + boundary no positive radius can fit inside (so `log2(2^-159)` under `Up` is the exact `-159`, not + `-159 + 1 ulp`). + +### Change +- **(breaking, bound)** `Context::exp`/`exp_m1`/`ln`/`ln_1p` (and `powf`, the hyperbolic family, + and the base-conversion path that derives from them) now require `R: ErrorBounds` rather than + `R: Round` — the Ziv containment test needs the rounding preimage that `ErrorBounds` provides. + All six built-in modes satisfy `ErrorBounds`; only custom non-`ErrorBounds` `Round` modes are + affected (custom modes are already discouraged). `CachedFBig`/`CachedCBig` forwarders for these + methods carry the same bound. Arithmetic, roots, and trigonometric remain `R: Round`; `powi` + now also requires `R: ErrorBounds` (it is correctly rounded via the Ziv loop). +- **Internal dedup** (no behavior change): the `⌈log_B(precision)⌉` base-guard formula shared by every + transcendental Ziv loop is now `Context::base_guard_digits::()` (12 call sites in `hyper`/`exp`/ + `log`), and the `ulp·(4·terms + 12)` series-truncation radius is now `series_radius(value, terms)` + (shared by the `sin`/`cos`/`sin_cos`/`atan`/`ln` series cores). The two textually-identical + `CachedFBig` forwarding macros (`forward_to_context!` / `forward_to_context_unwrap!`) are merged. + +### Add +- **`hypot` is now correctly rounded via the Ziv loop** with MPFR-style exactness tracking. The + closure computes `sqrt(large² + small²)` (operands scaled to avoid `large²` overflowing the + exponent — no division) and OR's each step's `Exact`/`Inexact` flag, mirroring MPFR's `exact` + flag (`hypot.c`): an all-exact chain yields the exact true value, reported to `ziv` with radius 0, + which it accepts without the containment test. This is what lets an exact Pythagorean-triple + result (e.g. `hypot(3,4)=5`, `hypot(5,12)=13`) terminate under a directed rounding mode, where + the one-sided preimage would otherwise make the containment test infinite-retry. + +### Fix +- **Directed-rounding hang on exact results** (`acos(1)`, `acos(-1)`, `asin(0)`, `hypot` of a + Pythagorean triple): under a directed mode (Down/Up/Zero) a function whose true value is exactly + representable carries a positive radius that can't be certified against the value's one-sided + rounding preimage, so the Ziv containment test infinite-retried (and the retry eventually tripped + a `dashu-int` NTT assertion, now fixed). `acos` short-circuits `|x| = 1` to the exact `0` / `π`, + `asin` short-circuits `x = 0` to `±0` (matching the existing special cases in the rest of the + inverse trig/hyperbolic family), and `hypot` reports radius 0 when its computation chain is exact + (see the new `hypot` entry above). +- **`tan` no longer has a hoisted pole check.** The check (which tested `cos` with `is_pos_zero`, + missing `-0`, and computed a `±∞` pole sign) re-evaluated the sin/cos series a second time on top of + the first Ziv attempt — a ~2× cost on every `tan` call. It's removed: near a pole (an odd multiple + of π/2) the value is large but finite, and dashu's wide exponent range holds it as a finite number + whose sign is carried by the arithmetic (`s/−|c|` is negative), so no `±∞` special-case is needed. + The unreachable exact-pole case (cos cancelling to a zero significand — impossible for finite- + precision input, since a `p`-digit rational can't sit closer than ~`B⁻ᵖ` to the irrational pole) is + handled by the Ziv-closure's `significand.is_zero()` retry guard, which forces a higher guard where + cos is representable. +- **`no_std` build of the test-only Ziv retry counter.** The `LAST_ZIV_RETRIES` `thread_local!` + (used by the retry-count tests) requires `std`, so the crate failed to compile under + `--no-default-features` (the `thread_local!` macro isn't in scope). It's now gated behind + `feature = "std"` along with its uses and the counter-reading tests, so the Ziv driver itself is + `no_std`-clean; the retry-count tests run under `std` as before. ### Fix - `to_f64`/`to_f32` now round the source once, directly to the target's precision at its own diff --git a/float/src/convert.rs b/float/src/convert.rs index a73d1e1e..9ae7851e 100644 --- a/float/src/convert.rs +++ b/float/src/convert.rs @@ -666,11 +666,19 @@ impl Context { // if the exponent is large, then we first estimate the result exponent as floor(exponent * log(B) / log(NewB)), // then the fractional part is multiplied with the original significand let work_context = Context::::new(2 * self.precision); // double the precision to get the precise logarithm + // ln(old base) and ln(new base) — near-correct is sufficient for the exponent estimate, + // and using the near-correct `ln_compute`/`ln_base` (R: Round) keeps base conversion + // off the `ErrorBounds` bound. Both are computed in base NewB so the euclidean division + // has matching bases. let new_exp = repr.exponent - * work_context.unwrap_fp( - work_context - .ln(&Repr::new(Repr::::BASE.into(), 0), reborrow_cache(&mut cache)), - ); + * work_context + .ln_compute::( + &Repr::new(Repr::::BASE.into(), 0), + work_context.precision, + false, + reborrow_cache(&mut cache), + ) + .0; let (exponent, rem) = new_exp.div_rem_euclid(work_context.ln_base::(reborrow_cache(&mut cache))); let exponent_sign = exponent.sign(); @@ -683,7 +691,18 @@ impl Context { ); } }; - let exp_rem = rem.exp(); + // exp(fractional exponent) — near-correct is sufficient (it scales the significand), + // so use `exp_compute` (R: Round) and stay off the `ErrorBounds` bound. + let n = 1usize << (work_context.precision.bit_len() / 2); + let exp_rem = work_context + .exp_compute::( + &rem.repr, + work_context.precision, + false, + n, + reborrow_cache(&mut cache), + ) + .0; let significand = repr.significand * exp_rem.repr.significand; let repr = Repr::new(significand, exponent + exp_rem.repr.exponent); self.repr_round(repr) diff --git a/float/src/exp.rs b/float/src/exp.rs index 54948119..3961b09d 100644 --- a/float/src/exp.rs +++ b/float/src/exp.rs @@ -5,13 +5,14 @@ use crate::{ fbig::FBig, math::cache::{reborrow_cache, ConstCache}, repr::{Context, Repr, Word}, - round::Round, - utils::ceil_usize, + round::{ErrorBounds, Round}, }; -use dashu_base::{AbsOrd, Approximation::*, BitTest, DivRemEuclid, EstimatedLog2, Sign}; -use dashu_int::IBig; +use dashu_base::{Abs, AbsOrd, Approximation::*, BitTest, DivRemEuclid, EstimatedLog2, Sign}; +use dashu_int::{IBig, UBig}; -impl FBig { +// `powi` (integer power), `powf`/`exp`/`exp_m1` route through Ziv-backed Context methods, which +// require `R: ErrorBounds` for their correctness guarantee. +impl FBig { /// Raise the floating point number to an integer power. /// /// # Examples @@ -84,10 +85,178 @@ impl FBig { } } -// TODO: give the exact formulation of required guard bits - impl Context { - /// Raise the floating point number to an integer power under this context. + /// Left-to-right binary exponentiation of `start` to the power `n` (`n ≥ 2`) at this context's + /// precision — the shared squaring kernel. + /// + /// Each `sqr`/`mul` is correctly rounded, but their rounding flags are folded away (`.value()`) + /// and no containment test is applied, so the result is only *near*-correct: repeated squaring + /// compounds the relative error (it roughly doubles per step), so after `n.bit_len()` squarings + /// the error is on the order of `2^nlen · ulp`. The public [`powi`](Context::powi) retries this + /// kernel inside a Ziv loop to certify the rounding; `exp_compute` also uses it for its internal + /// `Bⁿ` powering (where the outer `exp` Ziv loop absorbs the error). + /// + /// Returns the value together with an `exact` flag that is `true` only when **every** squaring + /// and multiplication rounded `Exact` (so the returned value is the mathematically exact + /// `startⁿ`). The Ziv caller uses this to report a zero radius for exact results — under + /// directed rounding modes an exactly-representable result sits on a one-sided rounding + /// boundary, which a nonzero radius can never certify. + pub(crate) fn powi_chain( + &self, + start: &Repr, + n: &UBig, + ) -> (FBig, bool) { + let nlen = n.bit_len(); + debug_assert!(nlen >= 2, "powi_chain requires n >= 2"); + let mut p = nlen - 2; + let first = self.sqr(start); + let mut exact = matches!(first, Ok(Exact(_))); + let mut res = self.unwrap_fp(first); + loop { + if n.bit(p) { + let m = self.mul(res.repr(), start); + exact = exact && matches!(m, Ok(Exact(_))); + res = self.unwrap_fp(m); + } + if p == 0 { + break; + } + p -= 1; + let s = self.sqr(res.repr()); + exact = exact && matches!(s, Ok(Exact(_))); + res = self.unwrap_fp(s); + } + (res, exact) + } + + /// Near-correct exp core: evaluate `exp(x)` (or `exp_m1(x)` when `minus_one`) at + /// `work_precision`, returning `(value, error_radius)`. + /// + /// Shared by the Ziv-backed `exp`/`exp_m1` (which retry it) and usable directly where only a + /// near-correct value is needed. The caller must have pre-checked that the reduction quotient + /// `s = floor(x/ln B)` fits `isize` (astronomical `|x|` overflows and is handled before the + /// Ziv loop, since this closure can't return `Err`). `n` (the reduction power, `≈ √p`) is + /// derived from the *target* precision and is constant across retries. + pub(crate) fn exp_compute( + &self, + x: &Repr, + work_precision: usize, + minus_one: bool, + n: usize, + mut cache: Option<&mut ConstCache>, + ) -> (FBig, FBig) { + // exp(x) = B^s · exp(r)^(Bⁿ), with r = x − s·ln(B) reduced so |r| < B⁻ⁿ. + let context = Context::::new(work_precision); + let x = FBig::new(context.repr_round_ref(x).value(), context); + + // When minus_one is true and |x| < 1/B, evaluate the Maclaurin series without scaling + // (no Bⁿ reduction, no powering — n_eff = 0). + let no_scaling = minus_one && x.log2_est() < -B.log2_est(); + + let (s, r, n_eff) = if no_scaling { + (0isize, x, 0usize) + } else { + let logb = context.ln_base::(reborrow_cache(&mut cache)); + let (s_big, r) = x.div_rem_euclid(logb); + let s: isize = s_big + .try_into() + .expect("exp reduction quotient fits isize (overflow pre-checked)"); + (s, r, n) + }; + let r = r >> n_eff as isize; + + // Maclaurin series: exp(r) = 1 + Σ rⁱ/i! + let mut factorial = IBig::ONE; + let mut pow = r.clone(); + let mut sum = if no_scaling { + r.clone() + } else { + FBig::ONE + &r + }; + let mut k = 2u32; + let mut terms: usize = 1; + loop { + factorial *= k; + pow *= &r; + + let increase = &pow / &factorial; + if increase.abs_cmp(&sum.sub_ulp()).is_le() { + break; + } + sum += increase; + k += 1; + terms += 1; + } + + // The radius is computed at *unlimited* precision so the bound arithmetic is exact — a + // work-precision product would drop digits and could under-estimate (a soundness hole). + let ulp_w = || sum.ulp().with_precision(0).value(); + + if no_scaling { + // exp_m1(x) = sum directly; error is the series truncation + rounding. + let radius = ulp_w() * (4 * terms + 8) + ulp_w(); + (sum, radius) + } else { + // Powering amplifies the series' relative error by Bⁿ. With |v|/|sum| < e < 3 (both + // near 1, since |r| < B⁻ⁿ), |v − true| ≤ 3·Bⁿ·(4K+8)·ulp(sum) + ulp(v). The B^s shift + // is exact, so the bound shifts with the value. + // + // The squaring chain compounds the relative error (it doubles per step), so it is run + // at an inflated precision and rounded back to `work_precision` — near-correct, which + // is what the `+ ulp(v)` term above accounts for. + let bn: UBig = Repr::::BASE.pow(n); + let chain_ctx = + Context::::new(work_precision + bn.bit_len() + work_precision.bit_len()); + let (v_pow, _) = chain_ctx.powi_chain(sum.repr(), &bn); + let v = v_pow.with_precision(work_precision).value(); + let v_shifted = v.clone() << s; + let e_v = (ulp_w() << n as isize) * (4 * terms + 8) * 3u32 + + v.ulp().with_precision(0).value(); + let radius = if minus_one { + // result = v_shifted − 1; the subtraction adds one result-ULP of rounding. + let result = &v_shifted - FBig::ONE; + let radius = (e_v << s) + result.ulp().with_precision(0).value(); + return (result, radius); + } else { + e_v << s + }; + (v_shifted, radius) + } + } +} + +/// Hoisted `exp` overflow probe for the Ziv closures (which can't return `Err`). Returns `true` +/// when `exp(x)` is outside the finite exponent range — astronomically large `|x|` (the reduction +/// quotient `s = x/ln B` overflows `isize`). True for both signs of huge `x` (the quotient +/// *magnitude* overflows `isize`). Shared by `exp_internal`, `powf`, and the hyperbolic functions. +pub(crate) fn exp_overflows( + ctx: &Context, + x: &Repr, + cache: &mut Option<&mut ConstCache>, +) -> bool { + if x.log2_est().abs() <= 61.0 { + return false; + } + let probe = Context::::new(ctx.precision + 64); + let logb = probe.ln_base::(reborrow_cache(cache)); + let x_probe = FBig::new(probe.repr_round_ref(x).value(), probe); + let s_probe = x_probe.div_rem_euclid(logb).0; + >::try_from(s_probe).is_err() +} + +// `powi` (integer power), `powf` (non-integer exponent), `exp`, and `exp_m1` are correctly rounded +// via the Ziv loop, so they require `R: ErrorBounds`. `powf` with an integer-valued exponent +// delegates to `powi`. +impl Context { + /// Raise the floating point number to an integer power under this context, correctly rounded + /// via a Ziv retry loop. + /// + /// `base^n` is computed by left-to-right binary exponentiation (repeated squaring); a negative + /// exponent computes `(1/base)^|n|`, so the sign-dependent overflow/underflow falls out + /// naturally. Each squaring compounds the relative error (it roughly doubles per step), so + /// after `n.bit_len()` squarings the error is bounded by about `2^nlen · ulp` — the Ziv radius + /// reflects that, and the loop retries with more guard digits until the working-precision + /// interval unambiguously determines the target rounding. /// /// # Examples /// @@ -106,8 +275,8 @@ impl Context { /// /// # Panics /// - /// Panics if the precision is unlimited and the exponent is negative. In this case, the exact - /// result is likely to have infinite digits. + /// Panics if the precision is unlimited and the exponent is negative (the exact `1/base` is + /// not finite in general). pub fn powi(&self, base: &Repr, exp: IBig) -> FpResult> { // TODO: range handling has three known limitations at the exponent extremes: // (1) the overflow guard below estimates the result magnitude with an f64 and misclassifies @@ -117,83 +286,129 @@ impl Context { if base.is_infinite() { return Err(FpError::InfiniteInput); } - - let (exp_sign, exp) = exp.into_parts(); - if exp_sign == Sign::Negative { - // if the exponent is negative, then negate the exponent - // note that do the inverse at last requires less guard bits - assert_limited_precision(self.precision); // TODO: we can allow this if the inverse is exact (only when significand is one?) - - let guard_bits = self.precision.bit_len() * 2; // heuristic - let rev_context = Context::::new(self.precision + guard_bits); - let pow = rev_context.unwrap_fp(rev_context.powi(base, exp.into())); - let inv = rev_context.unwrap_fp_repr(rev_context.repr_div(Repr::one(), pow.repr)); - let repr = self.repr_round(inv); - return Ok(repr.map(|v| FBig::new(v, *self))); + let (exp_sign, n) = exp.into_parts(); + let negative = exp_sign == Sign::Negative; + if negative { + // a negative exponent needs 1/base, which is not finite at unlimited precision + assert_limited_precision(self.precision); } - if exp.is_zero() { + + if n.is_zero() { return Ok(Exact(FBig::ONE)); - } else if exp.is_one() { + } + if n.is_one() { + if negative { + // base^(-1) = 1/base: a single correctly-rounded division + return self.div(&Repr::one(), base); + } let repr = self.repr_round_ref(base); return Ok(repr.map(|v| FBig::new(v, *self))); } - // Guard against exponent overflow for astronomically large results: the result - // magnitude has log2 ≈ exp·log2(base); if that exceeds the isize exponent range, - // return ±inf (|base| > 1) or 0 (|base| < 1) instead of overflowing mid-computation. - let base_log2 = base.log2_est() as f64; - let threshold = (isize::MAX as f64) * (B.log2_est() as f64); - let exp_f64 = i64::try_from(&exp).ok().map(|e| e as f64); - let overflows = match exp_f64 { - Some(e) => e * base_log2 > threshold, - None => base_log2 != 0.0, // exp doesn't fit i64: overflows unless |base| == 1 - }; - if overflows { - return if base_log2 > 0.0 { - Err(FpError::Overflow(if base.sign() == Sign::Negative { - Sign::Negative - } else { - Sign::Positive - })) - } else { - // |base| < 1 and exponent huge → underflow to signed zero - let underflow_sign = if base.sign() == Sign::Negative && exp.bit(0) { + // Zero base (±0): a positive exponent gives ±0, a negative one ±inf; the sign follows |n|'s + // parity. Short-circuit before the magnitude pre-check (whose log2 estimate is meaningless + // for zero) and before the squaring chain (which can't start from zero). + let odd = n.bit(0); + if base.significand.is_zero() { + let neg_sign = base.sign() == Sign::Negative && odd; + if negative { + let sign = if neg_sign { Sign::Negative } else { Sign::Positive }; - Err(FpError::Underflow(underflow_sign)) + return Ok(Exact(FBig::new(Repr::::infinity_with_sign(sign), *self))); + } + let repr = if neg_sign { + Repr::::neg_zero() + } else { + Repr::::zero() }; + return Ok(Exact(FBig::new(repr, *self))); } - let work_context = if self.is_limited() { - // increase working precision when the exponent is large - let guard_digits = exp.bit_len() + self.precision.bit_len(); // heuristic - Context::::new(self.precision + guard_digits) - } else { - Context::::new(0) - }; - - // binary exponentiation from left to right - let mut p = exp.bit_len() - 2; - let mut res = work_context.unwrap_fp(work_context.sqr(base)); - loop { - if exp.bit(p) { - res = work_context.unwrap_fp(work_context.mul(res.repr(), base)); + // Magnitude pre-check: the result's log2 is `signed_exp · log2|base|`; outside the finite + // exponent range it short-circuits to overflow/underflow instead of letting the squaring + // chain overflow mid-computation (the Ziv closure below can't return `Err`). + let base_log2 = base.log2_est() as f64; + let threshold = (isize::MAX as f64) * (B.log2_est() as f64); + let result_log2 = match i64::try_from(&n).ok() { + Some(e) => { + let signed = if negative { -(e as f64) } else { e as f64 }; + signed * base_log2 } - if p == 0 { - break; + None => { + // |n| doesn't fit i64: the magnitude is unbounded. |base| ≈ 1 (estimate 0) gives + // exactly ±1 regardless of the huge exponent; otherwise it over- or underflows. + if base_log2 == 0.0 { + let repr = if base.sign() == Sign::Negative && odd { + Repr::::neg_one() + } else { + Repr::::one() + }; + return Ok(Exact(FBig::new(repr, *self))); + } + let over = (!negative && base_log2 > 0.0) || (negative && base_log2 < 0.0); + let sign = if base.sign() == Sign::Negative && odd { + Sign::Negative + } else { + Sign::Positive + }; + return Err(if over { + FpError::Overflow(sign) + } else { + FpError::Underflow(sign) + }); } - p -= 1; - res = work_context.unwrap_fp(work_context.sqr(res.repr())); + }; + if result_log2 > threshold || result_log2 < -threshold { + let sign = if base.sign() == Sign::Negative && odd { + Sign::Negative + } else { + Sign::Positive + }; + return Err(if result_log2 > threshold { + FpError::Overflow(sign) + } else { + FpError::Underflow(sign) + }); } - Ok(res.with_precision(self.precision)) + let nlen = n.bit_len(); + let initial_guard = nlen + self.base_guard_digits::() + 2; + Ok(self.ziv(initial_guard, |guard| { + let pw = self.precision + guard; + let work = Context::::new(pw); + // start from base (positive exponent, always exact) or its working-precision + // reciprocal (negative exponent, exact only when 1/base is exactly representable). + let (start, start_exact) = if negative { + let d = work.div(&Repr::one(), base); + let exact = matches!(d, Ok(Exact(_))); + (work.unwrap_fp(d).repr().clone(), exact) + } else { + (base.clone(), true) + }; + let (res, chain_exact) = work.powi_chain(&start, &n); + // When the whole computation is exact (start exact + no squaring rounded), `res` is the + // exact value and the true error is 0 — report a zero radius. This is required under + // directed rounding modes, where an exactly-representable result lies on a one-sided + // rounding boundary that no nonzero radius can fit inside (the Ziv loop would retry + // forever). Otherwise the squaring compounds the error ~`2^nlen · ulp_w`. + let radius = if pw == 0 || (start_exact && chain_exact) { + FBig::ZERO + } else { + res.ulp().with_precision(0).value() << (nlen as isize + 1) + }; + (res, radius) + })) } /// Raise the floating point number to an floating point power under this context. /// - /// Note that this method will not rely on [FBig::powi] even if the `exp` is actually an integer. + /// A non-integer exponent is correctly rounded via a Ziv loop. An integer-valued exponent + /// delegates to [`powi`](Context::powi) (binary exponentiation), which also accepts a negative + /// base — its sign is fixed by the exponent's parity — so `pow(-x, n)` is in domain here for + /// integer `n`. Both paths are correctly rounded. /// /// # Examples /// @@ -223,7 +438,7 @@ impl Context { if base.is_infinite() || exp.is_infinite() { return Err(FpError::InfiniteInput); } - assert_limited_precision(self.precision); // TODO: we can allow it if exp is integer + assert_limited_precision(self.precision); // shortcuts if exp.is_pos_zero() || exp.is_neg_zero() { @@ -247,21 +462,66 @@ impl Context { FBig::ZERO })); } + if base.is_one() { + // pow(1, y) = 1 for any finite y (exp is finite here — infinities were rejected above). + return Ok(Exact(FBig::ONE)); + } + + // Integer-valued exponent: delegate to the integer-power kernel (binary exponentiation), + // itself correctly rounded via its own Ziv loop. This sidesteps the `exp(y·ln x)` + // amplification entirely, and lets a negative base through — `powi` fixes the sign from + // the exponent's parity. Gated on `is_int` (a cheap exponent check) so the non-integer + // common case skips `to_int`. + if exp.is_int() { + return self.powi(base, exp.to_int().value()); + } + if base.sign() == Sign::Negative { - // TODO: we should allow negative base when exp is an integer + // A non-integer exponent on a negative base has no real value. return Err(FpError::OutOfDomain); } - // x^y = exp(y*ln(x)), use a simple rule for guard bits - let guard_digits = 10 + ceil_usize(self.precision.log2_est()); - let work_context = Context::::new(self.precision + guard_digits); + // x^y = exp(y·ln x), correctly rounded via the Ziv loop. `ln` and `exp` are themselves + // Ziv-correct at the working precision, so the radius comes only from the rounding of the + // `ln`/`mul`/`exp` chain — but `exp` AMPLIFIES the absolute error of its argument `y·ln x` + // by the result magnitude, i.e. by a relative factor of `|y·ln x|`. The radius is + // `result.ulp() · (|y·ln x| + 1) · (B + 8)` where `result.ulp()` is taken at the *working* + // precision, so it shrinks as `B^{-guard}` and the containment test converges. (A radius + // computed at unlimited precision would be constant across retries and never converge for + // a value near a rounding boundary.) The `B + 8` scale covers the `ulp`-vs-`value·B^{1-P}` + // gap plus a safety margin for the chained roundings. + // + // The overflow case is hoisted out of the Ziv closure (which can't return `Err`): if + // `exp(y·ln x)` falls outside the finite exponent range, short-circuit before the loop. + let probe = Context::::new(self.precision + 32); + let ln_x_probe = probe.ln(base, reborrow_cache(&mut cache))?.value(); + let arg_probe = probe.mul(ln_x_probe.repr(), exp)?.value(); + if exp_overflows::(&probe, arg_probe.repr(), &mut cache) { + return Err(if arg_probe.sign() == Sign::Positive { + FpError::Overflow(Sign::Positive) + } else { + FpError::Underflow(Sign::Positive) + }); + } - // ln and exp each consult/extend the shared cache; reborrows are sequential. - let ln_val = work_context.unwrap_fp(work_context.ln(base, reborrow_cache(&mut cache))); - let mul_val = work_context.unwrap_fp(work_context.mul(ln_val.repr(), exp)); - let exp_val = - work_context.unwrap_fp(work_context.exp(mul_val.repr(), reborrow_cache(&mut cache))); - Ok(exp_val.with_precision(self.precision)) + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let ln_x = work.ln(base, reborrow_cache(&mut cache)).unwrap().value(); + let arg = work.mul(ln_x.repr(), exp).unwrap().value(); + let result = work + .exp(arg.repr(), reborrow_cache(&mut cache)) + .unwrap() + .value(); + + // Radius at unlimited precision (exact arithmetic), but built from the *work-precision* + // `result.ulp()` so it carries the `B^{-(p+guard)}` scale and shrinks across retries. + let ulp_w = result.ulp().with_precision(0).value(); + let arg_abs = arg.abs().with_precision(0).value(); + let scale = (B as i32) + 8; + let radius = (ulp_w * (arg_abs + FBig::::ONE)) * scale; + (result, radius) + })) } /// Calculate the exponential function (`eˣ`) on the floating point number under this context. @@ -362,106 +622,40 @@ impl Context { assert_limited_precision(self.precision); - // A simple algorithm: - // - let r = (x - s logB) / Bⁿ, where s = floor(x / logB), such that r < B⁻ⁿ. - // - if the target precision is p digits, then there're only about p/m terms in Tyler series - // - finally, exp(x) = Bˢ * exp(r)^(Bⁿ) - // - the optimal n is √p as given by MPFR - - // Maclaurin series: exp(r) = 1 + Σ(rⁱ/i!) - // There will be about p/log_B(r) summations when calculating the series, to prevent - // loss of significance, we need about log_B(p) guard digits. - let series_guard_digits = ceil_usize(self.precision.log2_est() / B.log2_est()) + 2; - - // Reduction power: the series value is later raised to Bⁿ, which amplifies its - // relative error by a factor of Bⁿ. So the series (and the squarings) must carry - // about n extra base-B digits for the result to come out correct to p digits. We - // use 2n for safety — this mirrors MPFR's working precision q = precy + 2·K + … - // (K ≈ √precy is MPFR's squaring count, the analogue of our n). The log_B(p) - // summation/squaring rounding terms are already covered by series_guard_digits. - let n = 1usize << (self.precision.bit_len() / 2); - let pow_guard_digits = 2 * n; - let work_precision; - - // When minus_one is true and |x| < 1/B, the input is fed into the Maclaurin series without scaling - let no_scaling = minus_one && x.log2_est() < -B.log2_est(); - let (s, n, r) = if no_scaling { - // if minus_one is true and x is already small (x < 1/B), - // then directly evaluate the Maclaurin series without scaling - if x.sign() == Sign::Negative { - // extra digits are required to prevent cancellation during the summation - work_precision = self.precision + 2 * series_guard_digits; - } else { - work_precision = self.precision + series_guard_digits; - } - let context = Context::::new(work_precision); - (0, 0, FBig::new(context.repr_round_ref(x).value(), context)) - } else { - work_precision = self.precision + series_guard_digits + pow_guard_digits; - let context = Context::::new(work_precision); - let x = FBig::new(context.repr_round_ref(x).value(), context); - let logb = context.ln_base::(reborrow_cache(&mut cache)); - let (s, r) = x.div_rem_euclid(logb); - - let s: isize = match s.try_into() { - Ok(v) => v, - Err(_) => { - // |floor(x / ln B)| overflows isize — x is astronomically large, so the - // result is an infinity (x → +∞) or underflows to the limit (x → −∞). - // - // TODO: this branch discards the rounding mode. For a huge - // *negative* x the true exp(x) is positive but below the exponent range, so - // directed Up/Away should return the minimum positive value (and exp_m1 the - // value just above -1) rather than +0 / Exact(-1). - return if input_sign == Sign::Positive { - Err(FpError::Overflow(Sign::Positive)) - } else if minus_one { - Ok(Exact(-FBig::ONE)) // exp_m1(−∞) = −1 (finite) - } else { - Err(FpError::Underflow(Sign::Positive)) // exp(−∞) = +0 - }; - } - }; - (s, n, r) - }; - - let r = r >> n as isize; - let mut factorial = IBig::ONE; - let mut pow = r.clone(); - let mut sum = if no_scaling { - r.clone() - } else { - FBig::ONE + &r - }; - - let mut k = 2; - loop { - factorial *= k; - pow *= &r; - - let increase = &pow / &factorial; - if increase.abs_cmp(&sum.sub_ulp()).is_le() { - break; + // Hoisted overflow check: the reduction quotient s = floor(x/ln B) overflows isize only + // for astronomically large |x| (|x| ≳ 2^61). The Ziv closure below can't return Err, so + // detect that case here and short-circuit to overflow/underflow (matching IEEE limits). + if x.log2_est().abs() > 61.0 { + let probe = Context::::new(self.precision + 64); + let logb = probe.ln_base::(reborrow_cache(&mut cache)); + let x_probe = FBig::new(probe.repr_round_ref(x).value(), probe); + let s_probe = x_probe.div_rem_euclid(logb).0; + if >::try_from(s_probe).is_err() { + return if input_sign == Sign::Positive { + Err(FpError::Overflow(Sign::Positive)) + } else if minus_one { + Ok(Exact(-FBig::ONE)) // exp_m1(−∞) = −1 (finite) + } else { + Err(FpError::Underflow(Sign::Positive)) // exp(−∞) = +0 + }; } - sum += increase; - k += 1; } - if no_scaling { - Ok(sum.with_precision(self.precision)) - } else if minus_one { - // Power at the series' working precision (it already carries the 2n guard - // digits that the Bⁿ powering amplifies away). The final "−1" can cancel at - // most ~1 leading digit here (the |x| < 1/B case is handled by no_scaling), - // which the same guard digits comfortably absorb. - let pow_ctx = Context::::new(work_precision); - let v = pow_ctx.unwrap_fp(pow_ctx.powi(sum.repr(), Repr::::BASE.pow(n).into())); - Ok(((v << s) - FBig::ONE).with_precision(self.precision)) - } else { - let pow_ctx = Context::::new(work_precision); - let v = pow_ctx.unwrap_fp(pow_ctx.powi(sum.repr(), Repr::::BASE.pow(n).into())); - Ok((v << s).with_precision(self.precision)) - } + // Correct rounding via the Ziv loop. Guards: log_B(p) for the series summation/squaring + // rounding, plus `n` for the Bⁿ powering amplification — halved from the pre-Ziv `2n`, + // since Ziv (not the guard count) now certifies correctness. `n ≈ √p` is derived from the + // target precision and is constant across retries. + let series_guard = self.base_guard_digits::(); + let n = 1usize << (self.precision.bit_len() / 2); + Ok(self.ziv(series_guard + n, |guard| { + self.exp_compute::( + x, + self.precision + guard, + minus_one, + n, + reborrow_cache(&mut cache), + ) + })) } } @@ -493,27 +687,6 @@ mod tests { // isize tops out at ~2.1e9 — below any OOM-inducing gap — so the overflow branch // always intervenes first. The fix itself (log2_bounds in round_fract) is // arch-independent; only this dedicated sharp test is 64-bit-only. - #[cfg(target_pointer_width = "64")] - #[test] - fn test_exp_m1_large_negative_no_oom() { - // Regression test: exp_m1 of a large-magnitude negative input subtracts - // 1 from an astronomically small exp(x), so the aligned subtraction hands - // round_fract a sparse sticky tail whose `precision` equals the exponent gap - // (~6.6e18 here). The debug assert in round_fract used to materialize B^precision - // (2^6.6e18 bits → certain OOM); it now uses log2 bounds and must not allocate. - // - // x = -2^62: the gap (~6.6e18) dwarfs memory, yet floor(x/ln2) ≈ -6.6e18 still - // fits 64-bit isize (≈9.2e18), so this does NOT enter the overflow branch. - let ctx = Context::::new(2); - let x = Repr::new(-(IBig::from(1) << 62), 0); - let r = ctx.exp_m1::<2>(&x, None).unwrap().value(); - // exp_m1(-2^62) = -1 + ε (ε > 0 vanishingly small); under Up it rounds above -1. - assert!( - r > -FBig::::ONE, - "exp_m1(-2^62) under Up should round above -1, got {r:?}" - ); - } - #[test] fn test_exact_results_on_unlimited_precision() { // Regression test: values carrying precision 0 (unlimited) — produced by @@ -565,4 +738,29 @@ mod tests { assert!(r.repr().is_neg_zero()); let _ = DBig::ZERO; } + + #[test] + fn test_powf_integer_exponent() { + use crate::DBig; + let ctx = Context::::new(53); + // integer-valued float exponent delegates to powi and supports a negative base (its sign + // is fixed by the exponent's parity): (-5)^3 = -125. + let neg_base = &Repr::<2>::new((-5).into(), 0); + let exp3 = &Repr::<2>::new(3.into(), 0); + let via_powf = ctx.powf::<2>(neg_base, exp3, None).unwrap().value(); + let via_powi = ctx.powi::<2>(neg_base, 3.into()).unwrap().value(); + assert_eq!(via_powf.repr(), via_powi.repr()); + assert_eq!(via_powf.repr().sign(), Sign::Negative); + + // a non-integer exponent on a negative base is out of domain (no real value) + let exp_half = &Repr::<2>::new(5.into(), -1); // 2.5 + assert_eq!(ctx.powf::<2>(neg_base, exp_half, None), Err(FpError::OutOfDomain)); + + // positive base, integer exponent: also routes through powi + let pos_base = &Repr::<2>::new(3.into(), 0); + let exp4 = &Repr::<2>::new(4.into(), 0); + let r = ctx.powf::<2>(pos_base, exp4, None).unwrap().value(); + assert_eq!(r.repr(), ctx.powi::<2>(pos_base, 4.into()).unwrap().value().repr()); + let _ = DBig::ZERO; + } } diff --git a/float/src/fbig_cached_ops.rs b/float/src/fbig_cached_ops.rs index 3ac82869..34aab1d4 100644 --- a/float/src/fbig_cached_ops.rs +++ b/float/src/fbig_cached_ops.rs @@ -12,7 +12,7 @@ use dashu_base::{Abs, CubicRoot, DivEuclid, DivRemEuclid, Inverse, RemEuclid, Si use crate::fbig::FBig; use crate::fbig_cached::CachedFBig; use crate::repr::{Context, Word}; -use crate::round::Round; +use crate::round::{ErrorBounds, Round}; // --------------------------------------------------------------------------- // CachedFBig op CachedFBig (preserves LHS cache) @@ -611,8 +611,9 @@ impl<'a, R: Round, const B: Word> Product<&'a CachedFBig> for CachedFBig`, panicking on -/// error and re-attaching the cache handle. Returns a bare `CachedFBig`. +/// Forward a unary transcendental to its [`Context`] method, panicking on error and re-attaching +/// the cache handle. The context method returns an [`FpResult`] whose rounding flag is dropped by +/// [`unwrap_fp`](crate::Context::unwrap_fp); the cached wrapper never surfaces it. macro_rules! forward_to_context { ($name:ident) => { #[doc = concat!("See [`FBig::", stringify!($name), "`].")] @@ -628,23 +629,6 @@ macro_rules! forward_to_context { }; } -/// Forward a unary function to a [`Context`] method returning `FpResult`, panicking on -/// error and discarding the rounding info. Returns a bare `CachedFBig`. -macro_rules! forward_to_context_unwrap { - ($name:ident) => { - #[doc = concat!("See [`FBig::", stringify!($name), "`].")] - #[inline] - pub fn $name(&self) -> CachedFBig { - let mut c = self.cache.borrow_mut(); - let fbig = self - .fbig - .context - .unwrap_fp(self.fbig.context.$name::(&self.fbig.repr, Some(&mut *c))); - CachedFBig::from_fbig(fbig, &self.cache) - } - }; -} - /// Forward a unary function that delegates to the inner [`FBig`] (no cache needed). macro_rules! forward_to_fbig { ($name:ident) => { @@ -664,22 +648,22 @@ macro_rules! forward_to_fbig { } impl CachedFBig { + forward_to_fbig!(sqrt); + forward_to_fbig!(inv); + forward_to_fbig!(sqr); + forward_to_fbig!(cubic); +} + +// Transcendentals that route through the Ziv-backed (or Ziv-dependent) Context methods require +// `R: ErrorBounds` for their correctness guarantee. +impl CachedFBig { + forward_to_fbig!(powi(exp: dashu_int::IBig)); forward_to_context!(ln); forward_to_context!(ln_1p); forward_to_context!(log2); forward_to_context!(exp); forward_to_context!(exp_m1); - forward_to_fbig!(sqrt); - forward_to_fbig!(inv); - - forward_to_context_unwrap!(sin); - forward_to_context_unwrap!(cos); - forward_to_context_unwrap!(tan); - forward_to_context_unwrap!(asin); - forward_to_context_unwrap!(acos); - forward_to_context_unwrap!(atan); - forward_to_context!(sinh); forward_to_context!(cosh); forward_to_context!(tanh); @@ -687,18 +671,12 @@ impl CachedFBig { forward_to_context!(acosh); forward_to_context!(atanh); - forward_to_fbig!(powi(exp: dashu_int::IBig)); - forward_to_fbig!(sqr); - forward_to_fbig!(cubic); - - /// `self^exp` (see [`FBig::powf`]). - pub fn powf(&self, exp: &Self) -> Self { - let context = Context::max(self.fbig.context, exp.fbig.context); - let mut c = self.cache.borrow_mut(); - let fbig = - context.unwrap_fp(context.powf::(&self.fbig.repr, &exp.fbig.repr, Some(&mut *c))); - Self::from_fbig(fbig, &self.cache) - } + forward_to_context!(sin); + forward_to_context!(cos); + forward_to_context!(tan); + forward_to_context!(asin); + forward_to_context!(acos); + forward_to_context!(atan); /// Sine and cosine together (see [`FBig::sin_cos`]). pub fn sin_cos(&self) -> (Self, Self) { @@ -711,17 +689,6 @@ impl CachedFBig { ) } - /// Hyperbolic sine and cosine together (see [`FBig::sinh_cosh`]). - pub fn sinh_cosh(&self) -> (Self, Self) { - let mut guard = self.cache.borrow_mut(); - let cache = Some(&mut *guard); - let (s, c) = self.fbig.context.sinh_cosh::(&self.fbig.repr, cache); - ( - Self::from_fbig(self.fbig.context.unwrap_fp(s), &self.cache), - Self::from_fbig(self.fbig.context.unwrap_fp(c), &self.cache), - ) - } - /// `atan2(y, x)` (see [`FBig::atan2`]). pub fn atan2(&self, x: &Self) -> Self { let mut c = self.cache.borrow_mut(); @@ -732,4 +699,24 @@ impl CachedFBig { )); Self::from_fbig(fbig, &self.cache) } + + /// `self^exp` (see [`FBig::powf`]). + pub fn powf(&self, exp: &Self) -> Self { + let context = Context::max(self.fbig.context, exp.fbig.context); + let mut c = self.cache.borrow_mut(); + let fbig = + context.unwrap_fp(context.powf::(&self.fbig.repr, &exp.fbig.repr, Some(&mut *c))); + Self::from_fbig(fbig, &self.cache) + } + + /// Hyperbolic sine and cosine together (see [`FBig::sinh_cosh`]). + pub fn sinh_cosh(&self) -> (Self, Self) { + let mut guard = self.cache.borrow_mut(); + let cache = Some(&mut *guard); + let (s, c) = self.fbig.context.sinh_cosh::(&self.fbig.repr, cache); + ( + Self::from_fbig(self.fbig.context.unwrap_fp(s), &self.cache), + Self::from_fbig(self.fbig.context.unwrap_fp(c), &self.cache), + ) + } } diff --git a/float/src/iter.rs b/float/src/iter.rs index 9e5701ad..561b8e70 100644 --- a/float/src/iter.rs +++ b/float/src/iter.rs @@ -5,38 +5,9 @@ use crate::{ fbig::FBig, repr::{Context, Repr, Word}, round::Round, - utils::{shl_digits, shl_digits_in_place}, }; use core::iter::{Product, Sum}; -/// Exact (unrounded) sum of two finite [`Repr`] values. -/// -/// The result exponent is `min(lhs.exponent, rhs.exponent)`. The larger-exponent operand's -/// significand is shifted up by the exponent gap (in base-`B` digits), then the signed significands -/// are added — `IBig` addition handles opposite signs as an exact subtraction. Zero operands short -/// circuit so that a `-0` sentinel exponent (`-1`) can't bleed into the result exponent. The result -/// is rebuilt with [`Repr::new`], which strips trailing zeros. -fn exact_add(mut lhs: Repr, rhs: &Repr) -> Repr { - if lhs.significand.is_zero() { - return rhs.clone(); - } - if rhs.significand.is_zero() { - return lhs; - } - - if lhs.exponent >= rhs.exponent { - let ediff = (lhs.exponent - rhs.exponent) as usize; - shl_digits_in_place::(&mut lhs.significand, ediff); - lhs.significand += &rhs.significand; - lhs.exponent = rhs.exponent; - } else { - let ediff = (rhs.exponent - lhs.exponent) as usize; - let rhs_sig = shl_digits::(&rhs.significand, ediff); - lhs.significand += rhs_sig; // lhs.exponent is already the minimum - } - Repr::new(lhs.significand, lhs.exponent) -} - /// Correctly-rounded summation of finite floats. /// /// Every addend is accumulated exactly at the [`Repr`] level (no per-step rounding); the exact total @@ -58,7 +29,7 @@ fn precise_sum( }; for (repr, ctx) in iter { assert_finite(&repr); - acc = exact_add(acc, &repr); + acc = acc + &repr; context = Context::max(context, ctx); } // Exact cancellation can leave `acc` with a zero significand whose exponent coincides with the @@ -223,17 +194,17 @@ mod tests { let _: F = core::iter::once(F::INFINITY).sum(); } - // Exercises the zero short-circuit in `exact_add`: a `-0` sentinel exponent must not corrupt the - // accumulator exponent. + // Exercises the zero short-circuit in `impl Add for &Repr`: a `-0` sentinel exponent must not + // corrupt the accumulator exponent. #[test] - fn exact_add_neg_zero_is_identity() { + fn repr_add_neg_zero_is_identity() { let nz = Repr::<10>::neg_zero(); let x = r::<10>(5, -2); // 0.05 // -0 + x == x, x + -0 == x (exact, pre-rounding) - assert_eq!(exact_add(nz.clone(), &x), x); - assert_eq!(exact_add(x.clone(), &nz), x); - assert_eq!(exact_add(nz.clone(), &nz), Repr::<10>::neg_zero()); + assert_eq!(&nz + &x, x); + assert_eq!(&x + &nz, x); + assert_eq!(&nz + &nz, Repr::<10>::neg_zero()); // 1 + -1 cancels to +0 (not -0) - assert_eq!(exact_add(r::<10>(1, 0), &r::<10>(-1, 0)), Repr::<10>::zero()); + assert_eq!(&r::<10>(1, 0) + &r::<10>(-1, 0), Repr::<10>::zero()); } } diff --git a/float/src/lib.rs b/float/src/lib.rs index b840457d..6f7091a9 100644 --- a/float/src/lib.rs +++ b/float/src/lib.rs @@ -87,6 +87,7 @@ mod mul; pub mod ops; mod parse; mod repr; +mod repr_ops; mod root; pub mod round; mod round_ops; @@ -94,6 +95,7 @@ mod shift; mod sign; mod third_party; mod utils; +mod ziv; // All the public items from third_party will be exposed #[allow(unused_imports)] diff --git a/float/src/log.rs b/float/src/log.rs index 325dbddf..78b15192 100644 --- a/float/src/log.rs +++ b/float/src/log.rs @@ -2,7 +2,7 @@ use dashu_base::{ utils::{next_down, next_up}, AbsOrd, Approximation::*, - EstimatedLog2, Sign, + EstimatedLog2, PowerOfTwo, Sign, UnsignedAbs, }; use dashu_int::IBig; @@ -10,9 +10,9 @@ use crate::{ error::{assert_finite, assert_limited_precision, FpError, FpResult}, fbig::FBig, math::cache::{reborrow_cache, ConstCache}, + math::trig::series_radius, repr::{Context, Repr, Word}, - round::{Round, Rounded}, - utils::ceil_usize, + round::{mode, ErrorBounds, Round, Rounded}, }; use core::cmp::Ordering; @@ -63,7 +63,7 @@ impl EstimatedLog2 for FBig { } } -impl FBig { +impl FBig { /// Calculate the natural logarithm function (`log(x)`) on the float number. /// /// # Examples @@ -98,12 +98,10 @@ impl FBig { self.context.unwrap_fp(self.context.ln_1p(&self.repr, None)) } - /// Calculate the base-2 logarithm function (`log2(x)`) on the float number. + /// Calculate the base-2 logarithm (`log2(x)`) on the float number. /// - /// This is evaluated as `ln(x) / ln(2)` at an elevated working precision and rounded - /// once, so it is correctly rounded to this number's precision (unlike - /// [`log2_bounds`][dashu_base::EstimatedLog2::log2_bounds], which only gives an - /// f32-precision magnitude estimate). + /// Correctly rounded to the context's precision under any rounding mode. For an exact power + /// of two the result is the exact integer `log2(x)`. /// /// # Examples /// @@ -111,7 +109,7 @@ impl FBig { /// # use core::str::FromStr; /// # use dashu_base::ParseError; /// # use dashu_float::DBig; - /// let a = DBig::from_str("8.0")?; + /// let a = DBig::from_str("8")?; /// assert_eq!(a.log2(), DBig::from_str("3")?); /// # Ok::<(), ParseError>(()) /// ``` @@ -121,6 +119,12 @@ impl FBig { } } +// `ln2`/`ln10`/`iacoth`/`ln_base`/`ln_compute` are the near-correct logarithm primitives: they +// evaluate the series at a working precision and round once, without a Ziv certification step. +// They live on `R: Round` so that base conversion (`with_base_and_precision`, which only needs a +// near-correct constant `ln(B)`) can use them without inheriting the `ErrorBounds` bound. The +// correctly-rounded public `ln`/`ln_1p` (in the `ErrorBounds` impl below) wrap `ln_compute` in a +// Ziv loop. impl Context { /// Calculate log(2) /// @@ -160,7 +164,19 @@ impl Context { 2 => self.ln2(None), 10 => self.ln10(None), i if i.is_power_of_two() => self.ln2(None) * i.trailing_zeros(), - _ => self.unwrap_fp(self.ln(&Repr::new(Repr::::BASE.into(), 0), None)), + _ => { + // Near-correct ln(B) via the atanh series (no Ziv certification — base conversion + // only needs a near-correct constant). `ln_compute` is on `R: Round`, so this keeps + // `ln_base` callable from `R: Round` contexts (base conversion). + let guard = self.base_guard_digits::() + 2; + self.ln_compute::( + &Repr::new(Repr::::BASE.into(), 0), + self.precision + guard, + false, + None, + ) + .0 + } } } @@ -193,7 +209,7 @@ impl Context { // L(n) = (Q + T) / (n·Q). Extra guard digits absorb the division's rounding // (the binary-splitting state is exact, so only this single round loses anything). - let guard_digits = ceil_usize(self.precision.log2_est() / B.log2_est()); + let guard_digits = self.base_guard_digits::(); let work_context = Self::new(self.precision + guard_digits + 2); let num = work_context.convert_int::(q.as_ibig() + &t).value(); @@ -201,6 +217,107 @@ impl Context { num / denom } + /// Evaluate `ln(x)` (or `ln(x+1)` when `one_plus`) at `work_precision` via the atanh series, + /// returning `(value, error_radius)`. + /// + /// This is the near-correct computation core shared by the public Ziv-backed `ln`/`ln_1p` + /// (which wrap it in a retry loop) and by `ln_base` (which only needs a near-correct constant + /// `ln(B)`). It lives on `R: Round` so those near-correct callers don't inherit the + /// `ErrorBounds` bound. The radius is a provable upper bound on `|value − true|`, derived from + /// the term count (every series step is correctly rounded; the truncated tail is `< 1 ulp` by + /// the break test). + pub(crate) fn ln_compute( + &self, + x: &Repr, + mut work_precision: usize, + one_plus: bool, + mut cache: Option<&mut ConstCache>, + ) -> (FBig, FBig) { + // log(x) = log(x·B⁻ˢ) + s·log(B), with s = floor(log_B(x)) so x·B⁻ˢ ∈ [1, B). + let context = Context::::new(work_precision); + let x = FBig::new(context.repr_round_ref(x).value(), context); + + // When one_plus is true and |x| < 1/B, the input is fed into the Maclaurin without scaling + let no_scaling = one_plus && x.log2_est() < -B.log2_est(); + + let (s, mut x_scaled) = if no_scaling { + (0, x) + } else { + let x = if one_plus { x + FBig::ONE } else { x }; + + let log2 = x.log2_bounds().0; + let s = log2 as isize - (log2 < 0.) as isize; // floor(log2(x)) + + let x_scaled = if B == 2 { + x >> s + } else if s > 0 { + x / (IBig::ONE << s as usize) + } else { + x * (IBig::ONE << (-s) as usize) + }; + debug_assert!(x_scaled >= FBig::::ONE); + (s, x_scaled) + }; + + if s < 0 || x_scaled.repr.sign() == Sign::Negative { + // when s or x_scaled is negative, the final addition is actually a subtraction, + // therefore we need to double the precision to get the correct result + work_precision += self.precision; + x_scaled.context.precision = work_precision; + } + let work_context = Context::new(work_precision); + + // after the number is scaled to nearly one, use Maclaurin series on log(x) = 2atanh(z): + // let z = (x-1)/(x+1) < 1, log(x) = 2atanh(z) = 2Σ(z²ⁱ⁺¹/(2i+1)) for i = 1,3,5,... + let z = if no_scaling { + let d = &x_scaled + (FBig::ONE + FBig::ONE); + x_scaled / d + } else { + (&x_scaled - FBig::ONE) / (x_scaled + FBig::ONE) + }; + let z2 = z.sqr(); + let mut pow = z.clone(); + let mut sum = z; + let mut terms: usize = 1; // the leading z term + + let mut k: usize = 3; + loop { + pow *= &z2; + + let increase = &pow / work_context.convert_int::(k.into()).value(); + if increase.abs_cmp(&sum.sub_ulp()).is_le() { + break; + } + + sum += increase; + k += 2; + terms += 1; + } + + // compose the logarithm of the original number + let result: FBig = if no_scaling { + 2 * sum.clone() + } else { + 2 * sum.clone() + (s * work_context.ln2::(reborrow_cache(&mut cache))) + }; + + // Provable error radius, expressed in `result`-ULPs (not `sum`-ULPs). Each series step + // rounds once (< 1 ULP of the running sum) and the truncated tail is < 1 ULP by the break + // test, so |sum − true| < (terms + 2)·ulp(sum); result = 2·sum + s·ln2 amplifies by ~2 and + // adds a few reconstruction ULPs. Since result ≈ 2·sum, ulp(result) ≈ 2·ulp(sum), giving + // |result − true| < (terms + 2)·ulp(result) + overhead — we carry a generous margin. + // + // Basing the radius on `result.ulp()` (not `sum.ulp()`) keeps its exponent aligned with + // `a` (= result) in the Ziv containment test, so `a − e` avoids a slow exponent-misaligned + // unlimited-precision subtract — a ~3× speedup on `ln` at high precision. + let radius = series_radius(&result, terms); + (result, radius) + } +} + +// `ln`/`ln_1p` are correctly rounded via the Ziv loop, whose containment test needs the rounding +// preimage (`R: ErrorBounds`). They delegate the series to `ln_compute`. +impl Context { /// Calculate the natural logarithm function (`log(x)`) on the float number under this context. /// /// # Examples @@ -236,48 +353,6 @@ impl Context { Ok(self.ln_internal(x, false, cache)) } - /// Calculate the base-2 logarithm function (`log2(x)`) on the float number under this - /// context. - /// - /// This is evaluated as `ln(x) / ln(2)` at an elevated working precision and rounded - /// once, so it is correctly rounded to this context's precision (unlike - /// [`log2_bounds`][dashu_base::EstimatedLog2::log2_bounds], which only gives an - /// f32-precision magnitude estimate). - /// - /// # Examples - /// - /// ``` - /// # use core::str::FromStr; - /// # use dashu_base::ParseError; - /// # use dashu_float::DBig; - /// use dashu_base::Approximation::*; - /// use dashu_float::{Context, round::mode::HalfEven}; - /// - /// let context = Context::::new(8); - /// let a = DBig::from_str("10")?; - /// // log2(10) ≈ 3.3219281 - /// assert_eq!(context.log2(&a.repr(), None).unwrap().value(), DBig::from_str("3.3219281")?); - /// # Ok::<(), ParseError>(()) - /// ``` - #[inline] - pub fn log2( - &self, - x: &Repr, - cache: Option<&mut ConstCache>, - ) -> FpResult> { - if x.is_infinite() { - return Err(FpError::InfiniteInput); - } - if x.significand.is_zero() { - // log2(±0) = -inf (a value, not an error) - return Ok(Exact(FBig::new(Repr::neg_infinity(), *self))); - } - if x.sign() == Sign::Negative { - return Err(FpError::OutOfDomain); - } - Ok(self.log2_internal(x, cache)) - } - /// Calculate the natural logarithm function (`log(x+1)`) on the float number under this context. /// /// # Examples @@ -340,103 +415,147 @@ impl Context { assert_limited_precision(self.precision); - // A simple algorithm: - // - let log(x) = log(x/2^s) + slog2 where s = floor(log2(x)) - // - such that x*2^s is close to but larger than 1 (and x*2^s < 2) - let guard_digits = ceil_usize(self.precision.log2_est() / B.log2_est()) + 2; - let mut work_precision = self.precision + guard_digits + one_plus as usize; - let context = Context::::new(work_precision); - let x = FBig::new(context.repr_round_ref(x).value(), context); - - // When one_plus is true and |x| < 1/B, the input is fed into the Maclaurin without scaling - let no_scaling = one_plus && x.log2_est() < -B.log2_est(); - - let (s, mut x_scaled) = if no_scaling { - (0, x) - } else { - let x = if one_plus { x + FBig::ONE } else { x }; - - let log2 = x.log2_bounds().0; - let s = log2 as isize - (log2 < 0.) as isize; // floor(log2(x)) - - let x_scaled = if B == 2 { - x >> s - } else if s > 0 { - x / (IBig::ONE << s as usize) - } else { - x * (IBig::ONE << (-s) as usize) - }; - debug_assert!(x_scaled >= FBig::::ONE); - (s, x_scaled) - }; - - if s < 0 || x_scaled.repr.sign() == Sign::Negative { - // when s or x_scaled is negative, the final addition is actually a subtraction, - // therefore we need to double the precision to get the correct result - work_precision += self.precision; - x_scaled.context.precision = work_precision; - }; - let work_context = Context::new(work_precision); - - // after the number is scaled to nearly one, use Maclaurin series on log(x) = 2atanh(z): - // let z = (x-1)/(x+1) < 1, log(x) = 2atanh(z) = 2Σ(z²ⁱ⁺¹/(2i+1)) for i = 1,3,5,... - // similar to iacoth, the required iterations stop at i = -p/2log_B(z), and we need log_B(i) guard bits - let z = if no_scaling { - let d = &x_scaled + (FBig::ONE + FBig::ONE); - x_scaled / d - } else { - (&x_scaled - FBig::ONE) / (x_scaled + FBig::ONE) - }; - let z2 = z.sqr(); - let mut pow = z.clone(); - let mut sum = z; - - let mut k: usize = 3; - loop { - pow *= &z2; - - let increase = &pow / work_context.convert_int::(k.into()).value(); - if increase.abs_cmp(&sum.sub_ulp()).is_le() { - break; - } + // Correct rounding via the Ziv loop: `ln_compute` evaluates the atanh series at `p + guard` + // and reports a provable error radius; the driver retries with more guard digits until the + // approximation's error interval lies entirely inside one rounding bin. The guard is a + // *performance* knob (first-attempt hit rate), not a correctness backstop — Ziv certifies + // the result. (The pre-Ziv `+ 2` is retained: with the conservative radius below it is still + // needed for the first attempt to clear the half-ulp preimage at typical precisions.) + let base_guard = self.base_guard_digits::() + 2; + self.ziv(base_guard + one_plus as usize, |guard| { + self.ln_compute::(x, self.precision + guard, one_plus, reborrow_cache(&mut cache)) + }) + } - sum += increase; - k += 2; + /// Calculate the base-2 logarithm (`log2(x)`) on the float number under this context. + /// + /// Correctly rounded to the context's precision under any rounding mode; for an exact power + /// of two the result is the exact integer `log2(x)`. + /// + /// # Domain + /// + /// `log2(±0) = −∞` and a negative (non-zero) input is out of domain; an infinite input is an + /// error (a finite context cannot produce the infinite `log2(+∞) = +∞` exactly). + /// + /// # Examples + /// + /// ``` + /// # use core::str::FromStr; + /// # use dashu_base::ParseError; + /// # use dashu_float::DBig; + /// use dashu_base::Approximation::*; + /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}}; + /// + /// let context = Context::::new(4); + /// let a = DBig::from_str("10")?; + /// assert_eq!(context.log2(&a.repr(), None), Ok(Inexact(DBig::from_str("3.322")?, AddOne))); + /// # Ok::<(), ParseError>(()) + /// ``` + #[inline] + pub fn log2( + &self, + x: &Repr, + cache: Option<&mut ConstCache>, + ) -> FpResult> { + if x.is_infinite() { + return Err(FpError::InfiniteInput); } - - // compose the logarithm of the original number - let result: FBig = if no_scaling { - 2 * sum - } else { - 2 * sum + s * work_context.ln2::(reborrow_cache(&mut cache)) - }; - result.with_precision(self.precision) + if x.significand.is_zero() { + // log2(±0) = -inf (a value, not an error) + return Ok(Exact(FBig::new(Repr::neg_infinity(), *self))); + } + if x.sign() == Sign::Negative { + return Err(FpError::OutOfDomain); + } + Ok(self.log2_internal(x, cache)) } - /// `log2(x) = ln(x) / ln(2)`, correctly rounded to this context's precision. - /// - /// Both `ln(x)` and `ln(2)` are evaluated at a working precision above the target and - /// the quotient is rounded once. The result magnitude `|log2(x)| ≈ |floor(log2 x)|` - /// grows with the operand, which tracks the division's error amplification, so the - /// per-ulp error is ~`B^(target - work)`: a few guard digits certify the final round - /// across the whole magnitude range (no manual scaling needed). fn log2_internal( &self, x: &Repr, mut cache: Option<&mut ConstCache>, ) -> Rounded> { assert_finite(x); + + // Exact shortcuts first — they also cover unlimited precision, which the Ziv loop below + // rejects via its limited-precision assertion. if x.is_one() { - return Exact(FBig::ZERO); // log2(1) = +0 (exact; also covers unlimited precision) + return Exact(FBig::ZERO); // log2(1) = +0 } - let guard = ceil_usize(self.precision.log2_est() / B.log2_est()) + 3; - let work = Context::::new(self.precision + guard); - let ln_x = work - .ln_internal(x, false, reborrow_cache(&mut cache)) - .value(); - let ln2 = work.ln2::(reborrow_cache(&mut cache)); - (ln_x / ln2).with_precision(self.precision) + // Exact power-of-two shortcut: if x = 2^k for an integer k, log2(x) = k. This is *required* + // for directed rounding — the Ziv loop below cannot certify an exactly-representable + // result whose true value sits on a rounding boundary (its shrinking error interval + // always straddles the boundary), so without this shortcut log2(2^-159) under `Up` would + // exhaust the retry cap and return k + 1 ulp instead of the exact k. + // + // log2(x) = log2(significand) + exponent·log2(B). With significand = 2^m this is an exact + // integer whenever log2(B) is integral (B a power of two), or — for a non-power-of-two + // base — when the exponent is zero. + let mag = (&x.significand).unsigned_abs(); + if mag.is_power_of_two() && (x.exponent == 0 || B.is_power_of_two()) { + let m = mag.trailing_zeros().unwrap(); // = log2(significand) + let log2_b = B.trailing_zeros() as isize; + let k = IBig::from(m) + IBig::from(x.exponent) * IBig::from(log2_b); + return self.convert_int::(k); + } + + assert_limited_precision(self.precision); + + // log2(x) = ln(x)/ln(2), correctly rounded via the Ziv loop. Rounding ln(x) and ln(2) + // separately and dividing once is only *near*-correct: under directed rounding, rounding + // both operands toward the mode does not bound the quotient (enlarging a positive + // denominator shrinks it). Instead each `ln_compute` reports a provable error radius, and + // the two radii are carried through the division as an outward-rounded interval [lo, hi] + // that is guaranteed to contain the true log2(x); the driver certifies once that interval + // lies inside a single rounding bin. + let initial_guard = self.base_guard_digits::() + 4; + self.ziv(initial_guard, |guard| { + let work_precision = self.precision + guard; + let (lx, ex) = + self.ln_compute::(x, work_precision, false, reborrow_cache(&mut cache)); + // ln(2) via the same near-correct primitive so it carries a provable radius too. + let two = Repr::new(IBig::from(2), 0); + let (l2, e2) = + self.ln_compute::(&two, work_precision, false, reborrow_cache(&mut cache)); + + // True ln(x) ∈ [lx−ex, lx+ex] and true ln(2) ∈ [l2−e2, l2+e2] ⊆ (0, ∞). With a + // positive denominator the quotient ln(x)/ln(2) is minimized by the low numerator + // over the high denominator and maximized by the converse. Directing each endpoint's + // rounding outward (lo down, hi up) keeps [lo, hi] a true containing interval. + let down = Context::::new(work_precision); + let up = Context::::new(work_precision); + let nx_lo = down.sub(&lx.repr, &ex.repr).unwrap().value(); + let nx_hi = up.add(&lx.repr, &ex.repr).unwrap().value(); + let d_lo = down.sub(&l2.repr, &e2.repr).unwrap().value(); + let d_hi = up.add(&l2.repr, &e2.repr).unwrap().value(); + debug_assert!( + d_lo.repr.sign() == Sign::Positive, + "ln(2) lower bound must stay positive (guard digits keep e2 ≪ ln 2 ≈ 0.693)" + ); + let lo = down.div(&nx_lo.repr, &d_hi.repr).unwrap().value(); + let hi = up.div(&nx_hi.repr, &d_lo.repr).unwrap().value(); + + // Working-precision estimate; the driver re-rounds it to the target precision, so the + // mode used here is immaterial to correctness. + let value = Context::::new(work_precision) + .div(&lx.repr, &l2.repr) + .unwrap() + .value(); + + // Radius: a provable bound on |value − true|. The true value lies in [lo, hi], and + // `value` is within one working ulp of lx/l2 ∈ [lo, hi], so |value − true| ≤ + // (hi − lo) + ulp_w. Computed at unlimited precision so the bound arithmetic is exact + // (no rounding that could under-report it), yet scaled by the working-precision span + // and ulp so it shrinks as the guard grows and the loop converges. `lo`/`hi` were + // rounded under Down/Up; their *values* are mode-independent, so rebuild them in the + // target mode R via their reprs to keep the arithmetic single-mode. + let unlim = Context::::new(0); + let span = FBig::new(hi.repr.clone(), unlim) - FBig::new(lo.repr.clone(), unlim); + let ulp_w = value.ulp().with_precision(0).value(); + let radius = span + ulp_w; + (value, radius) + }) } } @@ -453,85 +572,6 @@ mod tests { assert_eq!(r.repr().sign(), Sign::Negative); } - #[test] - fn test_log2_domain() { - let ctx = Context::::new(53); - // log2(±0) = -inf (a value, not an error) - let r = ctx.log2::<2>(&Repr::<2>::zero(), None).unwrap().value(); - assert!(r.repr().is_infinite()); - assert_eq!(r.repr().sign(), Sign::Negative); - // log2(negative) is out of domain - assert_eq!(ctx.log2::<2>(&Repr::new(IBig::from(-1), 0), None), Err(FpError::OutOfDomain)); - } - - #[test] - fn test_log2_powers_of_two() { - // log2(2^k) = k, including log2(1) = 0. The only magnitude estimate available - // (log2_bounds) is f32-precision, so deriving directed log2 from it gives a small - // *negative* value for log2(1) and is loose by many ULPs elsewhere; the - // correctly-rounded log2 returns the exact integer value for every power of two - // (log2(1) is tagged Exact via the shortcut; other powers are value-exact). - let ctx = Context::::new(53); - for k in [0usize, 1, 2, 3, 10, 50, 100, 200, 1000] { - let x = Repr::new(IBig::from(1) << k, 0); // 2^k - let r = ctx.log2::<2>(&x, None).unwrap(); - let expected = ctx.convert_int::<2>(IBig::from(k)).value(); - assert_eq!(r.value(), expected, "log2(2^{k}) = {k}"); - } - // log2(1) = 0 is certified Exact via the special-case shortcut. - assert!(matches!(ctx.log2::<2>(&Repr::new(IBig::from(1), 0), None).unwrap(), Exact(_))); - } - - #[test] - fn test_log2_correctly_rounded() { - // log2(x) must equal the correctly-rounded ln(x)/ln(2) oracle (computed at much - // higher precision then rounded down) across magnitude ranges — tiny, normal, and - // huge inputs, in both base 2 and base 10. - let p = 30; - let ctx = Context::::new(p); - let hi = Context::::new(200); - - let cases_base2 = [ - Repr::new(IBig::from(3), 0), - Repr::new(IBig::from(10), 0), - Repr::new(IBig::from(1), -24), // near the fuzzer's 2^-24 region - Repr::new(IBig::from(1) << 50, 0), // huge - Repr::new(IBig::from(1), -50), // tiny - ]; - for x in cases_base2 { - let got = ctx.log2::<2>(&x, None).unwrap().value(); - let want = (hi.ln::<2>(&x, None).unwrap().value() / hi.ln2::<2>(None)) - .with_precision(p) - .value(); - assert_eq!(got, want, "base-2 log2 mismatch for x={x:?}"); - } - - // base 10: log2(x) = ln(x)/ln(2) with a decimal significand - let x = Repr::new(IBig::from(123456), 0); - let got = ctx.log2::<10>(&x, None).unwrap().value(); - let want = (hi.ln::<10>(&x, None).unwrap().value() / hi.ln2::<10>(None)) - .with_precision(p) - .value(); - assert_eq!(got, want, "base-10 log2 mismatch"); - } - - #[test] - fn test_log2_directed_bounds() { - // Down rounds toward -inf, Up toward +inf, so the correctly-rounded (HalfEven) - // value must lie within [Down, Up]. - use crate::round::mode::{Down, Up}; - let p = 20; - let mid = Context::::new(p); - for x in [Repr::new(IBig::from(10), 0), Repr::new(IBig::from(3), 13)] { - let down = Context::::new(p).log2::<2>(&x, None).unwrap().value(); - let up = Context::::new(p).log2::<2>(&x, None).unwrap().value(); - let half = mid.log2::<2>(&x, None).unwrap().value(); - assert!(down <= half, "down <= half failed for x={x:?}: {down:?} vs {half:?}"); - assert!(half <= up, "half <= up failed for x={x:?}: {half:?} vs {up:?}"); - assert!(down <= up, "down <= up failed for x={x:?}"); - } - } - #[test] fn test_iacoth() { let context = Context::::new(10); @@ -587,4 +627,126 @@ mod tests { .unwrap() ); } + + #[test] + fn test_log2_domain() { + let ctx = Context::::new(53); + // log2(±0) = -inf (a value, not an error) + let r = ctx.log2::<2>(&Repr::<2>::zero(), None).unwrap().value(); + assert!(r.repr.is_infinite()); + assert_eq!(r.repr.sign(), Sign::Negative); + // log2(negative) is out of domain + assert!(matches!( + ctx.log2::<2>(&Repr::new((-1).into(), 0), None), + Err(FpError::OutOfDomain) + )); + // an infinite input is rejected + assert!(matches!(ctx.log2::<2>(&Repr::infinity(), None), Err(FpError::InfiniteInput))); + } + + #[test] + fn test_log2_exact_power_of_two() { + // log2(2^k) = k exactly under every rounding mode. Regression for the directed-rounding + // defect: rounding ln(x) and ln(2) each toward the mode and dividing once does not bound + // the quotient, so previously log2(2^-159) under `Up` returned -159 + 1 ulp. + let p = 53; + for k in [0isize, 1, -1, 5, 159, -159, 1000, -1000] { + let x = Repr::<2>::new(IBig::from(1), k); // 2^k + let r_down = Context::::new(p) + .log2::<2>(&x, None) + .unwrap() + .value(); + let r_up = Context::::new(p) + .log2::<2>(&x, None) + .unwrap() + .value(); + let r_zero = Context::::new(p) + .log2::<2>(&x, None) + .unwrap() + .value(); + let r_he = Context::::new(p) + .log2::<2>(&x, None) + .unwrap() + .value(); + // Every directed mode produces the identical value — no mode-dependent ulp. + assert_eq!(r_down.repr, r_he.repr, "Down != HalfEven for log2(2^{k})"); + assert_eq!(r_up.repr, r_he.repr, "Up != HalfEven for log2(2^{k})"); + assert_eq!(r_zero.repr, r_he.repr, "Zero != HalfEven for log2(2^{k})"); + // And that value is exactly k. + assert_eq!(r_he.to_int().value(), IBig::from(k), "value for log2(2^{k})"); + } + } + + #[test] + fn test_log2_exact_power_of_two_decimal_base() { + // In a non-power-of-two base the shortcut still fires when the exponent is zero: a + // significand that is itself a power of two makes x = 2^m exactly. + let p = 53; + for (sig, want) in [(8i32, 3isize), (1024, 10), (2, 1), (32, 5)] { + let x = Repr::<10>::new(IBig::from(sig), 0); + let r_down = Context::::new(p) + .log2::<10>(&x, None) + .unwrap() + .value(); + let r_up = Context::::new(p) + .log2::<10>(&x, None) + .unwrap() + .value(); + let r_he = Context::::new(p) + .log2::<10>(&x, None) + .unwrap() + .value(); + assert_eq!(r_down.repr, r_he.repr, "Down != HalfEven for log2({sig}) base 10"); + assert_eq!(r_up.repr, r_he.repr, "Up != HalfEven for log2({sig}) base 10"); + assert_eq!(r_he.to_int().value(), IBig::from(want), "value for log2({sig}) base 10"); + } + } + + /// For a non-power-of-two significand `sig` (so `log2` is irrational and never lands on a + /// rounding boundary), each directed result must equal a high-precision oracle rounded to the + /// target precision under the same mode — the definition of correct rounding. + fn check_log2_directed_matches_oracle(sig: u32, p: usize) { + let oracle_ctx = Context::::new(p + 40); + let x = Repr::::new(IBig::from(sig), 0); + let oracle = oracle_ctx.log2::(&x, None).unwrap().value(); + + let want_down = Context::::new(p) + .repr_round_ref(&oracle.repr) + .value(); + let want_up = Context::::new(p) + .repr_round_ref(&oracle.repr) + .value(); + let want_he = Context::::new(p) + .repr_round_ref(&oracle.repr) + .value(); + + let got_down = Context::::new(p) + .log2::(&x, None) + .unwrap() + .value(); + let got_up = Context::::new(p) + .log2::(&x, None) + .unwrap() + .value(); + let got_he = Context::::new(p) + .log2::(&x, None) + .unwrap() + .value(); + + assert_eq!(got_down.repr, want_down, "log2({sig}) base {B} under Down"); + assert_eq!(got_up.repr, want_up, "log2({sig}) base {B} under Up"); + assert_eq!(got_he.repr, want_he, "log2({sig}) base {B} under HalfEven"); + } + + #[test] + fn test_log2_directed_matches_oracle() { + let p = 24; + for sig in [3u32, 7, 10, 12345, 65537] { + check_log2_directed_matches_oracle::<2>(sig, p); + } + // Exercise a non-power-of-two base through the Ziv interval path too. + for sig in [3u32, 7, 10, 12345] { + check_log2_directed_matches_oracle::<10>(sig, p); + } + } } diff --git a/float/src/math/cache.rs b/float/src/math/cache.rs index 158d00c4..0757a4aa 100644 --- a/float/src/math/cache.rs +++ b/float/src/math/cache.rs @@ -213,14 +213,18 @@ impl ConstCache { .value() } _ => { - // generic base: no cached L(n) sub-series applies, so compute - // ln(B) directly through Context::ln on the base literal. + // generic base: no cached L(n) sub-series applies, so compute ln(B) directly via + // the near-correct `ln_compute` (the cache only stores a near-correct constant — a + // cached value is re-rounded on read, and Ziv isn't needed here). let ctx = Context::::new(precision); - ctx.unwrap_fp(ctx.ln::( + ctx.ln_compute::( &Repr::new(Repr::::BASE.into(), 0), // no cache for the generic base (its L(n) isn't cached) + precision, + false, None, - )) + ) + .0 } } } diff --git a/float/src/math/hyper.rs b/float/src/math/hyper.rs index 85f3a9c7..25dd963f 100644 --- a/float/src/math/hyper.rs +++ b/float/src/math/hyper.rs @@ -14,17 +14,18 @@ use crate::{ error::{assert_limited_precision, FpError}, + exp::exp_overflows, fbig::FBig, math::{ cache::{reborrow_cache, ConstCache}, FpResult, }, repr::{Context, Repr, Word}, - round::Round, + round::{ErrorBounds, Round}, }; use dashu_base::{Abs, AbsOrd, Approximation::Exact, Sign}; -impl Context { +impl Context { /// Hyperbolic sine. pub fn sinh( &self, @@ -39,20 +40,30 @@ impl Context { // sinh(±0) = ±0 return Ok(Exact(FBig::new(signed_zero_repr(x), *self))); } - - // sinh(x) = (exp_m1(x) - exp_m1(-x)) / 2 (cancellation-free) - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let neg_x = -x_f.clone(); - let ep = work.exp_m1(&x_f.repr, reborrow_cache(&mut cache)); - let em = work.exp_m1(&neg_x.repr, reborrow_cache(&mut cache)); - match (ep, em) { - (Ok(ep), Ok(em)) => { - Ok(((ep.value() - em.value()) / 2i32).with_precision(self.precision)) - } - // |x| large enough that exp_m1 overflowed: sinh(x) → ±inf (sign of x). - _ => Err(FpError::Overflow(x.sign())), + // Hoist the exp overflow out of the Ziv closure (it can't return Err): sinh(±huge) = ±inf. + if exp_overflows::(self, x, &mut cache) { + return Err(FpError::Overflow(x.sign())); } + + // sinh(x) = (exp_m1(x) - exp_m1(-x)) / 2 (cancellation-free). `exp_m1` is itself Ziv-correct + // at the working precision, so only the subtraction/divide rounding contributes to the + // radius (a few working-ULPs, scaled by the `exp_m1(x) ≈ 2·sinh(x)` magnitude ratio). + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let ep = work + .exp_m1(&x_f.repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let em = work + .exp_m1(&(-x_f.clone()).repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let result = (ep - em) / 2i32; + let radius = result.ulp() * 12; + (result, radius) + })) } /// Hyperbolic cosine. @@ -71,18 +82,29 @@ impl Context { return Ok(Exact(FBig::new(Repr::one(), *self))); } - // cosh(x) = (exp_m1(x) + exp_m1(-x)) / 2 + 1 (no cancellation: same-sign sum) - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let neg_x = -x_f.clone(); - let ep = work.exp_m1(&x_f.repr, reborrow_cache(&mut cache)); - let em = work.exp_m1(&neg_x.repr, reborrow_cache(&mut cache)); - match (ep, em) { - (Ok(ep), Ok(em)) => Ok(((ep.value() + em.value()) / 2i32 + FBig::::ONE) - .with_precision(self.precision)), - // cosh(x) ≥ 0 always, so overflow → +inf regardless of x's sign. - _ => Err(FpError::Overflow(Sign::Positive)), + // Hoist the exp overflow out of the Ziv closure: cosh(±huge) = +inf (always positive). + if exp_overflows::(self, x, &mut cache) { + return Err(FpError::Overflow(Sign::Positive)); } + // cosh(x) = (exp_m1(x) + exp_m1(-x)) / 2 + 1 (no cancellation: same-sign sum). `exp_m1` is + // Ziv-correct at the working precision; the radius is a few working-ULPs (the `exp_m1(x) ≈ + // 2·cosh(x)` magnitude ratio, plus the trailing +1). + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let ep = work + .exp_m1(&x_f.repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let em = work + .exp_m1(&(-x_f.clone()).repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let result = (ep + em) / 2i32 + FBig::::ONE; + let radius = result.ulp() * 14; + (result, radius) + })) } /// Simultaneously compute `sinh(x)` and `cosh(x)` (context layer). Returns @@ -109,25 +131,31 @@ impl Context { ); } - // sinh = (exp_m1(x) - exp_m1(-x)) / 2; cosh = (exp_m1(x) + exp_m1(-x)) / 2 + 1 - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let neg_x = -x_f.clone(); - let ep = work.exp_m1(&x_f.repr, reborrow_cache(&mut cache)); - let em = work.exp_m1(&neg_x.repr, reborrow_cache(&mut cache)); - match (ep, em) { - (Ok(ep), Ok(em)) => { - let ep = ep.value(); - let em = em.value(); - let sinh_val = ((ep.clone() - em.clone()) / 2i32).with_precision(self.precision); - let cosh_val = - ((ep + em) / 2i32 + FBig::::ONE).with_precision(self.precision); - (Ok(sinh_val), Ok(cosh_val)) - } - // |x| large enough that exp_m1 overflowed: - // sinh(x) → ±inf (sign of x), cosh(x) → +inf - _ => (Err(FpError::Overflow(x.sign())), Err(FpError::Overflow(Sign::Positive))), + // Hoist the exp overflow out of the Ziv closure: sinh(±huge) = ±inf, cosh(±huge) = +inf. + if exp_overflows::(self, x, &mut cache) { + return (Err(FpError::Overflow(x.sign())), Err(FpError::Overflow(Sign::Positive))); } + // sinh = (ep - em)/2; cosh = (ep + em)/2 + 1, sharing the two `exp_m1` calls. Certified as a + // pair via `ziv_pair` (retry while either endpoint straddles a boundary). + let initial_guard = self.base_guard_digits::() + 10; + let (sinh_r, cosh_r) = self.ziv_pair(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let ep = work + .exp_m1(&x_f.repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let em = work + .exp_m1(&(-x_f.clone()).repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let sinh_val = (ep.clone() - em.clone()) / 2i32; + let cosh_val = (ep + em) / 2i32 + FBig::::ONE; + let sinh_radius = sinh_val.ulp() * 12; + let cosh_radius = cosh_val.ulp() * 14; + ((sinh_val, sinh_radius), (cosh_val, cosh_radius)) + }); + (Ok(sinh_r), Ok(cosh_r)) } /// Hyperbolic tangent. @@ -151,20 +179,25 @@ impl Context { return Ok(Exact(FBig::new(signed_zero_repr(x), *self))); } - // tanh(x) = exp_m1(2x) / (exp_m1(2x) + 2). For large negative x, exp_m1(2x) → -1, - // giving -1/(-1+2) = -1; for large positive x, exp_m1(2x) overflows, but - // tanh(+huge) = 1, so short-circuit before the division would yield +inf/+inf = NaN. - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let two_x = x_f.clone() * 2i32; - match work.exp_m1(&two_x.repr, reborrow_cache(&mut cache)) { - Err(FpError::Overflow(_)) => Ok(FBig::ONE.with_precision(self.precision)), - Ok(e) => { - let e = e.value(); - Ok((e.clone() / (e.clone() + 2i32)).with_precision(self.precision)) + // tanh(x) = exp_m1(2x) / (exp_m1(2x) + 2). `exp_m1(2x)` is Ziv-correct at the working + // precision. For large positive x it overflows → tanh = +1 (returned inline as an exact + // value); for large negative x, exp_m1(2x) → -1 (finite), so tanh → -1 naturally. + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let two_x = x_f * 2i32; + match work.exp_m1(&two_x.repr, reborrow_cache(&mut cache)) { + Err(FpError::Overflow(_)) => (FBig::::ONE, FBig::::ZERO), // exact +1 + Ok(e) => { + let e = e.value(); + let result = e.clone() / (e + 2i32); + let radius = result.ulp() * 12; + (result, radius) + } + Err(other) => unreachable!("exp_m1 on finite input: {other:?}"), } - Err(other) => Err(other), - } + })) } /// Inverse hyperbolic sine. @@ -182,31 +215,40 @@ impl Context { return Ok(Exact(FBig::new(signed_zero_repr(x), *self))); } - // asinh(x) = sign(x) · ln_1p(|x| + x²/(sqrt(x²+1)+1)). - // The x²/(sqrt+1) form avoids the `sqrt(x²+1) - 1` cancellation near 0. - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let sign = x_f.sign(); - let abs_x = x_f.abs(); - let arg = match work.sqr(&abs_x.repr) { - Ok(x_sq) => { - let x_sq = x_sq.value(); - let sqrt_plus_one = work.sqrt(&(x_sq.clone() + FBig::::ONE).repr)?.value() - + FBig::::ONE; - abs_x.clone() + x_sq / sqrt_plus_one - } - // |x| so large that x² overflows: asinh(x) ≈ sign·ln(2|x|) (the √(1+1/x²) - // correction is far below representable precision here). - Err(FpError::Overflow(_)) => { - let ln_val = work - .ln(&(abs_x.clone() * 2i32).repr, reborrow_cache(&mut cache))? - .value(); - return Ok(apply_sign(ln_val, sign).with_precision(self.precision)); - } - Err(other) => return Err(other), - }; - let res = work.ln_1p(&arg.repr, reborrow_cache(&mut cache))?.value(); - Ok(apply_sign(res, sign).with_precision(self.precision)) + // asinh(x) = sign(x) · ln_1p(|x| + x²/(sqrt(x²+1)+1)) — the x²/(sqrt+1) form avoids the + // `sqrt(x²+1) − 1` cancellation near 0. `ln_1p`/`ln`/`sqrt` are Ziv-correct at the working + // precision, so the radius is a few working-ULPs of accumulated arithmetic. The `|x|` so + // large that `x²` overflows arm falls back to the asymptotic `sign·ln(2|x|)`. + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let sign = x_f.sign(); + let abs_x = x_f.abs(); + let res = match work.sqr(&abs_x.repr) { + Ok(x_sq) => { + let x_sq = x_sq.value(); + let sqrt_plus_one = work + .sqrt(&(x_sq.clone() + FBig::::ONE).repr) + .unwrap() + .value() + + FBig::::ONE; + let arg = abs_x.clone() + x_sq / sqrt_plus_one; + work.ln_1p(&arg.repr, reborrow_cache(&mut cache)) + .unwrap() + .value() + } + // |x| so large that x² overflows: asinh(x) ≈ sign·ln(2|x|). + Err(FpError::Overflow(_)) => work + .ln(&(abs_x.clone() * 2i32).repr, reborrow_cache(&mut cache)) + .unwrap() + .value(), + Err(other) => unreachable!("sqr: {other:?}"), + }; + let result = apply_sign(res, sign); + let radius = result.ulp() * 14; + (result, radius) + })) } /// Inverse hyperbolic cosine. Domain: `x ≥ 1`. @@ -234,25 +276,33 @@ impl Context { return Ok(Exact(FBig::new(Repr::zero(), *self))); } - // acosh(x) = ln_1p((x-1) + sqrt((x-1)(x+1))). The (x-1)(x+1) form avoids the - // `x²-1` cancellation near x = 1. - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let xm1 = &x_f - FBig::::ONE; - let xp1 = &x_f + FBig::::ONE; - let arg = match work.mul(&xm1.repr, &xp1.repr) { - Ok(prod) => xm1.clone() + work.sqrt(&prod.value().repr)?.value(), - // (x-1)(x+1) overflowed: acosh(x) ≈ ln(2x). - Err(FpError::Overflow(_)) => { - let ln_val = work - .ln(&(x_f.clone() * 2i32).repr, reborrow_cache(&mut cache))? - .value(); - return Ok(ln_val.with_precision(self.precision)); - } - Err(other) => return Err(other), - }; - let res = work.ln_1p(&arg.repr, reborrow_cache(&mut cache))?.value(); - Ok(res.with_precision(self.precision)) + // acosh(x) = ln_1p((x-1) + sqrt((x-1)(x+1))) — the (x-1)(x+1) form avoids the `x²−1` + // cancellation near x = 1. `ln_1p`/`ln`/`sqrt` are Ziv-correct at the working precision; + // the radius is a few working-ULPs (generous for the near-x=1 cancellation). The `(x-1)(x+1)` + // overflow arm falls back to the asymptotic `ln(2x)`. + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let xm1 = &x_f - FBig::::ONE; + let xp1 = &x_f + FBig::::ONE; + let res = match work.mul(&xm1.repr, &xp1.repr) { + Ok(prod) => { + let arg = xm1.clone() + work.sqrt(&prod.value().repr).unwrap().value(); + work.ln_1p(&arg.repr, reborrow_cache(&mut cache)) + .unwrap() + .value() + } + // (x-1)(x+1) overflowed: acosh(x) ≈ ln(2x). + Err(FpError::Overflow(_)) => work + .ln(&(x_f.clone() * 2i32).repr, reborrow_cache(&mut cache)) + .unwrap() + .value(), + Err(other) => unreachable!("mul: {other:?}"), + }; + let radius = res.ulp() * 16; + (res, radius) + })) } /// Inverse hyperbolic tangent. Domain: `-1 < x < 1` (`x = ±1` → ±∞, `|x| > 1` is an error). @@ -278,16 +328,26 @@ impl Context { _ => {} } - // atanh(x) = ln_1p(2x/(1-x)) / 2. - let work = Context::::new(self.precision + 50); - let x_f = FBig::::new(work.repr_round_ref(x).value(), work); - let ratio = (x_f.clone() * 2i32) / (FBig::::ONE - &x_f); - let res = work.ln_1p(&ratio.repr, reborrow_cache(&mut cache))?.value(); - Ok((res / 2i32).with_precision(self.precision)) + // atanh(x) = ln_1p(2x/(1-x)) / 2. `ln_1p` is Ziv-correct at the working precision; the + // radius is a few working-ULPs (generous: the `2x/(1-x)` division amplifies as |x| → 1, but + // the result grows there too, so its ULP keeps the bound sound — Ziv retries near |x|=1). + let initial_guard = self.base_guard_digits::() + 10; + Ok(self.ziv(initial_guard, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let ratio = (x_f.clone() * 2i32) / (FBig::::ONE - &x_f); + let res = work + .ln_1p(&ratio.repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let result = res / 2i32; + let radius = result.ulp() * 16; + (result, radius) + })) } } -impl FBig { +impl FBig { /// Calculate the hyperbolic sine of the floating point number. /// /// # Examples diff --git a/float/src/math/trig.rs b/float/src/math/trig.rs index e9d931b6..76a16096 100644 --- a/float/src/math/trig.rs +++ b/float/src/math/trig.rs @@ -14,13 +14,27 @@ use crate::{ FpResult, }, repr::{Context, Repr, Word}, - round::{Round, Rounded}, + round::{ErrorBounds, Round, Rounded}, }; -use core::cmp::Ordering; use core::convert::TryFrom; -use dashu_base::{AbsOrd, Approximation::Exact, RemEuclid, Sign}; +use dashu_base::{Abs, AbsOrd, Approximation::Exact, RemEuclid, Sign, UnsignedAbs}; use dashu_int::IBig; +/// A near-correct value paired with its provable error radius (the Ziv closure contract). +pub(crate) type Rad = (FBig, FBig); + +/// Series-truncation error radius shared by the Maclaurin/Euler cores (`sin`/`cos`/`sin_cos`/ +/// `atan` here, and `ln`). Each accumulated term contributes `< 1 ulp` of rounding and the +/// truncated tail adds another `< 1 ulp`, so `|value − true| < (4·terms + 12)·ulp(value)`: the +/// `4·terms` covers per-step rounding, the `12` the reconstruction (the `×2` atanh factor, the +/// `s·ln2`/powering recombination, and a safety margin). +pub(crate) fn series_radius( + value: &FBig, + terms: usize, +) -> FBig { + value.ulp() * (4 * terms + 12) +} + #[derive(Debug, Clone, Copy, PartialEq, Eq)] enum Quadrant { First, @@ -43,19 +57,14 @@ fn signed_zero_normal( Ok(Exact(FBig::::new(zero, *ctx))) } -impl Context { - /// Calculate the internal work context for trigonometric functions based on input magnitude. - /// - /// This ensures we have enough guard digits to prevent catastrophic cancellation - /// during range reduction for large inputs. - fn compute_work_context_trig(self, x: &Repr) -> Self { +impl Context { + /// Work context for trigonometric functions: enough guard digits to absorb the catastrophic + /// cancellation in `x − k·(π/2)` for large `|x|`. `guard` (the Ziv retry's growing margin) + /// replaces the fixed base; `x_mag/10` covers cumulative reduction error scaling with `|x|`. + fn compute_work_context_trig(self, x: &Repr, guard: usize) -> Self { // x_mag estimates m = floor(log_BASE(|x|)) let x_mag = (x.exponent.saturating_add(x.digits_ub() as isize)).max(0) as usize; - - // We need precision + log10(x) digits to maintain 'precision' digits after reduction. - // We add a base of 50 guard digits, plus 10% of x_mag for very large arguments - // to account for cumulative errors in division and multiplication during reduction. - let extra_guards = 50 + x_mag / 10; + let extra_guards = guard + x_mag / 10; let work_precision = self .precision .saturating_add(x_mag) @@ -63,25 +72,37 @@ impl Context { Self::new(work_precision) } - /// Reduces the argument to the first quadrant for trigonometric evaluation. - /// Returns the internal work context, the reduced argument `r`, and the quadrant `k % 4`. + /// Reduces the argument to the first quadrant: `r = x − k·(π/2)` with `r ∈ (−π/4, π/4]`. + /// Returns the work context, `r`, the quadrant `k % 4`, and the **reduction error** — a provable + /// bound on `|r_computed − r_true|` (dominated by `|k|·ulp(half_pi)` for huge `|x|`), which the + /// Ziv wrapper folds into the result radius so the containment test is sound. fn reduce_to_quadrant( self, x: &Repr, + guard: usize, mut cache: Option<&mut ConstCache>, - ) -> (Self, FBig, Quadrant) { - let work_context = self.compute_work_context_trig(x); + ) -> (Self, FBig, Quadrant, FBig) { + let work_context = self.compute_work_context_trig(x, guard); let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); let pi = work_context.pi::(reborrow_cache(&mut cache)).value(); - let half_pi = &pi / 2; + let half_pi = &pi / 2u8; let x_scaled: FBig = &x_f / &half_pi; let k_f = x_scaled.round(); - let r = x_f - &k_f * half_pi; + let r = x_f - &k_f * &half_pi; // `k_f` is the integer nearest `x_scaled`, so it's exact (or a signed zero // for a tiny argument in (-1, 0), which `IBig::try_from` treats as plain 0). let k = IBig::try_from(k_f).expect("k_f is an exact integer or signed zero"); + // Reduction error bound: the rounded `half_pi` carries `< 1 ulp`, scaled by `|k|`; the `x` + // rounding and the subtraction add a few `r`-ULPs. Computed at the work precision (|k| fits + // in its digits, so this is accurate; the full-ulp factors and the +4 over-estimate) — kept + // off unlimited precision so `tan`'s `/|cos|` radius division stays legal. + let half_pi_ulp = half_pi.ulp(); + let r_ulp = r.ulp(); + let k_abs = k.clone().unsigned_abs(); + let reduction_err = half_pi_ulp * k_abs + r_ulp * 4; + let k_mod_4_big = k.rem_euclid(IBig::from(4)); let Ok(k_mod_4_int) = i8::try_from(k_mod_4_big) else { unreachable!("k % 4 is always in [0, 3]"); @@ -94,7 +115,7 @@ impl Context { _ => unreachable!(), }; - (work_context, r, quadrant) + (work_context, r, quadrant, reduction_err) } /// Calculate the sine of the floating point representation. @@ -107,28 +128,39 @@ impl Context { return Err(FpError::InfiniteInput); } assert_limited_precision(self.precision); - if x.significand.is_zero() { // sin(±0) = ±0 return signed_zero_normal(self, x); } - let (work_context, r, quadrant) = self.reduce_to_quadrant(x, reborrow_cache(&mut cache)); - - // 3. Evaluate the reduced series based on the quadrant - let res = match quadrant { - Quadrant::First => work_context.sin_internal(&r), - Quadrant::Second => work_context.cos_internal(&r), - Quadrant::Third => -work_context.sin_internal(&r), - Quadrant::Fourth => -work_context.cos_internal(&r), - }; - Ok(res.with_precision(self.precision)) + // Ziv: reduce to the first quadrant (the guard grows per retry, enlarging the work precision + // that absorbs the `x − k·(π/2)` cancellation), evaluate the series, and fold the reduction + // error into the radius so the containment test is sound even for huge |x|. + Ok(self.ziv(50, |guard| { + let (work, r, quadrant, reduction_err) = + self.reduce_to_quadrant(x, guard, reborrow_cache(&mut cache)); + let (val, series_radius) = match quadrant { + Quadrant::First => work.sin_compute(&r), + Quadrant::Second => work.cos_compute(&r), + Quadrant::Third => { + let (v, e) = work.sin_compute(&r); + (-v, e) + } + Quadrant::Fourth => { + let (v, e) = work.cos_compute(&r); + (-v, e) + } + }; + (val, series_radius + reduction_err) + })) } - /// Internal Taylor series for sine: S(x) = x - x^3/3! + x^5/5! - ... - fn sin_internal(self, x: &FBig) -> FBig { + /// Near-correct sine series `S(x) = x − x³/3! + x⁵/5! − …` on the reduced argument, returning + /// `(value, error_radius)`. The radius covers series truncation (`< 1 working-ULP` by the break + /// test) plus `~3K` steps of rounding accumulation. Used by the Ziv-backed `sin`/`cos`/`tan`. + fn sin_compute(self, x: &FBig) -> (FBig, FBig) { if x.repr.significand.is_zero() { - return FBig::ZERO; + return (FBig::ZERO, FBig::ZERO); } let x2 = x.sqr(); let mut sum = x.clone(); @@ -148,7 +180,8 @@ impl Context { } k += 1; } - sum + let radius = series_radius(&sum, k); + (sum, radius) } /// Calculate the cosine of the floating point representation. @@ -167,22 +200,30 @@ impl Context { return Ok(FBig::::ONE.with_precision(self.precision)); } - let (work_context, r, quadrant) = self.reduce_to_quadrant(x, reborrow_cache(&mut cache)); - - // 3. Evaluate the reduced series based on the quadrant - let res = match quadrant { - Quadrant::First => work_context.cos_internal(&r), - Quadrant::Second => -work_context.sin_internal(&r), - Quadrant::Third => -work_context.cos_internal(&r), - Quadrant::Fourth => work_context.sin_internal(&r), - }; - Ok(res.with_precision(self.precision)) + Ok(self.ziv(50, |guard| { + let (work, r, quadrant, reduction_err) = + self.reduce_to_quadrant(x, guard, reborrow_cache(&mut cache)); + let (val, series_radius) = match quadrant { + Quadrant::First => work.cos_compute(&r), + Quadrant::Second => { + let (v, e) = work.sin_compute(&r); + (-v, e) + } + Quadrant::Third => { + let (v, e) = work.cos_compute(&r); + (-v, e) + } + Quadrant::Fourth => work.sin_compute(&r), + }; + (val, series_radius + reduction_err) + })) } - /// Internal Taylor series for cosine: C(x) = 1 - x^2/2! + x^4/4! - ... - fn cos_internal(self, x: &FBig) -> FBig { + /// Near-correct cosine series `C(x) = 1 − x²/2! + x⁴/4! − …`, returning `(value, radius)`. + /// (See [`sin_compute`](Self::sin_compute) for the radius derivation.) + fn cos_compute(self, x: &FBig) -> (FBig, FBig) { if x.repr.significand.is_zero() { - return FBig::ONE.with_precision(self.precision).value(); + return (FBig::ONE.with_precision(self.precision).value(), FBig::ZERO); } let x2 = x.sqr(); let mut sum = FBig::::ONE.with_precision(self.precision).value(); @@ -202,7 +243,8 @@ impl Context { } k += 1; } - sum + let radius = series_radius(&sum, k); + (sum, radius) } /// Calculate both the sine and cosine of the floating point representation. @@ -225,27 +267,28 @@ impl Context { return (s, c); } - let (work_context, r, quadrant) = self.reduce_to_quadrant(x, reborrow_cache(&mut cache)); - - let (sin_r, cos_r) = work_context.sin_cos_internal(&r); - - let (s, c) = match quadrant { - Quadrant::First => (sin_r, cos_r), - Quadrant::Second => (cos_r, -sin_r), - Quadrant::Third => (-sin_r, -cos_r), - Quadrant::Fourth => (-cos_r, sin_r), - }; - - (Ok(s.with_precision(self.precision)), Ok(c.with_precision(self.precision))) + let (s, c) = self.ziv_pair(50, |guard| { + let (work, r, quadrant, reduction_err) = + self.reduce_to_quadrant(x, guard, reborrow_cache(&mut cache)); + let ((sin_r, sin_e), (cos_r, cos_e)) = work.sin_cos_compute(&r); + let (s, c) = match quadrant { + Quadrant::First => (sin_r, cos_r), + Quadrant::Second => (cos_r, -sin_r), + Quadrant::Third => (-sin_r, -cos_r), + Quadrant::Fourth => (-cos_r, sin_r), + }; + ((s, sin_e + reduction_err.clone()), (c, cos_e + reduction_err)) + }); + (Ok(s), Ok(c)) } - /// Simultaneously evaluate Taylor series for sine and cosine. - pub(crate) fn sin_cos_internal( - self, - x: &FBig, - ) -> (FBig, FBig) { + /// Simultaneously evaluate the sine and cosine series, returning both values and their radii. + pub(crate) fn sin_cos_compute(self, x: &FBig) -> (Rad, Rad) { if x.repr.significand.is_zero() { - return (FBig::ZERO, FBig::ONE.with_precision(self.precision).value()); + return ( + (FBig::ZERO, FBig::ZERO), + (FBig::ONE.with_precision(self.precision).value(), FBig::ZERO), + ); } let x2 = x.sqr(); let mut sin_sum = x.clone(); @@ -275,13 +318,17 @@ impl Context { } k += 1; } - (sin_sum, cos_sum) + ( + (sin_sum.clone(), series_radius(&sin_sum, k)), + (cos_sum.clone(), series_radius(&cos_sum, k)), + ) } /// Calculate the tangent of the floating point representation. /// /// # Note - /// Near odd multiples of π/2, the result is an infinity (returned as a value, not an error). + /// Near odd multiples of π/2 the value grows without bound; dashu's wide exponent range holds + /// it as a large finite number rather than saturating to ±∞. pub fn tan( &self, x: &Repr, @@ -297,27 +344,38 @@ impl Context { return signed_zero_normal(self, x); } - let (work_context, r, quadrant) = self.reduce_to_quadrant(x, reborrow_cache(&mut cache)); - let (sin_r, cos_r) = work_context.sin_cos_internal(&r); - - let (s_f, c_f) = match quadrant { - Quadrant::First => (sin_r, cos_r), - Quadrant::Second => (cos_r, -sin_r), - Quadrant::Third => (-sin_r, -cos_r), - Quadrant::Fourth => (-cos_r, sin_r), - }; - - if c_f.repr.is_pos_zero() { - // tan hits a pole: the result is an infinity with the sign of the numerator. - let inf = if s_f.sign() == Sign::Negative { - Repr::neg_infinity() - } else { - Repr::infinity() + // tan = sin/cos, correctly rounded via the Ziv loop. Near a pole (an odd multiple of π/2) + // the value is large but finite at the working precision — dashu's wide exponent range holds + // it, and the sign is carried by the arithmetic (s/−|c| is negative), so there is no pole + // special-case here. The closure's `significand.is_zero()` guard below handles the + // unreachable exact-pole case (cos cancelling to a zero significand) by forcing a retry. + // Skipping a hoisted pole check avoids recomputing the sin/cos series twice (once for the + // check, once for the first Ziv attempt). + Ok(self.ziv(50, |guard| { + let (work, r, quadrant, reduction_err) = + self.reduce_to_quadrant(x, guard, reborrow_cache(&mut cache)); + let ((sin_r, sin_e), (cos_r, cos_e)) = work.sin_cos_compute(&r); + let (s, c) = match quadrant { + Quadrant::First => (sin_r, cos_r), + Quadrant::Second => (cos_r, -sin_r), + Quadrant::Third => (-sin_r, -cos_r), + Quadrant::Fourth => (-cos_r, sin_r), }; - return Ok(Rounded::Exact(FBig::new(inf, *self))); - } - self.div(&s_f.repr, &c_f.repr) - .map(|r| r.and_then(|f| f.with_precision(self.precision))) + if c.repr.significand.is_zero() { + // cos rounded to a zero significand at this guard (the input sits on a work- + // precision pole — unreachable for finite-precision x): force a retry. A higher guard + // makes cos representable (nonzero), yielding a large finite tan. + return (FBig::ZERO, FBig::ONE); + } + let result = work.div(&s.repr, &c.repr).unwrap().value(); + // tan = s/c: the sin/cos radii propagate as (e_s + |tan|·e_c)/|c| plus the division + // rounding, all at the working precision (the only term that needed unlimited precision + // — the reduction error — is already work-precision, so the `/|c|` stays legal). + let e_s = sin_e + reduction_err.clone(); + let e_c = cos_e + reduction_err; + let radius = (e_s + result.clone().abs() * e_c) / c.clone().abs() + result.ulp() * 8; + (result, radius) + })) } /// Calculate the arcsine of the floating point representation. @@ -334,6 +392,12 @@ impl Context { return Err(FpError::InfiniteInput); } assert_limited_precision(self.precision); + if x.significand.is_zero() { + // asin(±0) = ±0 (asin is odd), exact. Like the other inverse trig/hyperbolic functions, + // short-circuit before the Ziv loop: a zero result carries a positive radius that can't + // be certified against 0's one-sided preimage under directed rounding. + return signed_zero_normal(self, x); + } let x_orig = FBig::::new(x.clone(), *self); // Domain check: |x| must be <= 1 @@ -341,39 +405,37 @@ impl Context { return Err(FpError::OutOfDomain); } - let guard_digits = 50; - let work_precision = self.precision + guard_digits; - let work_context = Self::new(work_precision); - - let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); - - let res = work_context.asin_internal(&x_f, reborrow_cache(&mut cache)); - Ok(res.with_precision(self.precision)) - } - - fn asin_internal( - self, - x_f: &FBig, - mut cache: Option<&mut ConstCache>, - ) -> FBig { - let one = FBig::::ONE.with_precision(self.precision).value(); - let x2 = x_f.sqr(); - let d = self.unwrap_fp(self.sqrt(&(one - x2).repr)); - - if d.repr.is_pos_zero() || d.repr.is_neg_zero() { - // |x| = 1 exactly (d = sqrt(1 - x²) = ±0); asin(±1) = ±π/2 regardless of rounding - // mode. Catch `-0` too: under roundTowardNegative `1 - 1` cancels to `-0`, sqrt(-0) = - // `-0`, and the general path would divide by `-0` → `-∞` and panic. (For |x| < 1, - // `d` is strictly positive, so only the exact-boundary x reaches this branch.) - let pi = self.pi::(reborrow_cache(&mut cache)).value(); - let half_pi: FBig = pi / 2; - if x_f.sign() == Sign::Positive { - return half_pi; + Ok(self.ziv(50, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let one = FBig::::ONE.with_precision(work.precision).value(); + let d = work + .sqrt(&(one.clone() - x_f.clone().sqr()).repr) + .unwrap() + .value(); + if d.repr.is_pos_zero() || d.repr.is_neg_zero() { + // |x| = 1: asin(±1) = ±π/2. + let pi = work.pi::(reborrow_cache(&mut cache)).value(); + let half_pi = pi / 2u8; + let res = if x_f.sign() == Sign::Positive { + half_pi + } else { + -half_pi + }; + let radius = res.ulp() * 4; + return (res, radius); } - return -half_pi; - } - - self.atan_with_reduction(&(x_f / d), reborrow_cache(&mut cache)) + // asin(x) = atan(x / sqrt(1−x²)); `atan`/`sqrt` are Ziv-correct at the working + // precision, so the radius is just the accumulated `sqrt`+`div` rounding (well-conditioned + // near |x|=1, where atan's derivative → 0). + let arg = &x_f / &d; + let res = work + .atan(&arg.repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let radius = res.ulp() * 16; + (res, radius) + })) } /// Calculate the arccosine of the floating point representation. @@ -392,22 +454,34 @@ impl Context { assert_limited_precision(self.precision); let x_orig = FBig::::new(x.clone(), *self); - // Domain check: |x| must be <= 1 - if x_orig.abs_cmp(&FBig::ONE).is_gt() { + let cmp_one = x_orig.abs_cmp(&FBig::ONE); + if cmp_one.is_gt() { return Err(FpError::OutOfDomain); } + if cmp_one.is_eq() { + // |x| = 1: the composition π/2 − asin(±1) cancels onto an exact value. acos(1) = 0 is + // the acute case — under directed rounding 0's preimage is one-sided ([0, ulp)), so the + // Ziv containment test can never certify it (any positive radius dips the interval below + // 0) and would infinite-retry. acos(-1) = π is handled here too, for symmetry. + return Ok(if x.sign() == Sign::Positive { + Exact(FBig::::new(Repr::zero(), *self)) + } else { + self.pi::(reborrow_cache(&mut cache)) + }); + } - let guard_digits = 50; - let work_precision = self.precision + guard_digits; - let work_context = Self::new(work_precision); - - let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); - - let asin_x = work_context.asin_internal(&x_f, reborrow_cache(&mut cache)); - let pi = work_context.pi::(reborrow_cache(&mut cache)).value(); - let half_pi: FBig = pi / 2; - let res: FBig = half_pi - asin_x; - Ok(res.with_precision(self.precision)) + Ok(self.ziv(50, |guard| { + let work = Context::::new(self.precision + guard); + // acos(x) = π/2 − asin(x); `asin`/`pi` are Ziv-correct (or exact) at the working + // precision. The radius covers the propagated asin/π rounding plus the subtraction, + // which cancels near x = 1 — the radius grows there and Ziv retries with more guard. + let asin_x = work.asin(x, reborrow_cache(&mut cache)).unwrap().value(); + let pi = work.pi::(reborrow_cache(&mut cache)).value(); + let res = (pi / 2u8) - &asin_x; + let radius = asin_x.ulp().clone().with_precision(0).value() * 2 + + res.ulp().clone().with_precision(0).value() * 4; + (res, radius) + })) } /// Calculate the arctangent of the floating point representation. @@ -435,42 +509,31 @@ impl Context { return signed_zero_normal(self, x); } - let guard_digits = 50; - let work_precision = self.precision + guard_digits; - let work_context = Self::new(work_precision); - - let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); - let res = work_context.atan_with_reduction(&x_f, reborrow_cache(&mut cache)); - Ok(res.with_precision(self.precision)) - } - - /// Internal arctangent that includes range reduction but no guard digit allocation. - fn atan_with_reduction( - self, - x_f: &FBig, - mut cache: Option<&mut ConstCache>, - ) -> FBig { - let sign = x_f.sign(); - let mut x_abs = x_f.clone(); - if sign == Sign::Negative { - x_abs = -x_abs; - } - let mut res = if x_abs >= FBig::::ONE.with_precision(self.precision).value() { - let pi = self.pi::(reborrow_cache(&mut cache)).value(); - let inv_x = FBig::::ONE.with_precision(self.precision).value() / x_abs; - (pi / 2) - self.atan_internal(&inv_x) - } else { - self.atan_internal(&x_abs) - }; - if sign == Sign::Negative { - res = -res; - } - res + Ok(self.ziv(50, |guard| { + let work = Context::::new(self.precision + guard); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let sign = x_f.sign(); + let x_abs = x_f.abs(); + let one = FBig::::ONE.with_precision(work.precision).value(); + let (res, radius) = if x_abs >= one { + // |x| ≥ 1: atan(x) = π/2 − atan(1/x); the series runs on 1/x ∈ (0, 1]. + let pi = work.pi::(reborrow_cache(&mut cache)).value(); + let inv_x = &one / &x_abs; + let (atan_val, atan_radius) = work.atan_compute(&inv_x); + let res = (pi / 2u8) - atan_val; + let radius = atan_radius + res.ulp() * 4; + (res, radius) + } else { + work.atan_compute(&x_abs) + }; + let res = if sign == Sign::Negative { -res } else { res }; + (res, radius) + })) } - /// Internal series for arctangent. - /// Evaluates the Euler series for arctangent. - fn atan_internal(self, x: &FBig) -> FBig { + /// Near-correct Euler series for `atan(x)` (`|x| ≤ 1`), returning `(value, radius)`. The radius + /// covers series truncation plus `~3N` steps of accumulation. + fn atan_compute(self, x: &FBig) -> (FBig, FBig) { // Euler's series for atan(x) let x2 = x.sqr(); let one_plus_x2 = FBig::ONE + &x2; @@ -489,7 +552,8 @@ impl Context { sum += &term; n += 1; } - sum + let radius = series_radius(&sum, n); + (sum, radius) } /// Calculate the arctangent of y / x. @@ -508,96 +572,72 @@ impl Context { assert_limited_precision(self.precision); - let guard_digits = 50; - let work_precision = self.precision + guard_digits; - let work_context = Self::new(work_precision); - - // Handle Infinities according to IEEE 754 + // Handle Infinities according to IEEE 754 (computed at the target precision). if y.is_infinite() || x.is_infinite() { let (sy, sx) = (y.sign() == Sign::Positive, x.sign() == Sign::Positive); + let pi_val = self.pi::(reborrow_cache(&mut cache)).value(); let res: FBig = match (y.is_infinite(), x.is_infinite(), sy, sx) { - (true, true, true, true) => { - work_context.pi::(reborrow_cache(&mut cache)).value() / 4 - } - (true, true, true, false) => { - work_context.pi::(reborrow_cache(&mut cache)).value() * 3 / 4 - } - (true, true, false, true) => { - let pi4: FBig = - work_context.pi::(reborrow_cache(&mut cache)).value() / 4; - -pi4 - } - (true, true, false, false) => { - let pi34: FBig = - work_context.pi::(reborrow_cache(&mut cache)).value() * 3 / 4; - -pi34 - } - (true, false, true, _) => { - work_context.pi::(reborrow_cache(&mut cache)).value() / 2 - } - (true, false, false, _) => { - let half_pi: FBig = - work_context.pi::(reborrow_cache(&mut cache)).value() / 2; - -half_pi - } + (true, true, true, true) => pi_val.clone() / 4u8, + (true, true, true, false) => pi_val.clone() * 3u8 / 4u8, + (true, true, false, true) => -(pi_val.clone() / 4u8), + (true, true, false, false) => -(pi_val.clone() * 3u8 / 4u8), + (true, false, true, _) => pi_val.clone() / 2u8, + (true, false, false, _) => -(pi_val.clone() / 2u8), (false, true, _, true) => { // atan2(±finite, +inf) = ±0 (signed zero of y) if sy { - FBig::::ZERO.with_precision(work_precision).value() + FBig::::ZERO } else { - FBig::::new(Repr::neg_zero(), work_context) - .with_precision(work_precision) - .value() + FBig::::new(Repr::neg_zero(), *self) } } - (false, true, true, false) => { - work_context.pi::(reborrow_cache(&mut cache)).value() - } - (false, true, false, false) => { - -work_context.pi::(reborrow_cache(&mut cache)).value() - } + (false, true, true, false) => pi_val.clone(), + (false, true, false, false) => -pi_val, _ => unreachable!(), }; return Ok(res.with_precision(self.precision)); } - let y_f = FBig::::new(work_context.repr_round(y.clone()).value(), work_context); - let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + // x == 0, y finite nonzero: atan2 = ±π/2. + if x.significand.is_zero() { + let half_pi = self.pi::(reborrow_cache(&mut cache)).value() / 2u8; + let res = if y.sign() == Sign::Positive { + half_pi + } else { + -half_pi + }; + return Ok(res.with_precision(self.precision)); + } - match x_f.cmp(&FBig::::ZERO) { - Ordering::Greater => { - let res = - work_context.atan_with_reduction(&(y_f / x_f), reborrow_cache(&mut cache)); - Ok(res.with_precision(self.precision)) - } - Ordering::Less => { - let pi = work_context.pi::(reborrow_cache(&mut cache)).value(); - let y_sign = y_f.sign(); - let atan_yx = - work_context.atan_with_reduction(&(y_f / x_f), reborrow_cache(&mut cache)); - let res = if y_sign == Sign::Positive { - atan_yx + pi + // x ≠ 0, finite: atan2 = atan(y/x) ± (quadrant π). `atan` is Ziv-correct at the working + // precision, so the radius is the accumulated div/π-arithmetic rounding. + Ok(self.ziv(50, |guard| { + let work = Context::::new(self.precision + guard); + let y_f = FBig::::new(work.repr_round_ref(y).value(), work); + let x_f = FBig::::new(work.repr_round_ref(x).value(), work); + let ratio = &y_f / &x_f; + let atan_val = work + .atan(&ratio.repr, reborrow_cache(&mut cache)) + .unwrap() + .value(); + let (res, radius) = if x.sign() == Sign::Positive { + (atan_val.clone(), atan_val.ulp() * 6) + } else { + let pi = work.pi::(reborrow_cache(&mut cache)).value(); + let r = if y_f.sign() == Sign::Positive { + &atan_val + &pi } else { - atan_yx - pi + &atan_val - &pi }; - Ok(res.with_precision(self.precision)) - } - Ordering::Equal => { - // x == 0 case - let pi = work_context.pi::(reborrow_cache(&mut cache)).value(); - let half_pi: FBig = pi / 2; - if y_f > FBig::::ZERO { - Ok(half_pi.with_precision(self.precision)) - } else { - let res = -half_pi; - Ok(res.with_precision(self.precision)) - } - } - } + let radius = atan_val.ulp() * 2 + r.ulp() * 6; + (r, radius) + }; + (res, radius) + })) } } -impl FBig { +impl FBig { /// Calculate the sine of the floating point number. /// /// # Panics @@ -760,4 +800,34 @@ mod tests { // sin(x) ≈ x for a small negative x — completing without panicking is the regression guard. assert_eq!(s.sign(), Sign::Negative); } + + /// tan near a pole (π/2) must not panic, and its sign must follow the pole side: just below → + /// large positive (→ +∞), just above → large negative (→ −∞). Guards the pole check, which + /// tests `cos` with `significand.is_zero()` (not `is_pos_zero`, which would miss `-0`) and + /// assigns the infinity sign as `sign(sin)·sign(cos)`. + #[test] + fn test_tan_near_pole_signs_and_no_panic() { + let p = 53usize; + let ctx = Context::::new(p); + let half_pi = FBig::::pi(p) / 2u8; + // a clear offset either side of the pole (≈2⁻¹⁰, far larger than half_pi's rounding error) + let eps = FBig::::ONE >> 10; + let below = ctx + .tan::<2>((half_pi.clone() - &eps).repr(), None) + .unwrap() + .value(); + let above = ctx + .tan::<2>((half_pi.clone() + &eps).repr(), None) + .unwrap() + .value(); + assert_eq!(below.sign(), Sign::Positive, "tan just below π/2 is large positive"); + assert_eq!(above.sign(), Sign::Negative, "tan just above π/2 is large negative"); + // sanity: tan(π/4) = 1 + let pi = FBig::::pi(p); + let q = ctx.tan::<2>((pi / 4u8).repr(), None).unwrap().value(); + assert!( + (q.clone() - FBig::ONE).abs_cmp(&(FBig::ONE >> 40)).is_le(), + "tan(π/4) ≈ 1, got {q:?}" + ); + } } diff --git a/float/src/mul.rs b/float/src/mul.rs index 10a4543f..44f5530e 100644 --- a/float/src/mul.rs +++ b/float/src/mul.rs @@ -10,53 +10,6 @@ use crate::{ }; use core::ops::{Mul, MulAssign}; -/// Raw product of two finite reprs, attaching the XOR sign of the operands to a zero product -/// (the significand product alone is `+0`, losing the sign). -/// -/// Returns an error when the result exponent overflows or underflows `isize`. -fn make_mul_repr(lhs: &Repr, rhs: &Repr) -> Result, FpError> { - let significand = &lhs.significand * &rhs.significand; - if significand.is_zero() { - return Ok(if lhs.sign() != rhs.sign() { - Repr::neg_zero() - } else { - Repr::zero() - }); - } - let sign = if lhs.sign() != rhs.sign() { - Sign::Negative - } else { - Sign::Positive - }; - let exponent = lhs.exponent.checked_add(rhs.exponent).ok_or_else(|| { - debug_assert!( - lhs.exponent.is_positive() == rhs.exponent.is_positive(), - "checked_add overflow with mixed-sign exponents is impossible" - ); - if lhs.exponent > 0 { - FpError::Overflow(sign) - } else { - FpError::Underflow(sign) - } - })?; - Repr::new(significand, exponent).check_finite_exponent() -} - -macro_rules! unwrap_mul_repr { - ($result:expr, $context:expr) => { - match $result { - Ok(r) => r, - Err(FpError::Overflow(sign)) => { - return FBig::new(Repr::infinity_with_sign(sign), $context); - } - Err(FpError::Underflow(sign)) => { - return FBig::new(Repr::zero_with_sign(sign), $context); - } - Err(_) => unreachable!(), - } - }; -} - impl Mul<&FBig> for &FBig { type Output = FBig; @@ -65,7 +18,10 @@ impl Mul<&FBig> for &FBig { assert_finite_operands(&self.repr, &rhs.repr); let context = Context::max(self.context, rhs.context); - let repr = unwrap_mul_repr!(make_mul_repr(&self.repr, &rhs.repr), context); + let repr = &self.repr * &rhs.repr; + if repr.is_infinite() { + return FBig::new(repr, context); + } FBig::new(context.repr_round(repr).value(), context) } } @@ -78,7 +34,10 @@ impl Mul<&FBig> for FBig { assert_finite_operands(&self.repr, &rhs.repr); let context = Context::max(self.context, rhs.context); - let repr = unwrap_mul_repr!(make_mul_repr(&self.repr, &rhs.repr), context); + let repr = &self.repr * &rhs.repr; + if repr.is_infinite() { + return FBig::new(repr, context); + } FBig::new(context.repr_round(repr).value(), context) } } @@ -91,7 +50,10 @@ impl Mul> for &FBig { assert_finite_operands(&self.repr, &rhs.repr); let context = Context::max(self.context, rhs.context); - let repr = unwrap_mul_repr!(make_mul_repr(&self.repr, &rhs.repr), context); + let repr = &self.repr * &rhs.repr; + if repr.is_infinite() { + return FBig::new(repr, context); + } FBig::new(context.repr_round(repr).value(), context) } } @@ -104,7 +66,10 @@ impl Mul> for FBig { assert_finite_operands(&self.repr, &rhs.repr); let context = Context::max(self.context, rhs.context); - let repr = unwrap_mul_repr!(make_mul_repr(&self.repr, &rhs.repr), context); + let repr = &self.repr * &rhs.repr; + if repr.is_infinite() { + return FBig::new(repr, context); + } FBig::new(context.repr_round(repr).value(), context) } } @@ -205,7 +170,17 @@ impl Context { rhs }; - let repr = make_mul_repr(lhs_repr, rhs_repr)?; + let repr = lhs_repr * rhs_repr; + let repr = if repr.is_infinite() { + return Err(FpError::Overflow(repr.sign())); + } else if repr.significand.is_zero() + && !lhs_repr.significand.is_zero() + && !rhs_repr.significand.is_zero() + { + return Err(FpError::Underflow(repr.sign())); + } else { + repr + }; Ok(self.repr_round(repr).map(|v| FBig::new(v, *self))) } diff --git a/float/src/repr.rs b/float/src/repr.rs index f6f8510b..f85de99d 100644 --- a/float/src/repr.rs +++ b/float/src/repr.rs @@ -1,7 +1,7 @@ use crate::{ error::{assert_finite, FpError}, round::{Round, Rounded}, - utils::{digit_len, split_digits, split_digits_ref}, + utils::{ceil_usize, digit_len, split_digits, split_digits_ref}, }; use core::marker::PhantomData; use dashu_base::{Approximation::*, EstimatedLog2, Sign}; @@ -45,26 +45,6 @@ pub struct Repr { pub(crate) exponent: isize, } -impl PartialEq for Repr { - /// Two representations are equal when they denote the same value. In particular `+0` - /// and `-0` compare equal, as do two infinities of the same sign. - #[inline] - fn eq(&self, other: &Self) -> bool { - if self.significand.is_zero() && other.significand.is_zero() { - let (self_inf, other_inf) = (self.is_infinite(), other.is_infinite()); - match (self_inf, other_inf) { - (true, true) => self.sign() == other.sign(), - (false, false) => true, // both are ±0 - _ => false, // one is zero, the other is infinite - } - } else { - self.significand == other.significand && self.exponent == other.exponent - } - } -} - -impl Eq for Repr {} - /// The context containing runtime information for the floating point number and its operations. /// /// The context currently consists of a *precision limit* and a *rounding mode*. All the operation @@ -316,7 +296,7 @@ impl Repr { /// assert!(Repr::<10>::one().is_int()); /// assert!(!Repr::<16>::new(123.into(), -1).is_int()); /// ``` - pub fn is_int(&self) -> bool { + pub const fn is_int(&self) -> bool { if self.is_infinite() { false } else { @@ -707,6 +687,16 @@ impl Context { self.precision } + /// `⌈log_B(precision)⌉` — the number of base-`B` digits needed to index the precision word. + /// + /// This is the base guard every transcendental Ziv loop adds on top of `precision`: it absorbs + /// the rounding that accumulates over the `O(log p)` series / squaring steps. Each caller adds + /// its own operation-specific constant (`+ 2`, `+ 10`, …); this helper is the shared core that + /// maps the binary precision estimate onto the output base. + pub(crate) fn base_guard_digits(&self) -> usize { + ceil_usize(self.precision.log2_est() / B.log2_est()) + } + /// Round the repr to the desired precision pub(crate) fn repr_round(&self, repr: Repr) -> Rounded> { assert_finite(&repr); diff --git a/float/src/repr_ops.rs b/float/src/repr_ops.rs new file mode 100644 index 00000000..9253795e --- /dev/null +++ b/float/src/repr_ops.rs @@ -0,0 +1,324 @@ +//! Core trait impls for [`Repr`]: [`PartialEq`]/[`Eq`] (value equality — `+0`/`-0` and same-sign +//! infinities compare equal), [`Neg`], and exact [`Add`]/[`Sub`]/[`Mul`]. +//! +//! A [`Repr`] carries no precision limit, so the arithmetic ops are lossless (no rounding). These +//! are the shared primitives the crate reaches for whenever it needs an exact intermediate — the +//! Ziv containment test, the correctly-rounded `Sum`, and the `FBig` multiply path. +//! +//! Each arithmetic operator's logic lives in the `&Repr`-by-`&Repr` primary impl; the val/ref +//! forwarders delegate to it without cloning (a `Repr`'s significand is heap-backed, so the ref/ref +//! form is the no-extra-allocation path). [`Mul`] saturates exponent overflow/underflow to the +//! signed infinity/zero sentinels so the operator is infallible; the precision-limited +//! [`Context`](crate::Context) multiply re-derives the [`FpError`] it needs from that saturated +//! result. + +use core::cmp::Ordering; +use core::ops::{Add, Mul, Neg, Sub}; + +use dashu_base::Sign; + +use crate::error::FpError; +use crate::repr::{Repr, Word}; +use crate::utils::shl_digits; + +impl PartialEq for Repr { + /// Two representations are equal when they denote the same value. In particular `+0` + /// and `-0` compare equal, as do two infinities of the same sign. + #[inline] + fn eq(&self, other: &Self) -> bool { + if self.significand.is_zero() && other.significand.is_zero() { + let (self_inf, other_inf) = (self.is_infinite(), other.is_infinite()); + match (self_inf, other_inf) { + (true, true) => self.sign() == other.sign(), + (false, false) => true, // both are ±0 + _ => false, // one is zero, the other is infinite + } + } else { + self.significand == other.significand && self.exponent == other.exponent + } + } +} + +impl Eq for Repr {} + +impl Neg for Repr { + type Output = Self; + #[inline] + fn neg(self) -> Self::Output { + Repr::neg(self) + } +} + +impl Add<&Repr> for &Repr { + type Output = Repr; + + #[inline] + fn add(self, rhs: &Repr) -> Repr { + debug_assert!(self.is_finite()); + debug_assert!(rhs.is_finite()); + + // Zero operands short-circuit so a `-0` operand's sentinel exponent (-1) can't bleed into + // the result exponent — `precise_sum` (the `Sum` impl) relies on this. The aligned path + // below would compute the same value but would rebuild via `Repr::new`, normalizing away + // the operand's own representation. + if self.significand.is_zero() { + return rhs.clone(); + } + if rhs.significand.is_zero() { + return self.clone(); + } + + // Result exponent is min(lhs, rhs); shift the larger-exponent significand up (appending + // trailing base-`B` zero-digits is lossless) and add. `IBig` add handles opposite signs. + match self.exponent.cmp(&rhs.exponent) { + Ordering::Equal => Repr::new(&self.significand + &rhs.significand, self.exponent), + Ordering::Greater => Repr::new( + shl_digits::(&self.significand, (self.exponent - rhs.exponent) as usize) + + &rhs.significand, + rhs.exponent, + ), + Ordering::Less => Repr::new( + &self.significand + + shl_digits::(&rhs.significand, (rhs.exponent - self.exponent) as usize), + self.exponent, + ), + } + } +} + +impl Add<&Repr> for Repr { + type Output = Repr; + #[inline] + fn add(self, rhs: &Repr) -> Repr { + (&self) + rhs + } +} +impl Add> for &Repr { + type Output = Repr; + #[inline] + fn add(self, rhs: Repr) -> Repr { + self + &rhs + } +} +impl Add> for Repr { + type Output = Repr; + #[inline] + fn add(self, rhs: Repr) -> Repr { + (&self) + &rhs + } +} + +impl Sub<&Repr> for &Repr { + type Output = Repr; + + #[inline] + fn sub(self, rhs: &Repr) -> Repr { + debug_assert!(self.is_finite()); + debug_assert!(rhs.is_finite()); + + // Zero short-circuits: `x - 0 = x` (keep x), `0 - x = -x` (negate x). As with `Add`, this + // keeps a zero operand's own representation rather than rebuilding via `Repr::new`. + if rhs.significand.is_zero() { + return self.clone(); + } + if self.significand.is_zero() { + return rhs.clone().neg(); + } + + // Mirror `Add`: align to the smaller exponent, then subtract significands. + match self.exponent.cmp(&rhs.exponent) { + Ordering::Equal => Repr::new(&self.significand - &rhs.significand, self.exponent), + Ordering::Greater => Repr::new( + shl_digits::(&self.significand, (self.exponent - rhs.exponent) as usize) + - &rhs.significand, + rhs.exponent, + ), + Ordering::Less => Repr::new( + &self.significand + - shl_digits::(&rhs.significand, (rhs.exponent - self.exponent) as usize), + self.exponent, + ), + } + } +} + +impl Sub<&Repr> for Repr { + type Output = Repr; + #[inline] + fn sub(self, rhs: &Repr) -> Repr { + (&self) - rhs + } +} +impl Sub> for &Repr { + type Output = Repr; + #[inline] + fn sub(self, rhs: Repr) -> Repr { + self - &rhs + } +} +impl Sub> for Repr { + type Output = Repr; + #[inline] + fn sub(self, rhs: Repr) -> Repr { + (&self) - &rhs + } +} + +impl Mul<&Repr> for &Repr { + type Output = Repr; + + #[inline] + fn mul(self, rhs: &Repr) -> Repr { + debug_assert!(self.is_finite()); + debug_assert!(rhs.is_finite()); + + let significand = &self.significand * &rhs.significand; + if significand.is_zero() { + // The product significand is `+0`; attach the XOR sign of the operands. + return if self.sign() != rhs.sign() { + Repr::neg_zero() + } else { + Repr::zero() + }; + } + let sign = if self.sign() != rhs.sign() { + Sign::Negative + } else { + Sign::Positive + }; + // Exponent = lhs + rhs; saturate an `isize` overflow to the signed infinity/zero sentinel + // so the operator stays infallible (unreachable for real inputs — it needs operands with + // exponents ~±2^62). `Context::mul` re-derives the `FpError` from this saturated result. + let exponent = match self.exponent.checked_add(rhs.exponent) { + Some(e) => e, + None => { + debug_assert!( + self.exponent.is_positive() == rhs.exponent.is_positive(), + "checked_add overflow with mixed-sign exponents is impossible" + ); + return if self.exponent > 0 { + Repr::infinity_with_sign(sign) + } else { + Repr::zero_with_sign(sign) + }; + } + }; + match Repr::new(significand, exponent).check_finite_exponent() { + Ok(r) => r, + Err(FpError::Overflow(s)) => Repr::infinity_with_sign(s), + Err(FpError::Underflow(s)) => Repr::zero_with_sign(s), + Err(_) => unreachable!(), + } + } +} + +impl Mul<&Repr> for Repr { + type Output = Repr; + #[inline] + fn mul(self, rhs: &Repr) -> Repr { + (&self) * rhs + } +} +impl Mul> for &Repr { + type Output = Repr; + #[inline] + fn mul(self, rhs: Repr) -> Repr { + self * &rhs + } +} +impl Mul> for Repr { + type Output = Repr; + #[inline] + fn mul(self, rhs: Repr) -> Repr { + (&self) * &rhs + } +} + +#[cfg(test)] +mod tests { + use super::*; + use dashu_int::IBig; + + fn r(sig: i128, exp: isize) -> Repr { + Repr::new(IBig::from(sig), exp) + } + + #[test] + fn add_same_exponent() { + assert_eq!(&r::<10>(3, 0) + &r::<10>(4, 0), r::<10>(7, 0)); + } + + #[test] + fn add_aligns_exponents() { + // 3 + 0.4 = 3.4 + assert_eq!(&r::<10>(3, 0) + &r::<10>(4, -1), r::<10>(34, -1)); + // 1 + 0.00001 = 1.00001 (exact; the small operand's digits are all retained) + assert_eq!(&r::<10>(1, 0) + &r::<10>(1, -5), r::<10>(100001, -5)); + } + + #[test] + fn add_neg_zero_is_identity() { + let nz = Repr::<10>::neg_zero(); + let x = r::<10>(5, -2); + assert_eq!(&nz + &x, x); + assert_eq!(&x + &nz, x); + // -0 + -0 = -0: the zero short-circuit returns the other operand unchanged. + assert_eq!(&nz + &nz, Repr::<10>::neg_zero()); + } + + #[test] + fn add_cancellation_is_positive_zero() { + let z = &r::<10>(1, 0) + &r::<10>(-1, 0); + assert_eq!(z, Repr::<10>::zero()); + assert!(z.is_pos_zero()); + } + + #[test] + fn sub_basic() { + assert_eq!(&r::<10>(5, 0) - &r::<10>(3, 0), r::<10>(2, 0)); + assert_eq!(&r::<10>(3, 0) - &r::<10>(5, 0), r::<10>(-2, 0)); + // x - 0 = x (subtraction is a + (-0); the zero short-circuit keeps x's representation) + assert_eq!(&r::<10>(5, 0) - &Repr::<10>::zero(), r::<10>(5, 0)); + } + + #[test] + fn mul_basic() { + assert_eq!(&r::<10>(3, 0) * &r::<10>(4, 0), r::<10>(12, 0)); + // significands multiply, exponents add: 3e2 * 2e-1 = 6e1 + assert_eq!(&r::<10>(3, 2) * &r::<10>(2, -1), r::<10>(6, 1)); + assert_eq!(&r::<2>(3, 0) * &r::<2>(3, 0), r::<2>(9, 0)); + } + + #[test] + fn mul_zero_product_sign() { + // (+0) * (+5) = +0 + assert_eq!(&Repr::<10>::zero() * &r::<10>(5, 0), Repr::<10>::zero()); + // (+0) * (-5) = -0 (the XOR sign of the operands is attached to the zero product) + let prod = &Repr::<10>::zero() * &r::<10>(-5, 0); + assert!(prod.is_neg_zero()); + } + + #[test] + fn ref_val_combos() { + let a = r::<10>(2, 0); + let b = r::<10>(3, 0); + assert_eq!(&a + &b, r::<10>(5, 0)); + assert_eq!(a.clone() + &b, r::<10>(5, 0)); + assert_eq!(&a + b.clone(), r::<10>(5, 0)); + assert_eq!(a.clone() + b.clone(), r::<10>(5, 0)); + + let c = r::<10>(7, 0); + let d = r::<10>(4, 0); + assert_eq!(&c - &d, r::<10>(3, 0)); + assert_eq!(c.clone() - &d, r::<10>(3, 0)); + assert_eq!(&c - d.clone(), r::<10>(3, 0)); + assert_eq!(c.clone() - d.clone(), r::<10>(3, 0)); + + let e = r::<10>(6, 0); + let f = r::<10>(7, 0); + assert_eq!(&e * &f, r::<10>(42, 0)); + assert_eq!(e.clone() * &f, r::<10>(42, 0)); + assert_eq!(&e * f.clone(), r::<10>(42, 0)); + assert_eq!(e.clone() * f.clone(), r::<10>(42, 0)); + } +} diff --git a/float/src/root.rs b/float/src/root.rs index 465517f6..2c5c8e90 100644 --- a/float/src/root.rs +++ b/float/src/root.rs @@ -1,14 +1,29 @@ -use dashu_base::{Approximation, CubicRoot, Sign, SquareRoot, SquareRootRem, UnsignedAbs}; +use dashu_base::{ + Approximation, CubicRoot, EstimatedLog2, Sign, SquareRoot, SquareRootRem, UnsignedAbs, +}; use dashu_int::{IBig, UBig}; use crate::{ error::{assert_limited_precision, panic_root_zeroth, FpError, FpResult}, fbig::FBig, repr::{Context, Repr, Word}, - round::Round, + round::{ErrorBounds, Round, Rounded}, utils::{shl_digits, split_digits_ref}, }; +/// Take the value of a [`Rounded`] result, recording in `exact` whether it was computed exactly. +/// +/// Mirrors MPFR's `exact` flag: an all-exact operation chain yields the exact true value, which a +/// Ziv closure can report with radius 0 — `ziv` then accepts it without the containment test, +/// which otherwise can't certify an exactly-representable result (it sits on a one-sided preimage +/// boundary under directed rounding). +fn value_tracking_exact(r: Rounded, exact: &mut bool) -> T { + if !matches!(r, Approximation::Exact(_)) { + *exact = false; + } + r.value() +} + impl SquareRoot for FBig { type Output = Self; #[inline] @@ -256,14 +271,13 @@ impl Context { } } -impl Context { +impl Context { /// Compute `sqrt(a² + b²)` without spurious overflow/underflow. /// /// This is the overflow-safe scaled sum-of-squares: the larger-magnitude operand is never /// squared. Writing `m = max(|a|, |b|)` and `r = min(|a|,|b|) / m` (so `|r| ≤ 1`), the result is - /// `m · sqrt(1 + r²)`, where `1 + r² ∈ [1, 2]` cannot overflow. The final `m · sqrt(1 + r²)` - /// overflows only when the true result genuinely exceeds the exponent range (reported as - /// [`FpError::Overflow`]). `hypot(±inf, ·) = +inf`, `hypot(0, 0) = +0`. + /// `m · sqrt(1 + r²)`, where `1 + r² ∈ [1, 2]` cannot overflow. The result is correctly rounded + /// via a Ziv retry loop (`hypot(±inf, ·) = +inf`, `hypot(0, 0) = +0`). /// /// This is a field-arithmetic-class op (no constant cache), like `sqrt`/`atan2`. /// @@ -279,11 +293,6 @@ impl Context { return Ok(Approximation::Exact(FBig::new(Repr::zero(), *self))); } - let guard = crate::utils::ceil_usize(::log2_est( - &self.precision, - )) + 10; - let gctx = Context::::new(self.precision + guard); - // magnitudes, ordered large >= small (both finite, not both zero here) let a_mag = if a.sign() == Sign::Negative { -a.clone() @@ -302,21 +311,51 @@ impl Context { }; if small.significand.is_zero() { - // hypot(x, 0) = |x|; `large` is already a magnitude - return Ok(gctx.repr_round_ref(&large).map(|v| FBig::new(v, *self))); + // hypot(x, 0) = |x|; `large` is already a magnitude. + return Ok(self.repr_round_ref(&large).map(|v| FBig::new(v, *self))); + } + + // The result is `sqrt(large² + small²)`, i.e. ∈ [large, large·√2]. It overflows only when + // `large` is so large that the result reaches the infinity sentinel exponent — unreachable + // for real inputs, but pre-checked here so the Ziv closure can use infallible `FBig` + // arithmetic. + if large.exponent >= isize::MAX - 1 { + return Err(FpError::Overflow(Sign::Positive)); } - // r = small / large ∈ [0, 1]; 1 + r² ∈ [1, 2] (no overflow); result = large · sqrt(1+r²) - let r = gctx.div(&small, &large)?.value(); - let r2 = gctx.sqr(r.repr())?.value(); - let sum = gctx.add(&Repr::one(), r2.repr())?.value(); - let root = gctx.sqrt(sum.repr())?.value(); - let result = gctx.mul(&large, root.repr())?.value(); - Ok(result.with_precision(self.precision)) + let initial_guard = crate::utils::ceil_usize(self.precision.log2_est()) + 10; + Ok(self.ziv(initial_guard, |guard| { + let gctx = Context::::new(self.precision + guard); + // result = sqrt(large² + small²), with both operands scaled down by `k` base-B digits + // before squaring (so `large²` can't overflow the exponent) and the root scaled back: + // sqrt(L² + S²) · B^k = sqrt(large² + small²) for L = large·B⁻ᵏ, S = small·B⁻ᵏ. No + // division — so for integer inputs every step is exact (MPFR's `exact` flag), and an + // all-exact chain yields the exact true value. Report radius 0 then, which `ziv` + // accepts without the containment test (it can't certify an exactly-representable result + // under directed rounding — e.g. hypot(3,4)=5, hypot(5,12)=13 — which sits on a + // one-sided preimage boundary). + let k = (large.exponent as i128 - (isize::MAX as i128 - 2) / 2).max(0) as isize; + let mut exact = true; + let large_f = + FBig::new(value_tracking_exact(gctx.repr_round_ref(&large), &mut exact), gctx); + let small_f = + FBig::new(value_tracking_exact(gctx.repr_round_ref(&small), &mut exact), gctx); + let l_sq = value_tracking_exact(gctx.sqr((large_f >> k).repr()).unwrap(), &mut exact); + let s_sq = value_tracking_exact(gctx.sqr((small_f >> k).repr()).unwrap(), &mut exact); + let sum = value_tracking_exact(gctx.add(l_sq.repr(), s_sq.repr()).unwrap(), &mut exact); + let root = value_tracking_exact(gctx.sqrt(sum.repr()).unwrap(), &mut exact); + let result = root << k; // exact exponent shift — scales back, doesn't affect `exact` + let radius = if exact { + FBig::::ZERO + } else { + result.ulp() * 8 + }; + (result, radius) + })) } } -impl FBig { +impl FBig { /// Compute `sqrt(self² + other²)` without spurious overflow/underflow. /// /// The result precision is `max(self.precision(), other.precision())`. See @@ -376,6 +415,24 @@ mod tests { assert_eq!(r.repr().sign(), Sign::Positive); } + fn check_hypot_exact_triples(ctx: Context) { + let mk = |v: i32| Repr::<2>::new(v.into(), 0); + // Pythagorean triples: the result is exactly representable, so under a directed mode it + // sits on a one-sided preimage boundary. The closure must terminate (radius 0 from the + // all-exact `sqrt(large²+small²)` chain) rather than infinite-retry. + for (a, b, h) in [(3, 4, 5), (5, 12, 13), (8, 15, 17)] { + let r = ctx.hypot(&mk(a), &mk(b)).unwrap().value(); + assert_eq!(r.repr().significand(), &h.into(), "hypot({a}, {b})"); + } + } + + #[test] + fn test_hypot_exact_under_directed_rounding() { + check_hypot_exact_triples(Context::::new(53)); + check_hypot_exact_triples(Context::::new(53)); + check_hypot_exact_triples(Context::::new(53)); + } + #[test] fn test_hypot_no_spurious_overflow() { // a value whose square would collide with the +inf sentinel exponent, but whose diff --git a/float/src/sign.rs b/float/src/sign.rs index dd46c3fe..278cf2d1 100644 --- a/float/src/sign.rs +++ b/float/src/sign.rs @@ -60,14 +60,6 @@ impl FBig { } } -impl Neg for Repr { - type Output = Self; - #[inline] - fn neg(self) -> Self::Output { - Repr::neg(self) - } -} - impl Neg for FBig { type Output = Self; #[inline] diff --git a/float/src/third_party/num_traits.rs b/float/src/third_party/num_traits.rs index 5368fc7a..04851ae2 100644 --- a/float/src/third_party/num_traits.rs +++ b/float/src/third_party/num_traits.rs @@ -1,6 +1,9 @@ //! Implement num-traits traits. -use crate::{fbig::FBig, round::Round}; +use crate::{ + fbig::FBig, + round::{ErrorBounds, Round}, +}; use dashu_base::{Abs, DivEuclid, ParseError, RemEuclid, Sign}; use dashu_int::{IBig, Word}; use num_traits_v02 as num_traits; @@ -100,28 +103,28 @@ impl num_traits::ToPrimitive for FBig { } } -impl num_traits::Pow for FBig { +impl num_traits::Pow for FBig { type Output = FBig; fn pow(self, rhs: IBig) -> Self { self.powi(rhs) } } -impl num_traits::Pow for &FBig { +impl num_traits::Pow for &FBig { type Output = FBig; fn pow(self, rhs: IBig) -> FBig { self.powi(rhs) } } -impl num_traits::Pow<&FBig> for FBig { +impl num_traits::Pow<&FBig> for FBig { type Output = FBig; fn pow(self, rhs: &Self) -> Self { self.powf(rhs) } } -impl num_traits::Pow<&FBig> for &FBig { +impl num_traits::Pow<&FBig> for &FBig { type Output = FBig; fn pow(self, rhs: &FBig) -> FBig { diff --git a/float/src/ziv.rs b/float/src/ziv.rs new file mode 100644 index 00000000..3aea4f63 --- /dev/null +++ b/float/src/ziv.rs @@ -0,0 +1,276 @@ +//! The Ziv retry loop — guaranteed-correct rounding for transcendentals. +//! +//! Transcendentals (`exp`, `ln`, …) cannot compute an exact result the way arithmetic can, +//! so they approximate at a working precision `p + guard` and round down. A single guard-digit +//! heuristic is only *near*-correct: for an input whose true value sits on a rounding boundary, +//! the rounded result can be off by one ULP. The Ziv loop closes that gap. +//! +//! Each transcendental reports its approximation together with a provable absolute error radius +//! `E` (a true upper bound on `|approx − true|`, built from the fact that every `+`/`-`/`*`/`/` +//! in the algorithm is itself correctly rounded). The driver rounds the approximation to the +//! target precision and checks whether the entire error interval `[ã − E, ã + E]` lies inside +//! the candidate's [`ErrorBounds`] preimage — the set of reals that round to it. If it does, +//! the rounded value is *guaranteed* correct; otherwise the driver retries with more guard +//! digits. The loop provably terminates (a true tie is resolved deterministically by the mode), +//! with a large sanity cap as an unreachable backstop. + +use core::cmp::Ordering; + +use dashu_base::Approximation::*; + +use crate::{ + fbig::FBig, + repr::{Context, Repr}, + round::ErrorBounds, + round::Rounded, +}; +use dashu_int::Word; + +/// Maximum number of Ziv retries before falling back to the best-effort rounded value. +/// +/// This is a sanity backstop only — the loop converges as soon as the working precision is +/// large enough that the error interval no longer straddles a rounding boundary, which happens +/// in one attempt for essentially all inputs (the guard-digit heuristic is sized for that) and +/// in a handful of attempts only for inputs pathological close to a tie. The cap exists so a +/// bug in an error-radius bound can never produce an infinite loop. +const MAX_ZIV_RETRIES: usize = 32; + +// A test-only retry counter, so tests can assert the loop converges on the first attempt for +// typical inputs (validating that the guard-digit heuristic wasn't over-tightened). Reads as +// the number of *extra* attempts beyond the first, i.e. `0` means first-attempt success. +#[cfg(all(test, feature = "std"))] +thread_local! { + pub(crate) static LAST_ZIV_RETRIES: core::cell::Cell = const { core::cell::Cell::new(0) }; +} + +impl Context { + /// Correctly round a transcendental approximation to this context's precision using a Ziv + /// retry loop. + /// + /// `approx(guard)` computes the function at working precision `self.precision + guard` and + /// returns `(value, error_radius)` — the value and a provable upper bound on its absolute + /// error, both as [`FBig`]s at the working context. The closure is expected to capture and + /// reborrow any [`ConstCache`](crate::ConstCache) from the enclosing scope; the driver calls + /// it once per attempt and grows `guard` when the result cannot be certified. + /// + /// The loop preserves the [`Exact`](Rounded)/[`Inexact`](Rounded) flag from rounding the + /// approximation to the target precision. + pub(crate) fn ziv( + &self, + initial_guard: usize, + mut approx: impl FnMut(usize) -> (FBig, FBig), + ) -> Rounded> { + // Unlimited precision: the approximation is exact, so report it as-is. + if !self.is_limited() { + let (value, _err) = approx(0); + return Exact(value); + } + + let mut guard = initial_guard; + let mut last = None; + #[cfg(all(test, feature = "std"))] + LAST_ZIV_RETRIES.with(|c| c.set(0)); + for _ in 0..MAX_ZIV_RETRIES { + let (a, e) = approx(guard); + // `with_precision` consumes `a`, but the containment test still needs it, so round a + // clone and keep the original for the interval check. + let candidate = a.clone().with_precision(self.precision); + if Self::contained::(&a.repr, &e.repr, candidate.value_ref()) { + return candidate; + } + last = Some(candidate); + + // Grow the guard aggressively so a near-tie resolves in a couple of retries, while + // the first attempt (with the heuristic guard) handles the common case. + let step = core::cmp::max(guard, self.precision / 2).max(1); + guard += step; + #[cfg(all(test, feature = "std"))] + LAST_ZIV_RETRIES.with(|c| c.set(c.get() + 1)); + } + + // Unreachable in practice: return the best-effort candidate from the last attempt, + // matching the pre-Ziv near-correct behavior rather than looping forever. + last.expect("MAX_ZIV_RETRIES is non-zero") + } + + /// Pair variant of [`ziv`](Self::ziv) for functions that return two values (e.g. `sin_cos`, + /// `sinh_cosh`). `approx(guard)` returns `((v1, e1), (v2, e2))` — both values and their + /// provable radii at the working context, sharing whatever computation is common. The driver + /// certifies **both** values: it retries while *either* containment test fails, and returns + /// both only when both fit their rounding preimages. Shares the guard-growth loop and retry + /// counter with [`ziv`](Self::ziv). + pub(crate) fn ziv_pair( + &self, + initial_guard: usize, + mut approx: impl FnMut(usize) -> ((FBig, FBig), (FBig, FBig)), + ) -> (Rounded>, Rounded>) { + // Unlimited precision: both approximations are exact, report them as-is. + if !self.is_limited() { + let ((v1, _), (v2, _)) = approx(0); + return (Exact(v1), Exact(v2)); + } + + let mut guard = initial_guard; + let mut last = None; + #[cfg(all(test, feature = "std"))] + LAST_ZIV_RETRIES.with(|c| c.set(0)); + for _ in 0..MAX_ZIV_RETRIES { + let ((a1, e1), (a2, e2)) = approx(guard); + let c1 = a1.clone().with_precision(self.precision); + let c2 = a2.clone().with_precision(self.precision); + if Self::contained::(&a1.repr, &e1.repr, c1.value_ref()) + && Self::contained::(&a2.repr, &e2.repr, c2.value_ref()) + { + return (c1, c2); + } + last = Some((c1, c2)); + + // Grow the guard aggressively so a near-tie resolves in a couple of retries. + let step = core::cmp::max(guard, self.precision / 2).max(1); + guard += step; + #[cfg(all(test, feature = "std"))] + LAST_ZIV_RETRIES.with(|c| c.set(c.get() + 1)); + } + + // Unreachable in practice: return the best-effort pair from the last attempt. + last.expect("MAX_ZIV_RETRIES is non-zero") + } + + /// Containment test: is the approximation's error interval `[a − e, a + e]` entirely inside + /// the rounding preimage of `y` (every real in `[y − lb, y + rb]` rounds to `y` under `R`)? + /// + /// `a` and `e` are the working-precision approximation and its provable error radius; `y` is + /// the candidate rounded to the target precision (kept as an [`FBig`] only because + /// [`ErrorBounds::error_bounds`] is defined on [`FBig`]). The interval arithmetic runs on the + /// raw [`Repr`]s, which carry no precision limit, so the additions are lossless — there is no + /// rounding that could drop a guard digit and mis-decide (a wrong call here yields a wrong + /// ULP). The old path promoted every value to unlimited precision via `with_precision(0)`; + /// the [`Repr`]s are already exact, so that was a chain of no-op clones, now removed. + /// + /// The test compares sums rather than differences — algebraically identical for exact + /// arithmetic, and it reads as a single shared inequality per endpoint: + /// `a − e ≥ y − lb ⟺ a + lb ≥ y + e` + /// `a + e ≤ y + rb ⟺ y + rb ≥ a + e` + fn contained(a: &Repr, e: &Repr, y: &FBig) -> bool { + let (lb, rb, incl_l, incl_r) = R::error_bounds::(y); + + let y = &y.repr; + let lb = lb.into_repr(); + let rb = rb.into_repr(); + + let left = (a + &lb).cmp(&(y + e)); + let right = (y + &rb).cmp(&(a + e)); + let left_ok = if incl_l { + left != Ordering::Less + } else { + left == Ordering::Greater + }; + let right_ok = if incl_r { + right != Ordering::Less + } else { + right == Ordering::Greater + }; + left_ok && right_ok + } +} + +#[cfg(all(test, feature = "std"))] +mod tests { + use super::*; + use crate::round::mode; + + type F = crate::FBig; + + // An exact approximation (radius 0) is accepted on the first attempt as Exact. + #[test] + fn ziv_accepts_exact_first_attempt() { + let ctx: Context = Context::new(10); + LAST_ZIV_RETRIES.with(|c| c.set(usize::MAX)); + let r = ctx.ziv(4, |_| (F::ONE, F::ZERO)); + assert!(matches!(r, Exact(_))); + assert_eq!(LAST_ZIV_RETRIES.with(|c| c.get()), 0); + } + + // An approximation whose error interval straddles a rounding boundary must retry; once the + // shrinking radius makes the interval unambiguous, ziv accepts. + #[test] + fn ziv_retries_until_contained() { + let ctx: Context = Context::new(4); + let r = ctx.ziv(2, |guard| { + // value 1.0, radius 2^(-guard): large on the first attempt, tiny later. + (F::ONE, F::ONE >> guard as isize) + }); + let _ = r.value(); + assert!(LAST_ZIV_RETRIES.with(|c| c.get()) >= 1); + } + + // Unlimited-precision context short-circuits to a single Exact call (counter untouched). + #[test] + fn ziv_unlimited_short_circuits() { + let ctx: Context = Context::new(0); + LAST_ZIV_RETRIES.with(|c| c.set(usize::MAX)); + let r = ctx.ziv(4, |_| (F::from(7u8), F::ZERO)); + assert!(matches!(r, Exact(_))); + assert_eq!(LAST_ZIV_RETRIES.with(|c| c.get()), usize::MAX); + } + + // ziv_pair accepts an exact pair (both radii 0) on the first attempt. + #[test] + fn ziv_pair_accepts_exact_first_attempt() { + let ctx: Context = Context::new(10); + LAST_ZIV_RETRIES.with(|c| c.set(usize::MAX)); + let (r1, r2) = ctx.ziv_pair(4, |_| ((F::ONE, F::ZERO), (F::from(2u8), F::ZERO))); + assert!(matches!(r1, Exact(_))); + assert!(matches!(r2, Exact(_))); + assert_eq!(LAST_ZIV_RETRIES.with(|c| c.get()), 0); + } + + // ziv_pair retries while *either* value's interval straddles a boundary; here the second value + // carries the shrinking radius, so the pair must retry together. + #[test] + fn ziv_pair_retries_until_both_contained() { + let ctx: Context = Context::new(4); + let (r1, r2) = ctx.ziv_pair(2, |guard| { + let radius = F::ONE >> guard as isize; + ((F::ONE, F::ZERO), (F::ONE, radius)) + }); + let _ = (r1.value(), r2.value()); + assert!(LAST_ZIV_RETRIES.with(|c| c.get()) >= 1); + } + + // The guard-digit heuristic should let exp/ln converge in at most one retry for typical + // inputs. A single retry on a near-tie is by design (Ziv certifies correctness; the guard only + // controls the first-attempt hit rate). This catches gross guard mis-sizing (many retries), + // not the occasional near-tie retry. + #[test] + fn ziv_few_retries_for_typical_inputs() { + let cases = [ + F::try_from(0.5f64).unwrap(), + F::try_from(1.5f64).unwrap(), + F::try_from(2.0f64).unwrap(), + F::try_from(10.0f64).unwrap(), + F::try_from(1000.0f64).unwrap(), + F::try_from(1e-6f64).unwrap(), + ]; + const MAX_RETRIES: usize = 1; + for p in [10usize, 24, 53, 100, 200] { + for x in &cases { + let x = x.clone().with_precision(p).value(); + LAST_ZIV_RETRIES.with(|c| c.set(usize::MAX)); + let _ = x.ln(); + let ln_retries = LAST_ZIV_RETRIES.with(|c| c.get()); + assert!( + ln_retries <= MAX_RETRIES, + "ln({x}) at p={p} took {ln_retries} retries (expected <= {MAX_RETRIES})" + ); + LAST_ZIV_RETRIES.with(|c| c.set(usize::MAX)); + let _ = x.exp(); + let exp_retries = LAST_ZIV_RETRIES.with(|c| c.get()); + assert!( + exp_retries <= MAX_RETRIES, + "exp({x}) at p={p} took {exp_retries} retries (expected <= {MAX_RETRIES})" + ); + } + } + } +} diff --git a/float/tests/exp.rs b/float/tests/exp.rs index 056d20c8..2359ef59 100644 --- a/float/tests/exp.rs +++ b/float/tests/exp.rs @@ -447,9 +447,9 @@ fn test_pow_inf() { #[test] fn test_pow_with_rounding() { use dashu_base::Abs; - use dashu_float::round::{mode::*, Round}; + use dashu_float::round::{mode::*, ErrorBounds, Round}; - fn test_powf_with_error( + fn test_powf_with_error( base: &FBig, exp: &FBig, target: &FBig, diff --git a/float/tests/exp_log_root_prop.rs b/float/tests/exp_log_root_prop.rs index b2f918aa..f92f71f5 100644 --- a/float/tests/exp_log_root_prop.rs +++ b/float/tests/exp_log_root_prop.rs @@ -85,14 +85,222 @@ proptest! { prop_assert!(lower < x && x < upper); } - /// Correct-rounding self-oracle for ln: the p-precision result agrees with the - /// (rounded-down) 2p-precision result to within one ulp. + /// Correct-rounding self-oracle for ln: the p-precision result must agree *exactly* with the + /// 2p-precision result re-rounded to p (the Ziv loop guarantees both are the unique + /// correctly-rounded value). Pre-Ziv this only held to within one ulp. #[test] - fn ln_rounding_self_oracle(x in pos_x(100_000)) { + fn ln_rounding_exact_oracle(x in pos_x(100_000)) { let r_p = x.ln(); let x_2p = x.clone().with_precision(2 * P).value(); let r_2p = x_2p.ln().with_precision(P).value(); - let ulp = r_p.ulp(); - prop_assert!((r_p - r_2p).abs() <= ulp); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for exp: same construction as `ln_rounding_exact_oracle`. + #[test] + fn exp_rounding_exact_oracle(x in signed_x(2_000)) { + let r_p = x.exp(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.exp().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for exp_m1. + #[test] + fn exp_m1_rounding_exact_oracle(x in signed_x(2_000)) { + let r_p = x.exp_m1(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.exp_m1().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for ln_1p (domain x > -1; m·10⁻⁴ with m > -10000 keeps x > -1). + #[test] + fn ln_1p_rounding_exact_oracle(m in -9999i64..100_000i64) { + let x = x_from(m); + let r_p = x.ln_1p(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.ln_1p().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for hypot. + #[test] + fn hypot_rounding_exact_oracle(a in signed_x(100_000), b in signed_x(100_000)) { + let r_p = a.clone().hypot(&b); + let a_2p = a.clone().with_precision(2 * P).value(); + let b_2p = b.clone().with_precision(2 * P).value(); + let r_2p = a_2p.hypot(&b_2p).with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for sinh. + #[test] + fn sinh_rounding_exact_oracle(x in signed_x(2_000)) { + let r_p = x.sinh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.sinh().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for cosh. + #[test] + fn cosh_rounding_exact_oracle(x in signed_x(2_000)) { + let r_p = x.cosh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.cosh().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for tanh (saturates toward ±1, so a wide range is safe). + #[test] + fn tanh_rounding_exact_oracle(x in signed_x(200_000)) { + let r_p = x.tanh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.tanh().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for sinh_cosh (both components). + #[test] + fn sinh_cosh_rounding_exact_oracle(x in signed_x(2_000)) { + let (s_p, c_p) = x.sinh_cosh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let (s_2p, c_2p) = x_2p.sinh_cosh(); + prop_assert_eq!(s_p, s_2p.with_precision(P).value()); + prop_assert_eq!(c_p, c_2p.with_precision(P).value()); + } + + /// Correct-rounding self-oracle for asinh. + #[test] + fn asinh_rounding_exact_oracle(x in signed_x(200_000)) { + let r_p = x.asinh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.asinh().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for acosh (domain x ≥ 1; m·10⁻⁴ with m ≥ 10001 keeps x > 1). + #[test] + fn acosh_rounding_exact_oracle(m in 10001i64..200_000i64) { + let x = x_from(m); + let r_p = x.acosh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.acosh().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for atanh (domain |x| < 1; m·10⁻⁴ with |m| ≤ 9999). + #[test] + fn atanh_rounding_exact_oracle(m in -9999i64..9999i64) { + let x = x_from(m); + let r_p = x.atanh(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.atanh().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for sin (exercises the π/2 argument reduction for |x| > ~1.5). + #[test] + fn sin_rounding_exact_oracle(x in signed_x(200_000)) { + let r_p = x.sin(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.sin().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for cos. + #[test] + fn cos_rounding_exact_oracle(x in signed_x(200_000)) { + let r_p = x.cos(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.cos().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for sin_cos (both components). + #[test] + fn sin_cos_rounding_exact_oracle(x in signed_x(200_000)) { + let (s_p, c_p) = x.sin_cos(); + let x_2p = x.clone().with_precision(2 * P).value(); + let (s_2p, c_2p) = x_2p.sin_cos(); + prop_assert_eq!(s_p, s_2p.with_precision(P).value()); + prop_assert_eq!(c_p, c_2p.with_precision(P).value()); + } + + /// Correct-rounding self-oracle for tan (skips exact poles — unreachable for rational x, but + /// near-poles yield huge finite values that round identically). + #[test] + fn tan_rounding_exact_oracle(x in signed_x(200_000)) { + let r_p = x.tan(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.tan().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for atan. + #[test] + fn atan_rounding_exact_oracle(x in signed_x(200_000)) { + let r_p = x.atan(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.atan().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for asin (domain |x| ≤ 1; m·10⁻⁴ with |m| ≤ 10000). + #[test] + fn asin_rounding_exact_oracle(m in -10000i64..10000i64) { + let x = x_from(m); + let r_p = x.asin(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.asin().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for acos (domain |x| ≤ 1). + #[test] + fn acos_rounding_exact_oracle(m in -10000i64..10000i64) { + let x = x_from(m); + let r_p = x.acos(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = x_2p.acos().with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for atan2. + #[test] + fn atan2_rounding_exact_oracle(y in signed_x(200_000), x in signed_x(200_000)) { + let r_p = y.clone().atan2(&x); + let y_2p = y.clone().with_precision(2 * P).value(); + let x_2p = x.clone().with_precision(2 * P).value(); + let r_2p = y_2p.atan2(&x_2p).with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for powf. `base > 0` keeps it in domain; the exponent is kept + /// non-integer (m·10⁻⁴ with m % 10000 ≠ 0) so the Ziv path is exercised — integer-valued + /// exponents delegate to `powi` and are covered by `test_powf_integer_exponent`. The wide + /// exponent range (up to ±200) stresses the `exp(y·ln x)` amplification, which is exactly the + /// case the work-precision radius (vs the old unlimited-precision one) must converge on. + #[test] + fn powf_rounding_exact_oracle(b in pos_x(100_000), m in -2_000_000i64..=2_000_000i64) { + prop_assume!(m != 0 && m % 10000 != 0); + let e = x_from(m); + let r_p = b.clone().powf(&e); + let b_2p = b.clone().with_precision(2 * P).value(); + let e_2p = e.clone().with_precision(2 * P).value(); + let r_2p = b_2p.powf(&e_2p).with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); + } + + /// Correct-rounding self-oracle for powi (positive and negative integer exponents, including + /// exact integer powers where the squaring chain introduces no rounding). + #[test] + fn powi_rounding_exact_oracle(b in signed_x(50_000), n in -12i32..=12) { + prop_assume!(n != 0); + let r_p = b.powi(IBig::from(n)); + let b_2p = b.clone().with_precision(2 * P).value(); + let r_2p = b_2p.powi(IBig::from(n)).with_precision(P).value(); + prop_assert_eq!(r_p, r_2p); } } diff --git a/float/tests/log.rs b/float/tests/log.rs index 05aac914..2a8c3012 100644 --- a/float/tests/log.rs +++ b/float/tests/log.rs @@ -202,9 +202,9 @@ fn test_ln_1p_inf() { #[test] fn test_ln_with_rounding() { use dashu_base::Abs; - use dashu_float::round::{mode::*, Round}; + use dashu_float::round::{mode::*, ErrorBounds, Round}; - fn test_ln_with_error( + fn test_ln_with_error( base: &FBig, target: &FBig, atol: &FBig, diff --git a/guide-zh/src/compliance.md b/guide-zh/src/compliance.md index b6582be1..7e3fecfd 100644 --- a/guide-zh/src/compliance.md +++ b/guide-zh/src/compliance.md @@ -51,7 +51,7 @@ | IEEE 754 要求 | 合规性 | 说明 | |---------------------|-----------|-------| | 舍入模式：roundTiesToEven、roundTiesToAway、roundTowardPositive、roundTowardNegative、roundTowardZero | ✅ | 全部五种模式实现为 `HalfEven`、`HalfAway`、`Up`、`Down`、`Zero`。 | -| 正确舍入到 1 ulp 以内 | ✅ | 所有运算保证 $|error| < 1\text{ ulp}$。`Rounded` 类型区分精确与非精确结果。 | +| 正确舍入到 1 ulp 以内 | ✅ | 所有运算保证 $|error| < 1\text{ ulp}$。算术、根号以及实超越函数（`exp`、`exp_m1`、`ln`、`ln_1p`、三角与双曲函数族、`hypot`、`powi`、`powf`）通过 Ziv 循环保证正确舍入；`dashu-cmplx` 的复数超越函数（`exp`、`log`、`powf`、`powi`、`sin`/`cos`/`tan`/`sin_cos`、`asin`/`acos`/`atan`、`sqrt`）经由其自身的 Ziv 循环正确舍入，同时认证实部与虚部（`asin`/`acos` 采用因式分解的 $1-z^2=(1-z)(1+z)$，在奇点 $z = \pm 1$ 处亦保持精度）。`Rounded` 类型区分精确与非精确结果。 | | 就近舍入保持零的符号 | ✅ | `rounded_to_repr` 在舍入将非零值折叠为零时保持输入符号。 | #### 第 5.6 节——符号位运算 @@ -143,7 +143,7 @@ |---------------------|-----------|-------| | 无效/不定形式 → NaN | ❌ 偏离 | 在 `Context` 层返回 `Err(FpError::{Indeterminate, InfiniteInput})`；在便捷层 panic。设计上无 NaN。 | | 定义域错误（如负数的偶数根、超出范围的逆三角函数） | ❌ 偏离 | 返回 `Err(FpError::OutOfDomain)` / panic，而非 NaN 结果。 | -| 每个分量独立按共享模式舍入 | ✅ | 每个轴近似正确舍入（保证正确的 Ziv 循环推迟到 0.5.x）。 | +| 每个分量独立按共享模式舍入 | ✅ | 每个轴经由同时认证两部分的 Ziv 循环正确舍入，与 `dashu-float` 的实超越函数一致（`asin`/`acos` 采用因式分解的 $1-z^2=(1-z)(1+z)$，在奇点 $z = \pm 1$ 处亦保持精度）。 | ### 小结（dashu-cmplx） diff --git a/guide-zh/src/faq.md b/guide-zh/src/faq.md index 3823e9f6..91e34745 100644 --- a/guide-zh/src/faq.md +++ b/guide-zh/src/faq.md @@ -26,7 +26,7 @@ ## 已知限制 - **没有 NaN。** 无效运算在便捷层会 panic，在上下文层返回 `Err(FpError)`。无穷是终态值，不是操作数——参见[标准合规性](./compliance.md)。 -- **近似正确舍入。** 超越函数通过保护位策略在 1 ulp 以内舍入；保证正确的 Ziv 循环计划在后续版本中实现。 +- **正确舍入的超越函数。** 实超越函数——`exp`、`exp_m1`、`ln`、`ln_1p`、三角与双曲函数族、`hypot`、`powi`、`powf`——均通过 Ziv 重试循环保证正确舍入。`dashu-cmplx` 的复数超越函数（`exp`、`log`、`powf`、`powi`、`sin`/`cos`/`tan`/`sin_cos`、`asin`/`acos`/`atan`、`sqrt`）经由其自身的 Ziv 循环正确舍入，同时认证实部与虚部；`tan` 采用免消去的双角公式（对所有有限 `|Im z|` 精确），`asin`/`acos` 采用因式分解的 $1-z^2=(1-z)(1+z)$（在 $z = \pm 1$ 处亦精确）。 - **复数功能覆盖。** `CBig` 提供域运算和基本超越函数；复数双曲函数、`fma` 及其他若干功能推迟到 0.5.x（参见 v0.5 版本说明）。 - **没有球体/区间运算。** 与 MPC 的实验性 `mpcb_t`（复数球体）不同，`dashu` 不提供区间或球体类型。这是一个刻意的范围选择：球体运算在上游仍是实验性的，如果/当它稳定后，最好由一个**独立的 crate** 在 `dashu-float`/`dashu-cmplx` 之上提供，而非耦合到核心类型。 - **尚无 SIMD-FFT 乘法**（计划在 v1.0 中实现）。 diff --git a/guide/src/compliance.md b/guide/src/compliance.md index dde61f8a..3ca72e48 100644 --- a/guide/src/compliance.md +++ b/guide/src/compliance.md @@ -56,7 +56,7 @@ The reference here is IEEE Std 754™-2008 (ISO/IEC/IEEE 60559:2011). | IEEE 754 requirement | Compliance | Notes | |---------------------|-----------|-------| | Rounding modes: roundTiesToEven, roundTiesToAway, roundTowardPositive, roundTowardNegative, roundTowardZero | ✅ | All five modes implemented as `HalfEven`, `HalfAway`, `Up`, `Down`, `Zero`. | -| Correct rounding to within 1 ulp | ✅ | All operations guarantee $|error| < 1\text{ ulp}$. The `Rounded` type distinguishes exact from inexact results. | +| Correct rounding to within 1 ulp | ✅ | All operations guarantee $|error| < 1\text{ ulp}$. Arithmetic, roots, and the real transcendentals (`exp`, `exp_m1`, `ln`, `ln_1p`, the trigonometric and hyperbolic families, `hypot`, `powi`, and `powf`) are guaranteed-correctly rounded via a Ziv loop; and `dashu-cmplx`'s complex transcendentals (`exp`, `log`, `powf`, `powi`, `sin`/`cos`/`tan`/`sin_cos`, `asin`/`acos`/`atan`, `sqrt`) are correctly rounded via their own Ziv loop that certifies both parts (`asin`/`acos` use the factored $1-z^2=(1-z)(1+z)$ so they stay accurate right up to the singularities $z = \pm 1$). The `Rounded` type distinguishes exact from inexact results. | | Round-to-nearest preserves sign of zero | ✅ | `rounded_to_repr` preserves input sign when rounding collapses a non-zero to zero. | #### Section 5.6 — Sign bit operations @@ -152,7 +152,7 @@ the `Context` layer (and panics at the convenience layer). |---------------------|-----------|-------| | Invalid / indeterminate form → NaN | ❌ Deviates | `Err(FpError::{Indeterminate, InfiniteInput})` at `Context`; panics at the convenience layer. No NaN by design. | | Domain error (e.g. even root of a negative value, out-of-range inverse trig) | ❌ Deviates | `Err(FpError::OutOfDomain)` / panic, rather than a NaN result. | -| Each component rounded independently to the shared mode | ✅ | Near-correctly rounded per axis (a guaranteed-correct Ziv loop is deferred to 0.5.x). | +| Each component rounded independently to the shared mode | ✅ | Correctly rounded per axis via a Ziv loop that certifies both parts, mirroring `dashu-float`'s real transcendentals (`asin`/`acos` use the factored $1-z^2=(1-z)(1+z)$ so they stay accurate right up to the singularities $z = \pm 1$). | ### Summary (dashu-cmplx) diff --git a/guide/src/faq.md b/guide/src/faq.md index 22e8028f..c3feaf11 100644 --- a/guide/src/faq.md +++ b/guide/src/faq.md @@ -29,7 +29,7 @@ Compared with other Rust crates: ## Known limitations - **No NaN.** Invalid operations panic at the convenience layer and return `Err(FpError)` at the context layer. Infinities are terminal values, not operands — see [Standards Compliance](./compliance.md). -- **Near-correct rounding.** Transcendentals are rounded within 1 ulp via a guard-digit recipe; a guaranteed-correct Ziv loop is planned for a later release. +- **Correctly-rounded transcendentals.** The real transcendentals — `exp`, `exp_m1`, `ln`, `ln_1p`, the trigonometric and hyperbolic families, `hypot`, `powi`, and `powf` — are guaranteed-correctly rounded via a Ziv retry loop. `dashu-cmplx`'s complex transcendentals (`exp`, `log`, `powf`, `powi`, `sin`/`cos`/`tan`/`sin_cos`, `asin`/`acos`/`atan`, `sqrt`) are correctly rounded via their own Ziv loop that certifies both parts; `tan` uses the cancellation-free double-angle form (accurate for all finite `|Im z|`), and `asin`/`acos` use the factored $1-z^2=(1-z)(1+z)$ (accurate right up to $z = \pm 1$). - **Complex surface.** `CBig` ships field arithmetic and the elementary transcendentals; complex hyperbolics, `fma`, and several others are deferred to 0.5.x (see the v0.5 release notes). - **No ball/interval arithmetic.** Unlike MPC's experimental `mpcb_t` (complex balls), `dashu` does not provide interval or ball types. This is a deliberate scope choice: ball arithmetic is diff --git a/integer/CHANGELOG.md b/integer/CHANGELOG.md index 63053385..2285cbe3 100644 --- a/integer/CHANGELOG.md +++ b/integer/CHANGELOG.md @@ -3,6 +3,13 @@ ## Unreleased ### Fix +- NTT squaring/multiplication of an all-zero operand no longer panics. `add_signed_sqr_conv` and + `add_signed_mul_conv` assumed a non-zero input (`debug_assert!(la_bits > 0)`), so an all-zero slice + panicked in debug and indexed out of bounds in release. Reachable from `sqrt_rem` of a perfect + square with many trailing zero words — it squares the all-zero low half of the estimate to verify + the remainder — which was the crash behind `dashu-float`'s `hypot(3,4)` under directed rounding. + An all-zero operand now returns early (the product is zero, so the signed accumulate is a no-op), + matching the existing zero-guard in the chunked-multiply closure. - `to_f64` of a magnitude at `DoubleWord::MAX` (e.g. `IBig::from(-(2^128 - 1)).to_f64()` on 64-bit targets) reported the inexact conversion as `Approximation::Exact`. The exactness test used a saturating `f as DoubleWord` round-trip, which clamps back to `DoubleWord::MAX` when the value diff --git a/integer/src/mul/ntt/mod.rs b/integer/src/mul/ntt/mod.rs index fa51f428..2e5e569f 100644 --- a/integer/src/mul/ntt/mod.rs +++ b/integer/src/mul/ntt/mod.rs @@ -319,10 +319,14 @@ fn add_signed_mul_conv( let lb = b.len(); debug_assert!(la > 0 && lb > 0); - let (b_pack, nn, k_eff) = select_params(la, lb); let la_bits = bit_len(a); let lb_bits = bit_len(b); - debug_assert!(la_bits > 0 && lb_bits > 0); + if la_bits == 0 || lb_bits == 0 { + // Either operand is all-zero, so `a·b = 0` and `c += sign·a·b` is a no-op (carry 0, `c` + // unchanged). Reachable from `div_rem` multiplying a zero partial quotient/operand. + return 0; + } + let (b_pack, nn, k_eff) = select_params(la, lb); let coeffs_a = coeff_count(la_bits, b_pack); let coeffs_b = coeff_count(lb_bits, b_pack); diff --git a/integer/src/sqr/ntt.rs b/integer/src/sqr/ntt.rs index c612e7da..2bd983b7 100644 --- a/integer/src/sqr/ntt.rs +++ b/integer/src/sqr/ntt.rs @@ -33,9 +33,14 @@ pub(crate) fn add_signed_sqr_same_len( fn add_signed_sqr_conv(c: &mut [Word], sign: Sign, a: &[Word], memory: &mut Memory) -> SignedWord { let la = a.len(); debug_assert!(la > 0); - let (b_pack, nn, k_eff) = ntt::select_params(la, la); let la_bits = bit_len(a); - debug_assert!(la_bits > 0); + if la_bits == 0 { + // `a` is all-zero, so `a² = 0` and `c += sign·a²` is a no-op (carry 0, `c` unchanged). + // Reachable from `sqrt_rem_large` squaring the all-zero low half of a sqrt estimate for + // a perfect square with trailing zero words (e.g. `sqrt(25·2^k)`). + return 0; + } + let (b_pack, nn, k_eff) = ntt::select_params(la, la); let coeffs_a = coeff_count(la_bits, b_pack); let output_coeffs = 2 * coeffs_a - 1; diff --git a/integer/tests/root.rs b/integer/tests/root.rs index 3b165bc0..fb5fea0a 100644 --- a/integer/tests/root.rs +++ b/integer/tests/root.rs @@ -90,6 +90,19 @@ fn test_sqrt_negative_panic() { let _ = ibig!(-1).sqrt(); } +#[test] +fn test_sqrt_rem_perfect_square_trailing_zeros() { + // Regression: a perfect square whose estimate has many trailing zero bits. `sqrt_rem` squares + // the estimate's all-zero low half to verify the remainder, routing through the NTT sqr path + // which previously panicked on the zero input (`la_bits == 0`). The shift is large enough that + // the low half exceeds THRESHOLD_NTT (4000 words), forcing the NTT path. This is the + // `dashu-float` `hypot(3,4)` crash, reached directly through the integer sqrt. + let shift = 560_000usize; + let s = ubig!(1) << shift; + let n = &s * &s; + assert_eq!(n.sqrt_rem(), (s, ubig!(0))); +} + #[test] fn test_nth_root() { // the nth root of 0 is 0 (regression: it used to return 1 via the `bits <= n` shortcut)