diff --git a/.gitignore b/.gitignore index 861f28ef..cb6a98f3 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ /.vscode/settings.json benchmark/Cargo.lock benchmark/target +fuzz/target +fuzz/Cargo.lock \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml index cb56f25f..2f7761ca 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -26,7 +26,7 @@ members = [ "python", "rational", ] -exclude = ["benchmark"] +exclude = ["benchmark", "fuzz"] default-members = ["base", "integer", "float", "rational", "macros"] [features] diff --git a/float/CHANGELOG.md b/float/CHANGELOG.md index 6241ea4e..33f831e6 100644 --- a/float/CHANGELOG.md +++ b/float/CHANGELOG.md @@ -3,6 +3,7 @@ ## Unreleased - Fix rounding issues in `to_32()` and `to_f64()` (fixes [#53](https://github.com/cmpute/dashu/issues/53) and [#56](https://github.com/cmpute/dashu/issues/56)). +- Fix `FBig::fract()` inflating context precision for values smaller than one. ## 0.4.4 diff --git a/float/src/error.rs b/float/src/error.rs index 50ae8b15..33df1fcb 100644 --- a/float/src/error.rs +++ b/float/src/error.rs @@ -37,6 +37,26 @@ pub const fn panic_power_negative_base() -> ! { } /// Panics when taking an even order root of an negative number -pub(crate) fn panic_root_negative() -> ! { +pub fn panic_root_negative() -> ! { panic!("the root is a complex number!") } + +/// Panics when the result of an operation is NaN +pub fn panic_nan() -> ! { + panic!("the result of the operation is NaN!") +} + +/// Panics when the result of an operation overflows +pub fn panic_overflow() -> ! { + panic!("the result of the operation overflowed!") +} + +/// Panics when the result of an operation underflows +pub fn panic_underflow() -> ! { + panic!("the result of the operation underflowed!") +} + +/// Panics when the result of an operation is an exact infinity +pub fn panic_infinite() -> ! { + panic!("the result of the operation is an exact infinity!") +} diff --git a/float/src/lib.rs b/float/src/lib.rs index 1f49a8fe..1af2a3d1 100644 --- a/float/src/lib.rs +++ b/float/src/lib.rs @@ -76,6 +76,7 @@ mod fmt; mod helper_macros; mod iter; mod log; +pub mod math; mod mul; pub mod ops; mod parse; diff --git a/float/src/math/consts.rs b/float/src/math/consts.rs new file mode 100644 index 00000000..004ff0c4 --- /dev/null +++ b/float/src/math/consts.rs @@ -0,0 +1,105 @@ +use crate::{ + error::assert_limited_precision, + fbig::FBig, + repr::{Context, Word}, + round::{Round, Rounded}, +}; +use dashu_base::{BitTest, Sign, UnsignedAbs}; +use dashu_int::{IBig, UBig}; + +impl Context { + /// Calculate π using the Chudnovsky algorithm with binary splitting. + /// + /// The Chudnovsky algorithm is one of the most efficient methods for + /// high-precision π calculation, providing ~14.18 decimal digits per term. + /// + /// # Methodology + /// We use Binary Splitting to evaluate the series. This technique transforms + /// the linear-time summation into a recursive tree evaluation. By combining + /// terms into large products, it allows the library to leverage fast + /// multiplication algorithms (like Toom-3 or FFT) as the numbers grow, + /// leading to significant performance gains over simple iterative summation. + /// + /// // TODO: consider adding a static cache for π at common precisions. + #[must_use] + pub fn pi(&self) -> Rounded> { + assert_limited_precision(self.precision); + + // Calculate required bits based on target precision in base B. + // bits = ceil(precision * log2(B)) + let bits = if B.is_power_of_two() { + self.precision.saturating_mul(B.ilog2() as usize) + } else { + self.precision.saturating_mul(B.ilog2() as usize + 1) + }; + + let num_terms = (bits * 100 / 4708) + 1; + let guard_bits = num_terms.bit_len() + 32; + let work_bits = bits + guard_bits; + + // Evaluate the series components using binary splitting + let (_p, q, t) = chudnovsky_bs(0, num_terms); + + // Final formula: pi = (426880 * sqrt(10005) * Q) / T + + // Convert work bits back to base B precision. + // precision_B = ceil(work_bits / log2(B)) + let work_precision = if B == 2 { + work_bits + } else { + work_bits / B.ilog2() as usize + 1 + }; + let work_context = Self::new(work_precision); + + let q_f = work_context.convert_int::(q.into()).value(); + let t_f = work_context.convert_int::(t).value(); + + let sqrt_10005 = work_context + .sqrt(&work_context.convert_int::(10005.into()).value().repr) + .value(); + let constant = work_context.convert_int::(426_880.into()).value(); + + let pi = (constant * sqrt_10005 * q_f) / t_f; + pi.with_precision(self.precision) + } +} + +/// Binary splitting implementation for the Chudnovsky series. +/// Returns (P, Q, T) for the range [a, b). +fn chudnovsky_bs(a: usize, b: usize) -> (UBig, UBig, IBig) { + if b - a == 1 { + // Base case: calculate single term + if a == 0 { + return (UBig::ONE, UBig::ONE, IBig::from_parts_const(Sign::Positive, 13_591_409)); + } + + let k = a as u64; + let p = UBig::from(6 * k - 5) * (2 * k - 1) * (6 * k - 1); + let q = UBig::from(k).pow(3) * UBig::from(10_939_058_860_032_000_u64); + let t_val = IBig::from_parts_const(Sign::Positive, 13_591_409) + + IBig::from_parts_const(Sign::Positive, 545_140_134) * k; + let t_abs = &p * t_val.unsigned_abs(); + let t = IBig::from(t_abs) * Sign::from(a % 2 == 1); + return (p, q, t); + } + + // Recursive step + let mid = (a + b) / 2; + let (p_l, q_l, t_l) = chudnovsky_bs(a, mid); + let (p_r, q_r, t_r) = chudnovsky_bs(mid, b); + + let p = &p_l * &p_r; + let q = &q_l * &q_r; + // T = T_L * Q_R + T_R * P_L + let t = IBig::from(q_r) * t_l + IBig::from(p_l) * t_r; + (p, q, t) +} + +impl FBig { + /// Calculate π with the given precision and the default rounding mode. + #[inline] + #[must_use] + pub fn pi(precision: usize) -> Self { + Context::::new(precision).pi().value() + } +} diff --git a/float/src/math/mod.rs b/float/src/math/mod.rs index c642ffba..1c0e4313 100644 --- a/float/src/math/mod.rs +++ b/float/src/math/mod.rs @@ -1,31 +1,86 @@ -//! Implementations of advanced math functions +//! Advanced mathematical functions -// TODO: implement the math functions as associated methods, and add them to FBig through a trait -// REF: https://pkg.go.dev/github.com/ericlagergren/decimal +use crate::{ + error::{panic_infinite, panic_nan, panic_overflow, panic_underflow}, + fbig::FBig, + repr::{Context, Repr, Word}, + round::{Round, Rounded}, +}; -enum FpResult { - Normal(Repr), +pub mod consts; +pub mod trig; + +/// The result of an advanced mathematical operation. +/// +/// This enum is used to handle non-finite results (NaN, Infinite) and +/// boundary conditions (Overflow, Underflow) without panicking, +/// as the core [`FBig`] type only represents finite numbers. +/// +/// Finite results are wrapped in a [Rounded] to preserve rounding information. +#[derive(Clone, Debug, PartialEq, Eq)] +pub enum FpResult { + Normal(Rounded>), Overflow, Underflow, NaN, - /// An exact infinite result is obtained from finite inputs, such as - /// divide by zero, logarithm on zero. + /// divide by zero or logarithm of zero. Infinite, } -impl Context { - fn sin(&self, repr: Repr) -> FpResult { - todo!() +impl FpResult { + /// Convert the result into an [`FBig`] with the given context. + /// + /// # Panics + /// Panics if the result is not `Normal`. + #[inline] + #[must_use] + pub fn value(self, context: &Context) -> FBig { + match self { + Self::Normal(rounded) => FBig::new(rounded.value(), *context), + Self::NaN => panic_nan(), + Self::Infinite => panic_infinite(), + Self::Overflow => panic_overflow(), + Self::Underflow => panic_underflow(), + } + } + + /// Convert the result into an optional [`FBig`] with the given context. + /// Returns `None` if the result is not `Normal`. + #[inline] + #[must_use] + pub fn ok(self, context: &Context) -> Option>> { + match self { + Self::Normal(rounded) => Some(rounded.map(|repr| FBig::new(repr, *context))), + _ => None, + } + } + + /// Returns `true` if the result is `NaN`. + #[inline] + #[must_use] + pub const fn is_nan(&self) -> bool { + matches!(self, Self::NaN) + } + + /// Returns `true` if the result is `Infinite`. + #[inline] + #[must_use] + pub const fn is_infinite(&self) -> bool { + matches!(self, Self::Infinite) } -} -trait ContextOps { - fn context(&self) -> &Context; - fn repr(&self) -> &Repr; + /// Returns `true` if the result is a normal finite value. + #[inline] + #[must_use] + pub const fn is_normal(&self) -> bool { + matches!(self, Self::Normal(_)) + } + /// Returns `true` if the result is a finite value (Normal, Overflow, or Underflow). #[inline] - fn sin(&self) -> FpResult { - self.context().sin(self.repr()) + #[must_use] + pub const fn is_finite(&self) -> bool { + matches!(self, Self::Normal(_) | Self::Overflow | Self::Underflow) } -} \ No newline at end of file +} diff --git a/float/src/math/trig.rs b/float/src/math/trig.rs new file mode 100644 index 00000000..457e8cda --- /dev/null +++ b/float/src/math/trig.rs @@ -0,0 +1,587 @@ +use crate::{ + error::assert_limited_precision, + fbig::FBig, + math::FpResult, + repr::{Context, Repr, Word}, + round::Round, +}; +use core::cmp::Ordering; +use core::convert::TryFrom; +use dashu_base::{AbsOrd, RemEuclid, Sign}; +use dashu_int::IBig; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +enum Quadrant { + First, + Second, + Third, + Fourth, +} + +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 { + // 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 work_precision = self + .precision + .saturating_add(x_mag) + .saturating_add(extra_guards); + 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`. + fn reduce_to_quadrant(self, x: &Repr) -> (Self, FBig, Quadrant) { + let work_context = self.compute_work_context_trig(x); + let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + + let pi = work_context.pi::().value(); + let half_pi = &pi / 2; + let x_scaled: FBig = &x_f / &half_pi; + let k_f = x_scaled.round(); + let r = x_f - &k_f * half_pi; + let Ok(k) = IBig::try_from(k_f) else { + unreachable!( + "round() always returns an integer and trig functions ensure input is finite" + ); + }; + + 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]"); + }; + let quadrant = match k_mod_4_int { + 0 => Quadrant::First, + 1 => Quadrant::Second, + 2 => Quadrant::Third, + 3 => Quadrant::Fourth, + _ => unreachable!(), + }; + + (work_context, r, quadrant) + } + + /// Calculate the sine of the floating point representation. + #[must_use] + pub fn sin(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + assert_limited_precision(self.precision); + + if x.is_zero() { + let res = FBig::::ZERO.with_precision(self.precision); + return FpResult::Normal(res.map(|v| v.repr)); + } + + let (work_context, r, quadrant) = self.reduce_to_quadrant(x); + + // 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), + }; + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + + /// Internal Taylor series for sine: S(x) = x - x^3/3! + x^5/5! - ... + fn sin_internal(self, x: &FBig) -> FBig { + if x.repr.significand.is_zero() { + return FBig::ZERO; + } + let x2 = x.sqr(); + let mut sum = x.clone(); + let mut term = x.clone(); + let mut k = 1usize; + let threshold = sum.sub_ulp(); + loop { + term *= &x2; + term /= (2 * k) * (2 * k + 1); + if term.abs_cmp(&threshold).is_le() { + break; + } + if k % 2 == 1 { + sum -= &term; + } else { + sum += &term; + } + k += 1; + } + sum + } + + /// Calculate the cosine of the floating point representation. + #[must_use] + pub fn cos(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + assert_limited_precision(self.precision); + + if x.is_zero() { + let res = FBig::::ONE.with_precision(self.precision); + return FpResult::Normal(res.map(|v| v.repr)); + } + + let (work_context, r, quadrant) = self.reduce_to_quadrant(x); + + // 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), + }; + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + + /// Internal Taylor series for cosine: C(x) = 1 - x^2/2! + x^4/4! - ... + fn cos_internal(self, x: &FBig) -> FBig { + if x.repr.significand.is_zero() { + return FBig::ONE.with_precision(self.precision).value(); + } + let x2 = x.sqr(); + let mut sum = FBig::::ONE.with_precision(self.precision).value(); + let mut term = sum.clone(); + let mut k = 1usize; + let threshold = sum.sub_ulp(); + loop { + term *= &x2; + term /= (2 * k) * (2 * k - 1); + if term.abs_cmp(&threshold).is_le() { + break; + } + if k % 2 == 1 { + sum -= &term; + } else { + sum += &term; + } + k += 1; + } + sum + } + + /// Calculate both the sine and cosine of the floating point representation. + /// + /// This is more efficient than calling `sin` and `cos` separately. + #[must_use] + pub fn sin_cos(&self, x: &Repr) -> (FpResult, FpResult) { + if x.is_infinite() { + return (FpResult::NaN, FpResult::NaN); + } + assert_limited_precision(self.precision); + + if x.is_zero() { + let s = FBig::::ZERO.with_precision(self.precision); + let c = FBig::::ONE.with_precision(self.precision); + return (FpResult::Normal(s.map(|v| v.repr)), FpResult::Normal(c.map(|v| v.repr))); + } + + let (work_context, r, quadrant) = self.reduce_to_quadrant(x); + + 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), + }; + + ( + FpResult::Normal(s.with_precision(self.precision).map(|v| v.repr)), + FpResult::Normal(c.with_precision(self.precision).map(|v| v.repr)), + ) + } + + /// Simultaneously evaluate Taylor series for sine and cosine. + pub(crate) fn sin_cos_internal( + self, + x: &FBig, + ) -> (FBig, FBig) { + if x.repr.significand.is_zero() { + return (FBig::ZERO, FBig::ONE.with_precision(self.precision).value()); + } + let x2 = x.sqr(); + let mut sin_sum = x.clone(); + let mut cos_sum = FBig::::ONE.with_precision(self.precision).value(); + let mut sin_term = x.clone(); + let mut cos_term = cos_sum.clone(); + let mut k = 1usize; + let sin_threshold = sin_sum.sub_ulp(); + let cos_threshold = cos_sum.sub_ulp(); + loop { + cos_term *= &x2; + cos_term /= (2 * k) * (2 * k - 1); + sin_term *= &x2; + sin_term /= (2 * k) * (2 * k + 1); + + if sin_term.abs_cmp(&sin_threshold).is_le() && cos_term.abs_cmp(&cos_threshold).is_le() + { + break; + } + + if k % 2 == 1 { + cos_sum -= &cos_term; + sin_sum -= &sin_term; + } else { + cos_sum += &cos_term; + sin_sum += &sin_term; + } + k += 1; + } + (sin_sum, cos_sum) + } + + /// Calculate the tangent of the floating point representation. + /// + /// # Note + /// Near odd multiples of π/2, the result returns `Infinite`. + #[must_use] + pub fn tan(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + assert_limited_precision(self.precision); + + if x.is_zero() { + let res = FBig::::ZERO.with_precision(self.precision); + return FpResult::Normal(res.map(|v| v.repr)); + } + + let (work_context, r, quadrant) = self.reduce_to_quadrant(x); + 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_zero() { + return FpResult::Infinite; + } + FpResult::Normal(self.div(&s_f.repr, &c_f.repr).map(|v| v.repr)) + } + + /// Calculate the arcsine of the floating point representation. + /// + /// # Methodology + /// Uses the identity: `asin(x) = atan(x / sqrt(1 - x^2))` + /// Returns `NaN` if `|x| > 1`. + #[must_use] + pub fn asin(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + 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() { + return FpResult::NaN; + } + + 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); + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + + fn asin_internal(self, x_f: &FBig) -> FBig { + let one = FBig::::ONE.with_precision(self.precision).value(); + let x2 = x_f.sqr(); + let d = self.sqrt(&(one - x2).repr).value(); + + if d.repr.is_zero() { + let pi = self.pi::().value(); + let half_pi: FBig = pi / 2; + if x_f.sign() == Sign::Positive { + return half_pi; + } + return -half_pi; + } + + self.atan_with_reduction(&(x_f / d)) + } + + /// Calculate the arccosine of the floating point representation. + /// + /// # Methodology + /// Uses the identity: `acos(x) = pi/2 - asin(x)`. + /// Higher precision is used internally to avoid catastrophic cancellation near x ≈ 1. + #[must_use] + pub fn acos(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + 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() { + return FpResult::NaN; + } + + 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); + let pi = work_context.pi::().value(); + let half_pi: FBig = pi / 2; + let res: FBig = half_pi - asin_x; + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + + /// Calculate the arctangent of the floating point representation. + #[must_use] + pub fn atan(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + let pi = self.pi::().value(); + let half_pi: FBig = pi / 2; + let res: FBig = if x.sign() == Sign::Positive { + half_pi + } else { + -half_pi + }; + return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)); + } + + assert_limited_precision(self.precision); + + if x.is_zero() { + let res = FBig::::ZERO.with_precision(self.precision); + return FpResult::Normal(res.map(|v| v.repr)); + } + + 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); + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + + /// Internal arctangent that includes range reduction but no guard digit allocation. + fn atan_with_reduction(self, x_f: &FBig) -> 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::().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 + } + + /// Internal series for arctangent. + /// Evaluates the Euler series for arctangent. + fn atan_internal(self, x: &FBig) -> FBig { + // Euler's series for atan(x) + let x2 = x.sqr(); + let one_plus_x2 = FBig::ONE + &x2; + let mut term = x / &one_plus_x2; + let mut sum = term.clone(); + let factor = (2 * &x2) / one_plus_x2; + let mut n = 1usize; + let threshold = sum.sub_ulp(); + loop { + term *= &factor; + term *= n; + term /= 2 * n + 1; + if term.abs_cmp(&threshold).is_le() { + break; + } + sum += &term; + n += 1; + } + sum + } + + /// Calculate the arctangent of y / x. + /// + /// Handles signed infinities according to IEEE 754 standards. + /// Returns `NaN` if both arguments are zero. + #[must_use] + pub fn atan2(&self, y: &Repr, x: &Repr) -> FpResult { + if y.is_zero() && x.is_zero() { + return FpResult::NaN; + } + + 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 + if y.is_infinite() || x.is_infinite() { + let (sy, sx) = (y.sign() == Sign::Positive, x.sign() == Sign::Positive); + let res: FBig = match (y.is_infinite(), x.is_infinite(), sy, sx) { + (true, true, true, true) => work_context.pi::().value() / 4, + (true, true, true, false) => work_context.pi::().value() * 3 / 4, + (true, true, false, true) => { + let pi4: FBig = work_context.pi::().value() / 4; + -pi4 + } + (true, true, false, false) => { + let pi34: FBig = work_context.pi::().value() * 3 / 4; + -pi34 + } + (true, false, true, _) => work_context.pi::().value() / 2, + (true, false, false, _) => { + let half_pi: FBig = work_context.pi::().value() / 2; + -half_pi + } + (false, true, _, true) => FBig::::ZERO.with_precision(work_precision).value(), + (false, true, true, false) => work_context.pi::().value(), + (false, true, false, false) => -work_context.pi::().value(), + _ => unreachable!(), + }; + // Note: atan2(finite, +inf) returns unsigned ZERO. IEEE 754 requires signed zero, + // but `Repr` does not distinguish signed zero. + return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)); + } + + 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); + + match x_f.cmp(&FBig::::ZERO) { + Ordering::Greater => { + let res = work_context.atan_with_reduction(&(y_f / x_f)); + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + Ordering::Less => { + let pi = work_context.pi::().value(); + let y_sign = y_f.sign(); + let atan_yx = work_context.atan_with_reduction(&(y_f / x_f)); + let res = if y_sign == Sign::Positive { + atan_yx + pi + } else { + atan_yx - pi + }; + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + Ordering::Equal => { + // x == 0 case + let pi = work_context.pi::().value(); + let half_pi: FBig = pi / 2; + if y_f > FBig::::ZERO { + FpResult::Normal(half_pi.with_precision(self.precision).map(|v| v.repr)) + } else { + let res = -half_pi; + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) + } + } + } + } +} + +impl FBig { + /// Calculate the sine of the floating point number. + /// + /// # Panics + /// Panics if the input is infinite or the result is not representable as a normal value. + #[inline] + #[must_use] + pub fn sin(&self) -> Self { + self.context.sin(&self.repr).value(&self.context) + } + + /// Calculate the cosine of the floating point number. + /// + /// # Panics + /// Panics if the input is infinite or the result is not representable as a normal value. + #[inline] + #[must_use] + pub fn cos(&self) -> Self { + self.context.cos(&self.repr).value(&self.context) + } + + /// Calculate both the sine and cosine of the floating point number. + /// + /// This is more efficient than calling `sin` and `cos` separately. + /// + /// # Panics + /// Panics if the input is infinite or the results are not representable as normal values. + #[inline] + #[must_use] + pub fn sin_cos(&self) -> (Self, Self) { + let (s, c) = self.context.sin_cos(&self.repr); + (s.value(&self.context), c.value(&self.context)) + } + + /// Calculate the tangent of the floating point number. + /// + /// Returns `FpResult` to safely handle non-finite results (e.g., at singularities). + #[inline] + #[must_use] + pub fn tan(&self) -> FpResult { + self.context.tan(&self.repr) + } + + /// Calculate the arcsine of the floating point number. + /// + /// Returns `FpResult` to safely handle domain errors (e.g., |x| > 1). + #[inline] + #[must_use] + pub fn asin(&self) -> FpResult { + self.context.asin(&self.repr) + } + + /// Calculate the arccosine of the floating point number. + /// + /// Returns `FpResult` to safely handle domain errors (e.g., |x| > 1). + #[inline] + #[must_use] + pub fn acos(&self) -> FpResult { + self.context.acos(&self.repr) + } + + /// Calculate the arctangent of the floating point number. + /// + /// # Panics + /// Panics if the result is not representable as a normal value. + #[inline] + #[must_use] + pub fn atan(&self) -> Self { + self.context.atan(&self.repr).value(&self.context) + } + + /// Calculate the arctangent of y / x. + /// + /// Returns `FpResult` to safely handle special cases like (0,0) or infinities. + #[inline] + #[must_use] + pub fn atan2(&self, x: &Self) -> FpResult { + self.context.atan2(&self.repr, &x.repr) + } +} diff --git a/float/src/round_ops.rs b/float/src/round_ops.rs index c41f1604..c8c0a5d0 100644 --- a/float/src/round_ops.rs +++ b/float/src/round_ops.rs @@ -46,17 +46,27 @@ impl FBig { } // Split the float number at the radix point, assuming it exists (the number is not a integer). - // The method returns (integral part, fractional part, fraction precision). + // The method returns (integral part, fractional part, fractional scale). // // Different from the public `split_at_point()` API, this method doesn't take the ownership of // this number. pub(crate) fn split_at_point_internal(&self) -> (IBig, IBig, usize) { debug_assert!(self.repr.exponent < 0); + let shift = (-self.repr.exponent) as usize; if self.repr.smaller_than_one() { - return (IBig::ZERO, self.repr.significand.clone(), self.context.precision); + // For numbers smaller than 1, the integral part is zero and the stored + // significand is the whole fractional payload. + // + // The third return value is the fractional scale, i.e. the number of + // radix digits after the point. It must be -exponent, because callers + // such as round_fract use it as the denominator exponent B^scale. + // + // This is intentionally not self.context.precision: context precision is + // the significant-digit precision of the float, while this value describes + // the positional scale of the fractional part. + return (IBig::ZERO, self.repr.significand.clone(), shift); } - let shift = (-self.repr.exponent) as usize; let (hi, lo) = split_digits_ref::(&self.repr.significand, shift); (hi, lo, shift) } @@ -119,10 +129,13 @@ impl FBig { /// /// Panics if the number is infinte #[inline] + #[must_use] pub fn fract(&self) -> Self { assert_finite(&self.repr); if self.repr.exponent >= 0 { return Self::ZERO; + } else if self.repr.smaller_than_one() { + return self.clone(); } let (_, lo, precision) = self.split_at_point_internal(); diff --git a/float/tests/round.rs b/float/tests/round.rs index 88b16f59..38ebba22 100644 --- a/float/tests/round.rs +++ b/float/tests/round.rs @@ -1,4 +1,7 @@ -use dashu_float::DBig; +use core::str::FromStr; + +use dashu_base::{Approximation::*, ParseError}; +use dashu_float::{round::Rounding::NoOp, DBig}; mod helper_macros; @@ -124,6 +127,39 @@ fn test_trunc_fract() { assert_eq!(dbig!(12e-2).fract().precision(), 2); } +/// Numbers with |self| < 1 can be stored with a single significand digit and a +/// negative exponent (e.g. 0.009 = 9e-3). Rounding must use `-exponent` as the +/// fractional scale, not `context.precision`; otherwise 9 / 10^1 = 0.9 rounds up. +#[test] +fn test_round_smaller_than_one_uses_exponent_scale() -> Result<(), ParseError> { + let a = DBig::from_str("0.009")?.with_precision(1).unwrap(); + assert_eq!(a.round(), DBig::ZERO); + + let b = DBig::from_str("0.09")?.with_precision(1).unwrap(); + assert_eq!(b.round(), DBig::ZERO); + + let c = DBig::from_str("1e-5")?.with_precision(3).unwrap(); + assert_eq!(c.round(), DBig::ZERO); + assert_eq!(c.to_int(), Inexact(ibig!(0), NoOp)); + + Ok(()) +} + +/// `.fract()` must not inflate context precision to `-exponent` when trailing +/// zeros were normalized away from the significand. +#[test] +fn test_fract_preserves_context_precision() -> Result<(), ParseError> { + let a = DBig::from_str("1e-5")?.with_precision(3).unwrap(); + assert_eq!(a.fract(), a); + assert_eq!(a.fract().precision(), 3); + + let b = DBig::from_str("9e-5")?.with_precision(1).unwrap(); + assert_eq!(b.fract(), b); + assert_eq!(b.fract().precision(), 1); + + Ok(()) +} + #[test] #[should_panic] fn test_floor_inf() { diff --git a/float/tests/trig.rs b/float/tests/trig.rs new file mode 100644 index 00000000..e7274b9e --- /dev/null +++ b/float/tests/trig.rs @@ -0,0 +1,142 @@ +use dashu_float::ops::Abs; +use dashu_float::{round, DBig, FBig, Repr}; + +#[test] +fn test_pi() { + let pi = DBig::pi(10); + assert_eq!(pi.to_string(), "3.141592654"); + + let pi20 = DBig::pi(20); + assert_eq!(pi20.to_string(), "3.1415926535897932385"); + + let pi100 = DBig::pi(100); + assert_eq!(pi100.to_string(), "3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068"); + + let pi500 = DBig::pi(500); + assert_eq!(pi500.to_string(), "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491"); + + let pi_bin = FBig::::pi(100); + assert_eq!(pi_bin.to_string(), "11.00100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111"); +} + +#[test] +fn test_sin_cos() { + let x = DBig::ZERO.with_precision(30).value(); + let (s, c) = x.sin_cos(); + assert_eq!(s, DBig::ZERO); + assert_eq!(c, DBig::ONE.with_precision(30).value()); + + let pi = DBig::pi(30); + let (s, c) = pi.sin_cos(); + assert!(s.abs() < DBig::from_parts(1.into(), -29)); + let neg_one = -DBig::ONE.with_precision(30).value(); + assert!((c - neg_one).abs() < DBig::from_parts(1.into(), -29)); +} + +#[test] +fn test_tan() { + let x = DBig::ZERO.with_precision(30).value(); + assert_eq!(x.tan().value(&x.context()), DBig::ZERO); + + let pi = DBig::pi(30); + let pi4: DBig = pi / 4; + let tan_pi4 = pi4.tan().value(&pi4.context()); + assert!((tan_pi4 - DBig::ONE).abs() < DBig::from_parts(1.into(), -29)); +} + +#[test] +fn test_atan() { + let x = DBig::ZERO.with_precision(30).value(); + assert_eq!(x.atan(), DBig::ZERO); + + let one = DBig::ONE.with_precision(30).value(); + let pi = DBig::pi(30); + let pi4: DBig = pi / 4; + let atan_one = one.atan(); + assert!((atan_one - pi4).abs() < DBig::from_parts(1.into(), -29)); +} + +#[test] +fn test_asin_acos() { + let x = DBig::ZERO.with_precision(30).value(); + assert_eq!(x.asin().value(&x.context()), DBig::ZERO); + + let pi = DBig::pi(30); + let half_pi: DBig = &pi / 2; + assert!((x.acos().value(&x.context()) - half_pi).abs() < DBig::from_parts(1.into(), -29)); + + let half = DBig::from_parts(5.into(), -1).with_precision(30).value(); + let asin_half = half.asin().value(&half.context()); + // asin(0.5) = pi/6 + let pi6: DBig = &pi / 6; + assert!((asin_half - pi6).abs() < DBig::from_parts(1.into(), -29)); + + // Domain error test + let two = DBig::from_parts(2.into(), 0).with_precision(10).value(); + assert!(two.asin().is_nan()); +} + +#[test] +fn test_atan2() { + let zero = DBig::ZERO.with_precision(30).value(); + let one = DBig::ONE.with_precision(30).value(); + let neg_one = -one.clone(); + let pi = DBig::pi(30); + + // atan2(0, 1) = 0 + assert_eq!(zero.atan2(&one).value(&zero.context()), zero); + + // atan2(1, 0) = pi/2 + let half_pi: DBig = &pi / 2; + assert!( + (one.atan2(&zero).value(&one.context()) - half_pi.clone()).abs() + < DBig::from_parts(1.into(), -29) + ); + + // atan2(0, -1) = pi + assert!( + (zero.atan2(&neg_one).value(&zero.context()) - &pi).abs() < DBig::from_parts(1.into(), -29) + ); + + // atan2(-1, 0) = -pi/2 + let m_half_pi: DBig = -half_pi; + assert!( + (neg_one.atan2(&zero).value(&neg_one.context()) - m_half_pi).abs() + < DBig::from_parts(1.into(), -29) + ); + + // Undefined case + let z0 = DBig::ZERO.with_precision(10).value(); + assert!(z0.atan2(&z0).is_nan()); +} + +#[test] +fn test_atan2_infinities() { + let x = DBig::ZERO.with_precision(30).value(); + let ctx = x.context(); + let inf = Repr::infinity(); + let neg_inf = Repr::neg_infinity(); + let pi = ctx.pi::<10>().value(); + let pi_4 = &pi / 4; + let pi_3_4 = &pi * 3 / 4; + + // atan2(+inf, +inf) = pi/4 + let res: DBig = ctx.atan2(&inf, &inf).value(&ctx); + let diff: DBig = res - &pi_4; + assert!(diff.abs() < DBig::from_parts(1.into(), -29)); + + // atan2(+inf, -inf) = 3pi/4 + let res: DBig = ctx.atan2(&inf, &neg_inf).value(&ctx); + let diff: DBig = res - &pi_3_4; + assert!(diff.abs() < DBig::from_parts(1.into(), -29)); + + // atan2(-inf, +inf) = -pi/4 + let res: DBig = ctx.atan2(&neg_inf, &inf).value(&ctx); + let diff: DBig = res + &pi_4; + assert!(diff.abs() < DBig::from_parts(1.into(), -29)); + + // atan2(-inf, -inf) = -3pi/4 + let res: DBig = ctx.atan2(&neg_inf, &neg_inf).value(&ctx); + let diff: DBig = res + &pi_3_4; + assert!(diff.abs() < DBig::from_parts(1.into(), -29)); +} diff --git a/fuzz/Cargo.toml b/fuzz/Cargo.toml new file mode 100644 index 00000000..d065e166 --- /dev/null +++ b/fuzz/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "fuzz" +version = "0.0.0" +edition = "2024" +publish = false + +[dependencies] +dashu-float = { path = "../float" } +dashu-base = { path = "../base" } +rand = "0.10.1" +rug = "1.24" diff --git a/fuzz/src/lib.rs b/fuzz/src/lib.rs new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/fuzz/src/lib.rs @@ -0,0 +1 @@ + diff --git a/fuzz/tests/trig_random.rs b/fuzz/tests/trig_random.rs new file mode 100644 index 00000000..9637369d --- /dev/null +++ b/fuzz/tests/trig_random.rs @@ -0,0 +1,356 @@ +use core::str::FromStr; +use dashu_float::math::FpResult; +use dashu_float::ops::Abs; +use dashu_float::round::mode::HalfEven; +use dashu_float::{DBig, FBig}; +use rand::prelude::*; +use rug::Float; + +/// Reproduction case for a bug discovered during fuzzing where very small +/// numbers with many digits triggered an assertion failure in the rounding logic. +#[test] +#[ignore] +fn test_reproduce_assertion_failure() { + let x_str = "-5.525474318981006776603409487767135633516667011547942409467e-3"; + let prec = 100; + let x_dashu = DBig::from_str(x_str).unwrap().with_rounding::(); + let dashu_ctx = dashu_float::Context::::new(prec); + let _sin_d = dashu_ctx.sin(x_dashu.repr()).value(&dashu_ctx); +} + +#[test] +#[ignore] +fn test_pi_fuzz() { + for prec in (10..1000).step_by(53) { + let pi_dashu = DBig::pi(prec).with_rounding::(); + let bits = (prec * 3322).div_ceil(1000) + 32; + let pi_rug = Float::with_val(bits as u32, rug::float::Constant::Pi); + let s_r_val = DBig::from_str(&pi_rug.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (pi_dashu.clone() - s_r_val).abs() + <= DBig::from_parts(10.into(), -(isize::try_from(prec).unwrap())), + "Pi mismatch at prec={prec}: dashu={pi_dashu}, rug={pi_rug}" + ); + } +} + +/// Generates a truly arbitrary `DBig` value for testing. +fn random_dbig(rng: &mut R, large_exp: bool) -> DBig { + let sign = if rng.random_bool(0.5) { 1 } else { -1 }; + let num_digits = rng.random_range(1..100); + let mut s = String::new(); + if sign == -1 { + s.push('-'); + } + for _ in 0..num_digits { + s.push(char::from_digit(rng.random_range(0..10), 10).unwrap()); + } + let exponent = if large_exp { + rng.random_range(-2000..2000) + } else { + rng.random_range(-10..10) + }; + s.push_str(&format!("e{exponent}")); + DBig::from_str(&s).unwrap_or(DBig::ZERO) +} + +#[test] +#[ignore] +fn test_trig_fuzz_comprehensive() { + let mut rng = StdRng::seed_from_u64(42); + let precisions = [10, 20, 50, 100]; + + for i in 0..2000 { + let x_dashu = random_dbig(&mut rng, true).with_rounding::(); + let x_str = format!("{x_dashu:e}"); + + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + let x_f_repr = x_dashu.repr().clone(); + + // Sin + let sin_d = + match std::panic::catch_unwind(|| dashu_ctx.sin(&x_f_repr).value(&dashu_ctx)) { + Ok(v) => v, + Err(_) => { + panic!("PANIC at iteration {i}, prec {prec}, x = {x_str}"); + } + }; + + // Rug baseline + let x_bits = ((x_dashu.repr().exponent().abs() as f64 * 3.322).ceil() as u32) + 500; + let bits = (((prec as f64).max(100.0) * 3.322).ceil() as u32) + x_bits; + let x_rug = match Float::parse(&x_str) { + Ok(parsed) => Float::with_val(bits, parsed), + Err(_) => continue, + }; + + let sin_r = x_rug.clone().sin(); + let s_r_val = DBig::from_str(&sin_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (sin_d.clone() - s_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Sin mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={sin_d}, rug={sin_r}" + ); + } + } +} + +#[test] +#[ignore] +fn test_atan2_fuzz_comprehensive() { + let mut rng = StdRng::seed_from_u64(45); + let precisions = [20, 50]; + + for i in 0..500 { + let y_dashu = random_dbig(&mut rng, true).with_rounding::(); + let x_dashu = random_dbig(&mut rng, true).with_rounding::(); + let y_str = format!("{y_dashu:e}"); + let x_str = format!("{x_dashu:e}"); + + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + + let atan2_d = std::panic::catch_unwind(|| { + dashu_ctx + .atan2(y_dashu.repr(), x_dashu.repr()) + .value(&dashu_ctx) + }) + .unwrap_or_else(|_| { + panic!("PANIC at iteration {i}, prec {prec}, y = {y_str}, x = {x_str}"); + }); + + let bits = (u32::try_from(prec).unwrap() * 4) + 1000; + let y_rug = Float::with_val(bits, Float::parse(&y_str).unwrap()); + let x_rug = Float::with_val(bits, Float::parse(&x_str).unwrap()); + let atan2_r = y_rug.atan2(&x_rug); + + let a_r_val = DBig::from_str(&atan2_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (atan2_d.clone() - a_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Atan2 mismatch at iteration {i}, y={y_str}, x={x_str}, prec={prec}: dashu={atan2_d}, rug={atan2_r}" + ); + } + } +} + +/// Generates a random `DBig` within [min, max] range. +fn random_dbig_range(rng: &mut R, min: f64, max: f64) -> DBig { + let val: f64 = rng.random_range(min..max); + DBig::from_str(&format!("{val:.15}")).unwrap() +} + +#[test] +#[ignore] +fn test_inv_trig_fuzz() { + let mut rng = StdRng::seed_from_u64(43); + let precisions = [20, 50]; + + for i in 0..200 { + // Test asin/acos within [-1, 1] + let x_dashu = random_dbig_range(&mut rng, -1.0, 1.0).with_rounding::(); + let x_str = format!("{x_dashu:e}"); + + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + + // Asin + let asin_d = dashu_ctx.asin(x_dashu.repr()).value(&dashu_ctx); + let bits = (u32::try_from(prec).unwrap() * 4) + 128; + let x_rug = Float::with_val(bits, Float::parse(&x_str).unwrap()); + let asin_r = x_rug.clone().asin(); + let a_r_val = DBig::from_str(&asin_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (asin_d.clone() - a_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Asin mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={asin_d}, rug={asin_r}" + ); + + // Acos + let acos_d = dashu_ctx.acos(x_dashu.repr()).value(&dashu_ctx); + let acos_r = x_rug.acos(); + let a_r_val = DBig::from_str(&acos_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (acos_d.clone() - a_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Acos mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={acos_d}, rug={acos_r}" + ); + } + } +} + +#[test] +#[ignore] +fn test_edge_cases_fuzz() { + let mut rng = StdRng::seed_from_u64(46); + let precisions = [30, 100]; + + for _ in 0..50 { + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + + // Numbers very close to 1.0 (test asin/acos precision) + let epsilon = 10.0f64.powi(-(rng.random_range(1..15))); + let x_val = 1.0 - epsilon; + let x_dashu = DBig::from_str(&format!("{x_val:.16}")) + .unwrap() + .with_rounding::(); + let x_str = format!("{x_dashu:e}"); + + let asin_d = dashu_ctx.asin(x_dashu.repr()).value(&dashu_ctx); + let bits = (u32::try_from(prec).unwrap() * 4) + 256; + let x_rug = Float::with_val(bits, Float::parse(&x_str).unwrap()); + let asin_r = x_rug.asin(); + let a_r_val = DBig::from_str(&asin_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (asin_d.clone() - a_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Edge Asin mismatch: x={x_str}, prec={prec}" + ); + } + } +} + +#[test] +#[ignore] +fn test_tan_large_exponent_regression() { + let x_str = "-3.67225387623341113999117300261402819219640608e511"; + for prec in [20usize, 50] { + let x_dashu = DBig::from_str(x_str).unwrap().with_rounding::(); + let dashu_ctx = dashu_float::Context::::new(prec); + let tan_d = dashu_ctx.tan(x_dashu.repr()).value(&dashu_ctx); + + let bits = (u32::try_from(prec).unwrap() * 4) + 512 + 1700; // extra bits for large exponent + let x_rug = Float::with_val(bits, Float::parse(x_str).unwrap()); + let tan_r = x_rug.tan(); + let t_r_val = DBig::from_str(&tan_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (tan_d.clone() - t_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Large-exponent tan regression failed at prec={prec}: dashu={tan_d}, rug={tan_r}" + ); + } +} + +#[test] +#[ignore] +fn test_pythagorean_identity_fuzz() { + let mut rng = StdRng::seed_from_u64(99); + let precisions = [20usize, 50, 100]; + + for i in 0..1000 { + let x_dashu = random_dbig(&mut rng, true).with_rounding::(); + + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + let (s, c) = dashu_ctx.sin_cos(x_dashu.repr()); + if let (FpResult::Normal(s_r), FpResult::Normal(c_r)) = (s, c) { + let s_f = FBig::from_repr(s_r.value(), dashu_ctx); + let c_f = FBig::from_repr(c_r.value(), dashu_ctx); + let sum = s_f.clone() * &s_f + c_f.clone() * &c_f; + let one = DBig::ONE + .with_precision(prec) + .value() + .with_rounding::(); + assert!( + (sum.clone() - one).abs() + <= DBig::from_parts(1000.into(), -(isize::try_from(prec).unwrap())), + "sin²+cos²≠1 at iteration {i}, prec={prec}, x={x_dashu:e}, sum={sum}" + ); + } + } + } +} + +#[test] +#[ignore] +fn test_cos_fuzz_comprehensive() { + let mut rng = StdRng::seed_from_u64(47); + let precisions = [10usize, 20, 50, 100]; + + for i in 0..2000 { + let x_dashu = random_dbig(&mut rng, true).with_rounding::(); + let x_str = format!("{x_dashu:e}"); + + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + let cos_d = + std::panic::catch_unwind(|| dashu_ctx.cos(x_dashu.repr()).value(&dashu_ctx)) + .unwrap_or_else(|_| panic!("PANIC at iteration {i}, prec {prec}, x = {x_str}")); + + let x_bits = ((x_dashu.repr().exponent().abs() as f64 * 3.322).ceil() as u32) + 500; + let bits = (((prec as f64).max(100.0) * 3.322).ceil() as u32) + x_bits; + let x_rug = match Float::parse(&x_str) { + Ok(parsed) => Float::with_val(bits, parsed), + Err(_) => continue, + }; + let cos_r = x_rug.cos(); + let c_r_val = DBig::from_str(&cos_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + assert!( + (cos_d.clone() - c_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Cos mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={cos_d}, rug={cos_r}" + ); + } + } +} + +#[test] +#[ignore] +fn test_tan_fuzz_strict() { + let mut rng = StdRng::seed_from_u64(44); + let precisions = [20usize, 50]; + + for i in 0..500 { + let x_dashu = random_dbig(&mut rng, true).with_rounding::(); + let x_str = format!("{x_dashu:e}"); + + for &prec in &precisions { + let dashu_ctx = dashu_float::Context::::new(prec); + + // Only skip if we can verify it's actually near a singularity (|cos| < 10^-5) + let cos_d = dashu_ctx.cos(x_dashu.repr()).value(&dashu_ctx); + if cos_d.abs() < DBig::from_parts(1.into(), -5).with_rounding::() { + continue; + } + + let tan_d = + std::panic::catch_unwind(|| dashu_ctx.tan(x_dashu.repr()).value(&dashu_ctx)) + .unwrap_or_else(|_| panic!("PANIC at iteration {i}, prec {prec}, x = {x_str}")); + + let x_bits = ((x_dashu.repr().exponent().abs() as f64 * 3.322).ceil() as u32) + 500; + let bits = (((prec as f64).max(100.0) * 3.322).ceil() as u32) + x_bits; + let x_rug = match Float::parse(&x_str) { + Ok(parsed) => Float::with_val(bits, parsed), + Err(_) => continue, + }; + let tan_r = x_rug.tan(); + let t_r_val = DBig::from_str(&tan_r.to_string_radix(10, Some(prec))) + .unwrap() + .with_rounding::(); + + assert!( + (tan_d.clone() - t_r_val).abs() + <= DBig::from_parts(100.into(), -(isize::try_from(prec).unwrap())), + "Tan mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={tan_d}, rug={tan_r}" + ); + } + } +}