From 53c7aa421b83bd4cc29f561b2a0adc6f874f0d5d Mon Sep 17 00:00:00 2001 From: Pedro Martins Date: Thu, 23 Apr 2026 04:46:10 +0100 Subject: [PATCH 1/2] feat(float): implement arbitrary-precision trigonometric functions and constants --- float/Cargo.toml | 1 + float/src/lib.rs | 1 + float/src/math/consts.rs | 105 ++++++++ float/src/math/mod.rs | 131 +++++++-- float/src/math/trig.rs | 529 +++++++++++++++++++++++++++++++++++++ float/src/round_ops.rs | 7 +- float/tests/trig.rs | 112 ++++++++ float/tests/trig_random.rs | 119 +++++++++ 8 files changed, 989 insertions(+), 16 deletions(-) create mode 100644 float/src/math/consts.rs create mode 100644 float/src/math/trig.rs create mode 100644 float/tests/trig.rs create mode 100644 float/tests/trig_random.rs diff --git a/float/Cargo.toml b/float/Cargo.toml index d990864..1b32639 100644 --- a/float/Cargo.toml +++ b/float/Cargo.toml @@ -62,6 +62,7 @@ serde_json = { version = "1.0" } postgres = { version = "0.19.4" } criterion = { version = "0.5.1", features = ["html_reports"] } +rug = "1.24" [[test]] name = "random" diff --git a/float/src/lib.rs b/float/src/lib.rs index 1f49a8f..1af2a3d 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 0000000..b2c2cdb --- /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, +}; +use dashu_base::{BitTest, EstimatedLog2}; +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. + pub fn pi(&self) -> FBig { + assert_limited_precision(self.precision); + + // Calculate required bits based on target precision in base B. + // bits = ceil(precision * log2(B)) + let bits = if B == 2 { + self.precision + } else { + let (_, log2_b_ub) = B.log2_bounds(); + (self.precision as f64 * log2_b_ub as f64) 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 { + let (log2_b_lb, _) = B.log2_bounds(); + (work_bits as f64 / log2_b_lb as f64) as usize + 1 + }; + let work_context = Context::::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::(426880.into()).value(); + + let pi = (constant * sqrt_10005 * q_f) / t_f; + pi.with_precision(self.precision).value() + } +} + +/// 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(13591409)); + } + + 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(10939058860032000u64); + let t_val = IBig::from(13591409) + IBig::from(545140134u64) * k; + let t = if a % 2 == 1 { + -(IBig::from(p.clone()) * t_val) + } else { + IBig::from(p.clone()) * t_val + }; + 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] + pub fn pi(precision: usize) -> Self { + Context::::new(precision).pi() + } +} diff --git a/float/src/math/mod.rs b/float/src/math/mod.rs index c642ffb..0c3ecd3 100644 --- a/float/src/math/mod.rs +++ b/float/src/math/mod.rs @@ -1,31 +1,132 @@ -//! 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::{ + repr::{Context, Repr, Word}, + round::Round, + fbig::FBig, +}; -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. +#[derive(Clone, Debug, PartialEq, Eq)] +pub enum FpResult { + Normal(Repr), 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] + pub fn value(self, context: &Context) -> FBig { + match self { + FpResult::Normal(repr) => FBig::new(repr, *context), + FpResult::NaN => panic!("The result of the operation is NaN"), + FpResult::Infinite => panic!("The result of the operation is Infinite"), + FpResult::Overflow => panic!("The result of the operation overflowed"), + FpResult::Underflow => panic!("The result of the operation underflowed"), + } + } + + /// Returns `true` if the result is `NaN`. + #[inline] + pub fn is_nan(&self) -> bool { + matches!(self, FpResult::NaN) + } + + /// Returns `true` if the result is `Infinite`. + #[inline] + pub fn is_infinite(&self) -> bool { + matches!(self, FpResult::Infinite) + } + + /// Returns `true` if the result is a normal finite value. + #[inline] + pub fn is_normal(&self) -> bool { + matches!(self, FpResult::Normal(_)) + } + + /// Returns `true` if the result is a finite value (Normal, Overflow, or Underflow). + #[inline] + pub fn is_finite(&self) -> bool { + matches!(self, FpResult::Normal(_) | FpResult::Overflow | FpResult::Underflow) } } -trait ContextOps { - fn context(&self) -> &Context; - fn repr(&self) -> &Repr; +/// Operations that can be performed on floating point numbers via their context. +pub trait ContextOps { + fn context(&self) -> &Context; + fn repr(&self) -> &Repr; + /// Calculate the sine of the number. #[inline] - fn sin(&self) -> FpResult { + fn sin(&self) -> FpResult { self.context().sin(self.repr()) } -} \ No newline at end of file + + /// Calculate the cosine of the number. + #[inline] + fn cos(&self) -> FpResult { + self.context().cos(self.repr()) + } + + /// Calculate both the sine and cosine of the number. + #[inline] + fn sin_cos(&self) -> (FpResult, FpResult) { + self.context().sin_cos(self.repr()) + } + + /// Calculate the tangent of the number. + #[inline] + fn tan(&self) -> FpResult { + self.context().tan(self.repr()) + } + + /// Calculate the arcsine of the number. + #[inline] + fn asin(&self) -> FpResult { + self.context().asin(self.repr()) + } + + /// Calculate the arccosine of the number. + #[inline] + fn acos(&self) -> FpResult { + self.context().acos(self.repr()) + } + + /// Calculate the arctangent of the number. + #[inline] + fn atan(&self) -> FpResult { + self.context().atan(self.repr()) + } + + /// Calculate the 2-argument arctangent of the number (`y`) and `x`. + #[inline] + fn atan2(&self, x: &Repr) -> FpResult { + self.context().atan2(self.repr(), x) + } +} + +impl ContextOps for FBig { + #[inline] + fn context(&self) -> &Context { + &self.context + } + #[inline] + fn repr(&self) -> &Repr { + &self.repr + } +} diff --git a/float/src/math/trig.rs b/float/src/math/trig.rs new file mode 100644 index 0000000..bf7c991 --- /dev/null +++ b/float/src/math/trig.rs @@ -0,0 +1,529 @@ +use crate::{ + error::assert_limited_precision, + fbig::FBig, + math::FpResult, + repr::{Context, Repr, Word}, + round::Round, +}; +use core::convert::TryFrom; +use dashu_base::{Abs, AbsOrd, Sign}; +use dashu_int::IBig; + +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(&self, x: &Repr) -> Context { + // x_mag estimates m = floor(log_BASE(|x|)) + let x_f_approx = FBig::::new(x.clone(), Context::new(10)); + let x_mag = (x_f_approx.repr.exponent + x_f_approx.repr.digits() as isize).max(0) as usize; + + // We need precision + log10(x) digits to maintain 'precision' digits after reduction. + let guard_digits = x_mag + 50; + Context::::new(self.precision.saturating_add(guard_digits)) + } + + /// Calculate the sine of the floating point representation. + pub fn sin(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + assert_limited_precision(self.precision); + + if x.significand.is_zero() { + let res = FBig::::ZERO.with_precision(self.precision).value(); + return FpResult::Normal(res.repr); + } + + let work_context = self.compute_work_context(x); + + // Round x to work_precision to satisfy internal division constraints and improve performance. + let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + + // 1. Angle reduction: x = k * (pi/2) + r, where |r| <= pi/4 + let pi = work_context.pi::(); + let half_pi = &pi / 2; + let x_scaled: FBig = &x_f / &half_pi; + let k_f = x_scaled.round(); + let k = IBig::try_from(k_f.clone()).expect("rounded x_scaled is not representable as an integer"); + let r = x_f - k_f * half_pi; + + // 2. Identify the quadrant + let k_mod_4_big = &k % IBig::from(4); + let mut k_mod_4 = i8::try_from(k_mod_4_big).unwrap(); + if k_mod_4 < 0 { + k_mod_4 += 4; + } + + // 3. Evaluate the reduced series based on the quadrant + let res = match k_mod_4 { + 0 => work_context.sin_internal(&r), + 1 => work_context.cos_internal(&r), + 2 => -work_context.sin_internal(&r), + 3 => -work_context.cos_internal(&r), + _ => unreachable!(), + }; + FpResult::Normal(res.with_precision(self.precision).value().repr) + } + + /// Calculate the cosine of the floating point representation. + pub fn cos(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + assert_limited_precision(self.precision); + + if x.significand.is_zero() { + let res = FBig::::ONE.with_precision(self.precision).value(); + return FpResult::Normal(res.repr); + } + + let work_context = self.compute_work_context(x); + let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + + let pi = work_context.pi::(); + let half_pi = &pi / 2; + let x_scaled: FBig = &x_f / &half_pi; + let k_f = x_scaled.round(); + let k = IBig::try_from(k_f.clone()).expect("rounded x_scaled is not representable as an integer"); + let r = x_f - k_f * half_pi; + + // 2. Identify the quadrant + let k_mod_4_big = &k % IBig::from(4); + let mut k_mod_4 = i8::try_from(k_mod_4_big).unwrap(); + if k_mod_4 < 0 { + k_mod_4 += 4; + } + + // 3. Evaluate the reduced series based on the quadrant + let res = match k_mod_4 { + 0 => work_context.cos_internal(&r), + 1 => -work_context.sin_internal(&r), + 2 => -work_context.cos_internal(&r), + 3 => work_context.sin_internal(&r), + _ => unreachable!(), + }; + FpResult::Normal(res.with_precision(self.precision).value().repr) + } + + /// Calculate both the sine and cosine of the floating point representation. + /// + /// This is more efficient than calling `sin` and `cos` separately. + 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.significand.is_zero() { + let s = FBig::::ZERO.with_precision(self.precision).value(); + let c = FBig::::ONE.with_precision(self.precision).value(); + return (FpResult::Normal(s.repr), FpResult::Normal(c.repr)); + } + + let work_context = self.compute_work_context(x); + let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + + let pi = work_context.pi::(); + let half_pi = &pi / 2; + let x_scaled: FBig = &x_f / &half_pi; + let k_f = x_scaled.round(); + let k = IBig::try_from(k_f.clone()).expect("rounded x_scaled is not representable as an integer"); + let r = x_f - k_f * half_pi; + + let (sin_r, cos_r) = work_context.sin_cos_internal(&r); + + let k_mod_4_big = &k % IBig::from(4); + let mut k_mod_4 = i8::try_from(k_mod_4_big).unwrap(); + if k_mod_4 < 0 { + k_mod_4 += 4; + } + + let (s, c) = match k_mod_4 { + 0 => (sin_r, cos_r), + 1 => (cos_r, -sin_r), + 2 => (-sin_r, -cos_r), + 3 => (-cos_r, sin_r), + _ => unreachable!(), + }; + + ( + FpResult::Normal(s.with_precision(self.precision).value().repr), + FpResult::Normal(c.with_precision(self.precision).value().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 k % 2 == 1 { + sum -= &term; + } else { + sum += &term; + } + if term.abs_cmp(&threshold).is_le() { + break; + } + k += 1; + } + sum + } + + /// 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 k % 2 == 1 { + sum -= &term; + } else { + sum += &term; + } + if term.abs_cmp(&threshold).is_le() { + break; + } + k += 1; + } + sum + } + + /// 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); + if k % 2 == 1 { + cos_sum -= &cos_term; + } else { + cos_sum += &cos_term; + } + sin_term *= &x2; + sin_term /= (2 * k) * (2 * k + 1); + if k % 2 == 1 { + sin_sum -= &sin_term; + } else { + sin_sum += &sin_term; + } + if sin_term.abs_cmp(&sin_threshold).is_le() + && cos_term.abs_cmp(&cos_threshold).is_le() + { + break; + } + k += 1; + } + (sin_sum, cos_sum) + } + + /// Calculate the tangent of the floating point representation. + /// + /// # Note + /// Near odd multiples of π/2, the result returns `Infinite`. + pub fn tan(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + let (s_res, c_res) = self.sin_cos(x); + if let (FpResult::Normal(s), FpResult::Normal(c)) = (s_res, c_res) { + let s_f = FBig::::new(s, *self); + let c_f = FBig::::new(c, *self); + if c_f.repr.is_zero() { + return FpResult::Infinite; + } + return FpResult::Normal((s_f / c_f).repr); + } + FpResult::NaN + } + + /// 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`. + pub fn asin(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + let guard_digits = 50; + let work_precision = self.precision + guard_digits; + let work_context = Context::::new(work_precision); + + let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + let one = FBig::::ONE.with_precision(work_precision).value(); + + // Domain check: |x| must be <= 1 + if x_f.clone().abs() > one { + return FpResult::NaN; + } + + let x2 = x_f.sqr(); + let d = work_context.sqrt(&(one - x2).repr).value(); + + // If x is exactly 1 or -1, d is zero. + if d.repr.is_zero() { + let pi = work_context.pi::(); + let half_pi: FBig = pi / 2; + let res: FBig = if x_f.sign() == Sign::Positive { half_pi } else { -half_pi }; + return FpResult::Normal(res.with_precision(self.precision).value().repr); + } + + let atan_res = work_context.atan(&(x_f / d).repr); + if let FpResult::Normal(r) = atan_res { + let res = FBig::::new(r, work_context); + return FpResult::Normal(res.with_precision(self.precision).value().repr); + } + FpResult::NaN + } + + /// 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. + pub fn acos(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + return FpResult::NaN; + } + // Use higher precision for BOTH asin and the final subtraction + let guard_digits = 50; + let work_precision = self.precision + guard_digits; + let work_context = Context::::new(work_precision); + + let asin_res = work_context.asin(x); + match asin_res { + FpResult::Normal(asin_repr) => { + let pi = work_context.pi::(); + let half_pi: FBig = pi / 2; + let asin_x = FBig::::new(asin_repr, work_context); + + let res: FBig = half_pi - asin_x; + FpResult::Normal(res.with_precision(self.precision).value().repr) + } + res => res, + } + } + + /// Calculate the arctangent of the floating point representation. + pub fn atan(&self, x: &Repr) -> FpResult { + if x.is_infinite() { + let pi = self.pi::(); + let half_pi: FBig = pi / 2; + let res: FBig = if x.sign() == Sign::Positive { half_pi } else { -half_pi }; + return FpResult::Normal(res.repr); + } + assert_limited_precision(self.precision); + + if x.significand.is_zero() { + let res = FBig::::ZERO.with_precision(self.precision).value(); + return FpResult::Normal(res.repr); + } + + let guard_digits = 50; + let work_precision = self.precision + guard_digits; + let work_context = Context::::new(work_precision); + + let mut x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); + let sign = x_f.sign(); + if sign == Sign::Negative { + x_f = -x_f; + } + + // Range reduction for atan: if |x| >= 1, use atan(x) = pi/2 - atan(1/x) + let mut res = if x_f >= FBig::::ONE.with_precision(work_precision).value() { + let pi = work_context.pi::(); + let inv_x = FBig::ONE / x_f; + (pi / 2) - work_context.atan_internal(&inv_x) + } else { + work_context.atan_internal(&x_f) + }; + + if sign == Sign::Negative { + res = -res; + } + FpResult::Normal(res.with_precision(self.precision).value().repr) + } + + /// Internal series for arctangent using high-stability formula. + 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. + pub fn atan2(&self, y: &Repr, x: &Repr) -> FpResult { + if y.significand.is_zero() && x.significand.is_zero() && y.exponent == 0 && x.exponent == 0 { + return FpResult::NaN; + } + + let guard_digits = 50; + let work_precision = self.precision + guard_digits; + let work_context = Context::::new(work_precision); + + // Handle Infinities according to IEEE 754 + if y.is_infinite() || x.is_infinite() { + let pi = work_context.pi::(); + let half_pi: FBig = pi.clone() / 2; + let res = if y.is_infinite() { + if x.is_infinite() { + // atan2(inf, inf) + if y.sign() == Sign::Positive { + if x.sign() == Sign::Positive { pi.clone() / 4 } else { pi.clone() * 3 / 4 } + } else { + if x.sign() == Sign::Positive { + let pi4: FBig = pi.clone() / 4; + -pi4 + } else { + let pi34: FBig = pi.clone() * 3 / 4; + -pi34 + } + } + } else { + // atan2(inf, finite) + if y.sign() == Sign::Positive { half_pi } else { -half_pi } + } + } else { + // atan2(finite, inf) + if x.sign() == Sign::Positive { + FBig::::ZERO.with_precision(work_precision).value() + } else { + if y.sign() == Sign::Positive { pi } else { -pi } + } + }; + return FpResult::Normal(res.with_precision(self.precision).value().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); + + if x_f > FBig::::ZERO { + self.atan(&(y_f / x_f).repr) + } else if x_f < FBig::::ZERO { + let pi = work_context.pi::(); + let y_sign = y_f.sign(); + let atan_yx_res = self.atan(&(y_f.clone() / x_f).repr); + if let FpResult::Normal(atan_yx_repr) = atan_yx_res { + let atan_yx = FBig::::new(atan_yx_repr, work_context); + let res = if y_sign != Sign::Negative { + atan_yx + pi + } else { + atan_yx - pi + }; + return FpResult::Normal(res.with_precision(self.precision).value().repr); + } + FpResult::NaN + } else { + // x == 0 case + let pi = work_context.pi::(); + let half_pi: FBig = pi / 2; + if y_f > FBig::::ZERO { + FpResult::Normal(half_pi.with_precision(self.precision).value().repr) + } else { + let res = -half_pi; + FpResult::Normal(res.with_precision(self.precision).value().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] 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] 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] 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] 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] 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] 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] 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] 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 c41f160..041b2d3 100644 --- a/float/src/round_ops.rs +++ b/float/src/round_ops.rs @@ -53,7 +53,12 @@ impl FBig { pub(crate) fn split_at_point_internal(&self) -> (IBig, IBig, usize) { debug_assert!(self.repr.exponent < 0); if self.repr.smaller_than_one() { - return (IBig::ZERO, self.repr.significand.clone(), self.context.precision); + // FIX: when the number is smaller than 1.0 (exponent < 0), the fractional part + // starts from the radix point. The number of digits required to represent this + // part is exactly the absolute value of the exponent. + // Returning self.context.precision here was causing assertion failures in round_fract + // when the significand had more digits than the context allowed. + return (IBig::ZERO, self.repr.significand.clone(), (-self.repr.exponent) as usize); } let shift = (-self.repr.exponent) as usize; diff --git a/float/tests/trig.rs b/float/tests/trig.rs new file mode 100644 index 0000000..d732ef9 --- /dev/null +++ b/float/tests/trig.rs @@ -0,0 +1,112 @@ +use dashu_float::{DBig, Repr}; +use dashu_float::ops::Abs; + +#[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"); +} + +#[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_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 ctx = DBig::ZERO.with_precision(30).value().context().clone(); + let inf = Repr::infinity(); + let neg_inf = Repr::neg_infinity(); + let pi = ctx.pi::<10>(); + 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/float/tests/trig_random.rs b/float/tests/trig_random.rs new file mode 100644 index 0000000..338831c --- /dev/null +++ b/float/tests/trig_random.rs @@ -0,0 +1,119 @@ +use dashu_float::DBig; +use dashu_float::ops::Abs; +use dashu_float::round::mode::HalfEven; +use rug::Float; +use rand_v08::prelude::*; +use core::str::FromStr; + +mod helper_macros; + +/// Reproduction case for a bug discovered during fuzzing where very small +/// numbers with many digits triggered an assertion failure in the rounding logic. +#[test] +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] +fn test_pi_fuzz() { + for prec in (10..1000).step_by(53) { + let pi_dashu = DBig::pi(prec).with_rounding::(); + let bits = (prec * 3322 + 999) / 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(), -(prec as isize)), + "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.gen_bool(0.5) { 1 } else { -1 }; + let num_digits = rng.gen_range(1..100); + let mut s = String::new(); + if sign == -1 { s.push('-'); } + for _ in 0..num_digits { + s.push(char::from_digit(rng.gen_range(0..10), 10).unwrap()); + } + let exponent = if large_exp { rng.gen_range(-500..500) } else { rng.gen_range(-10..10) }; + s.push_str(&format!("e{exponent}")); + DBig::from_str(&s).unwrap_or(DBig::ZERO) +} + +#[test] +fn test_trig_fuzz_comprehensive() { + let mut rng = StdRng::seed_from_u64(42); + let precisions = [10, 20, 50, 100]; + + for i in 0..1000 { + 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(), -(prec as isize)), + "Sin mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={sin_d}, rug={sin_r}"); + } + } +} + +#[test] +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 = match std::panic::catch_unwind(|| { + dashu_ctx.atan2(y_dashu.repr(), x_dashu.repr()).value(&dashu_ctx) + }) { + Ok(v) => v, + Err(_) => { + panic!("PANIC at iteration {i}, prec {prec}, y = {y_str}, x = {x_str}"); + } + }; + + let bits = (prec as u32 * 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(), -(prec as isize)), + "Atan2 mismatch at iteration {i}, y={y_str}, x={x_str}, prec={prec}: dashu={atan2_d}, rug={atan2_r}"); + } + } +} From adce9166bb0970183649dc991398ae65538beaf5 Mon Sep 17 00:00:00 2001 From: Pedro Martins Date: Thu, 23 Apr 2026 19:50:59 +0100 Subject: [PATCH 2/2] fix: resolve trigonometric precision loss and optimize performance - Fix catastrophic precision loss in argument reduction for large inputs (x > 10^500) by scaling guard digits with magnitude. - Optimize memory by implementing a zero-clone range reduction logic and avoiding UBig clones in Chudnovsky's algorithm. - Fix 1-ulp errors in tan() by computing sin/cos division within the high-precision work context. - Refactor atan2() to handle signed infinities correctly and standardize FpResult error reporting. - Hardened test suite with aggressive Rug-based fuzzing (2000 iterations) and regression cases for massive exponents. --- float/src/error.rs | 20 +++ float/src/math/consts.rs | 51 +++--- float/src/math/mod.rs | 27 ++- float/src/math/trig.rs | 331 ++++++++++++++++++++++--------------- float/tests/trig.rs | 37 +++-- float/tests/trig_random.rs | 292 ++++++++++++++++++++++++++++---- 6 files changed, 545 insertions(+), 213 deletions(-) diff --git a/float/src/error.rs b/float/src/error.rs index 50ae8b1..aad99b8 100644 --- a/float/src/error.rs +++ b/float/src/error.rs @@ -40,3 +40,23 @@ pub const fn panic_power_negative_base() -> ! { pub(crate) 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/math/consts.rs b/float/src/math/consts.rs index b2c2cdb..3555c92 100644 --- a/float/src/math/consts.rs +++ b/float/src/math/consts.rs @@ -2,26 +2,26 @@ use crate::{ error::assert_limited_precision, fbig::FBig, repr::{Context, Word}, - round::Round, + round::{Round, Rounded}, }; -use dashu_base::{BitTest, EstimatedLog2}; +use dashu_base::{BitTest, EstimatedLog2, 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. + /// 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, + /// 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. - pub fn pi(&self) -> FBig { + pub fn pi(&self) -> Rounded> { assert_limited_precision(self.precision); // Calculate required bits based on target precision in base B. @@ -32,16 +32,16 @@ impl Context { let (_, log2_b_ub) = B.log2_bounds(); (self.precision as f64 * log2_b_ub as f64) 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 { @@ -51,15 +51,17 @@ impl Context { (work_bits as f64 / log2_b_lb as f64) as usize + 1 }; let work_context = Context::::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 sqrt_10005 = work_context + .sqrt(&work_context.convert_int::(10005.into()).value().repr) + .value(); let constant = work_context.convert_int::(426880.into()).value(); - + let pi = (constant * sqrt_10005 * q_f) / t_f; - pi.with_precision(self.precision).value() + pi.with_precision(self.precision) } } @@ -71,24 +73,25 @@ fn chudnovsky_bs(a: usize, b: usize) -> (UBig, UBig, IBig) { if a == 0 { return (UBig::ONE, UBig::ONE, IBig::from(13591409)); } - + 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(10939058860032000u64); let t_val = IBig::from(13591409) + IBig::from(545140134u64) * k; + let t_abs = &p * t_val.unsigned_abs(); let t = if a % 2 == 1 { - -(IBig::from(p.clone()) * t_val) + -IBig::from(t_abs) } else { - IBig::from(p.clone()) * t_val + IBig::from(t_abs) }; 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 @@ -100,6 +103,6 @@ impl FBig { /// Calculate π with the given precision and the default rounding mode. #[inline] pub fn pi(precision: usize) -> Self { - Context::::new(precision).pi() + Context::::new(precision).pi().value() } } diff --git a/float/src/math/mod.rs b/float/src/math/mod.rs index 0c3ecd3..73d236d 100644 --- a/float/src/math/mod.rs +++ b/float/src/math/mod.rs @@ -1,8 +1,9 @@ //! Advanced mathematical functions use crate::{ + error::{panic_infinite, panic_nan, panic_overflow, panic_underflow}, repr::{Context, Repr, Word}, - round::Round, + round::{Round, Rounded}, fbig::FBig, }; @@ -14,9 +15,11 @@ pub mod trig; /// 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(Repr), + Normal(Rounded>), Overflow, Underflow, NaN, @@ -33,11 +36,21 @@ impl FpResult { #[inline] pub fn value(self, context: &Context) -> FBig { match self { - FpResult::Normal(repr) => FBig::new(repr, *context), - FpResult::NaN => panic!("The result of the operation is NaN"), - FpResult::Infinite => panic!("The result of the operation is Infinite"), - FpResult::Overflow => panic!("The result of the operation overflowed"), - FpResult::Underflow => panic!("The result of the operation underflowed"), + FpResult::Normal(rounded) => FBig::new(rounded.value(), *context), + FpResult::NaN => panic_nan(), + FpResult::Infinite => panic_infinite(), + FpResult::Overflow => panic_overflow(), + FpResult::Underflow => panic_underflow(), + } + } + + /// Convert the result into an optional [FBig] with the given context. + /// Returns `None` if the result is not `Normal`. + #[inline] + pub fn ok(self, context: &Context) -> Option>> { + match self { + FpResult::Normal(rounded) => Some(rounded.map(|repr| FBig::new(repr, *context))), + _ => None, } } diff --git a/float/src/math/trig.rs b/float/src/math/trig.rs index bf7c991..5802eaa 100644 --- a/float/src/math/trig.rs +++ b/float/src/math/trig.rs @@ -6,7 +6,7 @@ use crate::{ round::Round, }; use core::convert::TryFrom; -use dashu_base::{Abs, AbsOrd, Sign}; +use dashu_base::{AbsOrd, RemEuclid, Sign}; use dashu_int::IBig; impl Context { @@ -16,12 +16,17 @@ impl Context { /// during range reduction for large inputs. fn compute_work_context(&self, x: &Repr) -> Context { // x_mag estimates m = floor(log_BASE(|x|)) - let x_f_approx = FBig::::new(x.clone(), Context::new(10)); - let x_mag = (x_f_approx.repr.exponent + x_f_approx.repr.digits() as isize).max(0) as usize; - + let x_mag = (x.exponent + x.digits() as isize).max(0) as usize; + // We need precision + log10(x) digits to maintain 'precision' digits after reduction. - let guard_digits = x_mag + 50; - Context::::new(self.precision.saturating_add(guard_digits)) + // 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); + Context::::new(work_precision) } /// Calculate the sine of the floating point representation. @@ -32,29 +37,34 @@ impl Context { assert_limited_precision(self.precision); if x.significand.is_zero() { - let res = FBig::::ZERO.with_precision(self.precision).value(); - return FpResult::Normal(res.repr); + let res = FBig::::ZERO.with_precision(self.precision); + return FpResult::Normal(res.map(|v| v.repr)); } let work_context = self.compute_work_context(x); - + // Round x to work_precision to satisfy internal division constraints and improve performance. let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); - + // 1. Angle reduction: x = k * (pi/2) + r, where |r| <= pi/4 - let pi = work_context.pi::(); + 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 k = IBig::try_from(k_f.clone()).expect("rounded x_scaled is not representable as an integer"); - let r = x_f - k_f * half_pi; + let r = x_f - &k_f * half_pi; + let k = match IBig::try_from(k_f) { + Ok(k) => k, + Err(_) => { + unreachable!("round() always returns an integer and sin() ensures input is finite") + } + }; // 2. Identify the quadrant - let k_mod_4_big = &k % IBig::from(4); - let mut k_mod_4 = i8::try_from(k_mod_4_big).unwrap(); - if k_mod_4 < 0 { - k_mod_4 += 4; - } + let k_mod_4_big = k.rem_euclid(IBig::from(4)); + let k_mod_4 = match i8::try_from(k_mod_4_big) { + Ok(v) => v, + Err(_) => unreachable!("k % 4 is always in [0, 3]"), + }; // 3. Evaluate the reduced series based on the quadrant let res = match k_mod_4 { @@ -64,7 +74,7 @@ impl Context { 3 => -work_context.cos_internal(&r), _ => unreachable!(), }; - FpResult::Normal(res.with_precision(self.precision).value().repr) + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) } /// Calculate the cosine of the floating point representation. @@ -75,26 +85,31 @@ impl Context { assert_limited_precision(self.precision); if x.significand.is_zero() { - let res = FBig::::ONE.with_precision(self.precision).value(); - return FpResult::Normal(res.repr); + let res = FBig::::ONE.with_precision(self.precision); + return FpResult::Normal(res.map(|v| v.repr)); } let work_context = self.compute_work_context(x); let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); - - let pi = work_context.pi::(); + + 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 k = IBig::try_from(k_f.clone()).expect("rounded x_scaled is not representable as an integer"); - let r = x_f - k_f * half_pi; + let r = x_f - &k_f * half_pi; + let k = match IBig::try_from(k_f) { + Ok(k) => k, + Err(_) => { + unreachable!("round() always returns an integer and cos() ensures input is finite") + } + }; // 2. Identify the quadrant - let k_mod_4_big = &k % IBig::from(4); - let mut k_mod_4 = i8::try_from(k_mod_4_big).unwrap(); - if k_mod_4 < 0 { - k_mod_4 += 4; - } + let k_mod_4_big = k.rem_euclid(IBig::from(4)); + let k_mod_4 = match i8::try_from(k_mod_4_big) { + Ok(v) => v, + Err(_) => unreachable!("k % 4 is always in [0, 3]"), + }; // 3. Evaluate the reduced series based on the quadrant let res = match k_mod_4 { @@ -104,7 +119,7 @@ impl Context { 3 => work_context.sin_internal(&r), _ => unreachable!(), }; - FpResult::Normal(res.with_precision(self.precision).value().repr) + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) } /// Calculate both the sine and cosine of the floating point representation. @@ -117,28 +132,33 @@ impl Context { assert_limited_precision(self.precision); if x.significand.is_zero() { - let s = FBig::::ZERO.with_precision(self.precision).value(); - let c = FBig::::ONE.with_precision(self.precision).value(); - return (FpResult::Normal(s.repr), FpResult::Normal(c.repr)); + 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 = self.compute_work_context(x); let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); - - let pi = work_context.pi::(); + + 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 k = IBig::try_from(k_f.clone()).expect("rounded x_scaled is not representable as an integer"); - let r = x_f - k_f * half_pi; + let r = x_f - &k_f * half_pi; + let k = match IBig::try_from(k_f) { + Ok(k) => k, + Err(_) => unreachable!( + "round() always returns an integer and sin_cos() ensures input is finite" + ), + }; let (sin_r, cos_r) = work_context.sin_cos_internal(&r); - let k_mod_4_big = &k % IBig::from(4); - let mut k_mod_4 = i8::try_from(k_mod_4_big).unwrap(); - if k_mod_4 < 0 { - k_mod_4 += 4; - } + let k_mod_4_big = k.rem_euclid(IBig::from(4)); + let k_mod_4 = match i8::try_from(k_mod_4_big) { + Ok(v) => v, + Err(_) => unreachable!("k % 4 is always in [0, 3]"), + }; let (s, c) = match k_mod_4 { 0 => (sin_r, cos_r), @@ -149,8 +169,8 @@ impl Context { }; ( - FpResult::Normal(s.with_precision(self.precision).value().repr), - FpResult::Normal(c.with_precision(self.precision).value().repr), + FpResult::Normal(s.with_precision(self.precision).map(|v| v.repr)), + FpResult::Normal(c.with_precision(self.precision).map(|v| v.repr)), ) } @@ -163,7 +183,7 @@ impl Context { let mut sum = x.clone(); let mut term = x.clone(); let mut k = 1usize; - let threshold = sum.sub_ulp(); + let threshold = sum.sub_ulp(); loop { term *= &x2; term /= (2 * k) * (2 * k + 1); @@ -207,12 +227,12 @@ impl Context { } /// Simultaneously evaluate Taylor series for sine and cosine. - pub(crate) fn sin_cos_internal(&self, x: &FBig) -> (FBig, FBig) { + 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(), - ); + return (FBig::ZERO, FBig::ONE.with_precision(self.precision).value()); } let x2 = x.sqr(); let mut sin_sum = x.clone(); @@ -237,8 +257,7 @@ impl Context { } else { sin_sum += &sin_term; } - if sin_term.abs_cmp(&sin_threshold).is_le() - && cos_term.abs_cmp(&cos_threshold).is_le() + if sin_term.abs_cmp(&sin_threshold).is_le() && cos_term.abs_cmp(&cos_threshold).is_le() { break; } @@ -255,14 +274,16 @@ impl Context { if x.is_infinite() { return FpResult::NaN; } - let (s_res, c_res) = self.sin_cos(x); - if let (FpResult::Normal(s), FpResult::Normal(c)) = (s_res, c_res) { - let s_f = FBig::::new(s, *self); - let c_f = FBig::::new(c, *self); - if c_f.repr.is_zero() { - return FpResult::Infinite; - } - return FpResult::Normal((s_f / c_f).repr); + + let work_context = self.compute_work_context(x); + let (s_res, c_res) = work_context.sin_cos(x); + if let (FpResult::Normal(s_rounded), FpResult::Normal(c_rounded)) = (s_res, c_res) { + let s_f = FBig::::new(s_rounded.value(), work_context); + let c_f = FBig::::new(c_rounded.value(), work_context); + if c_f.repr.is_zero() { + return FpResult::Infinite; + } + return FpResult::Normal((s_f / c_f).with_precision(self.precision).map(|v| v.repr)); } FpResult::NaN } @@ -276,33 +297,37 @@ impl Context { if x.is_infinite() { return FpResult::NaN; } - let guard_digits = 50; + let guard_digits = 50; let work_precision = self.precision + guard_digits; let work_context = Context::::new(work_precision); let x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); let one = FBig::::ONE.with_precision(work_precision).value(); - + // Domain check: |x| must be <= 1 - if x_f.clone().abs() > one { - return FpResult::NaN; + if x_f.abs_cmp(&one).is_gt() { + return FpResult::NaN; } - + let x2 = x_f.sqr(); let d = work_context.sqrt(&(one - x2).repr).value(); - + // If x is exactly 1 or -1, d is zero. if d.repr.is_zero() { - let pi = work_context.pi::(); - let half_pi: FBig = pi / 2; - let res: FBig = if x_f.sign() == Sign::Positive { half_pi } else { -half_pi }; - return FpResult::Normal(res.with_precision(self.precision).value().repr); + let pi = work_context.pi::().value(); + let half_pi: FBig = pi / 2; + let res = if x_f.sign() == Sign::Positive { + half_pi + } else { + -half_pi + }; + return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)); } - + let atan_res = work_context.atan(&(x_f / d).repr); - if let FpResult::Normal(r) = atan_res { - let res = FBig::::new(r, work_context); - return FpResult::Normal(res.with_precision(self.precision).value().repr); + if let FpResult::Normal(r_rounded) = atan_res { + let res = FBig::::new(r_rounded.value(), work_context); + return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)); } FpResult::NaN } @@ -323,56 +348,60 @@ impl Context { let asin_res = work_context.asin(x); match asin_res { - FpResult::Normal(asin_repr) => { - let pi = work_context.pi::(); + FpResult::Normal(asin_rounded) => { + let pi = work_context.pi::().value(); let half_pi: FBig = pi / 2; - let asin_x = FBig::::new(asin_repr, work_context); - + let asin_x = FBig::::new(asin_rounded.value(), work_context); + let res: FBig = half_pi - asin_x; - FpResult::Normal(res.with_precision(self.precision).value().repr) + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) } res => res, } } - /// Calculate the arctangent of the floating point representation. pub fn atan(&self, x: &Repr) -> FpResult { if x.is_infinite() { - let pi = self.pi::(); + 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.repr); + 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.significand.is_zero() { - let res = FBig::::ZERO.with_precision(self.precision).value(); - return FpResult::Normal(res.repr); + 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 = Context::::new(work_precision); - + let mut x_f = FBig::::new(work_context.repr_round(x.clone()).value(), work_context); let sign = x_f.sign(); if sign == Sign::Negative { x_f = -x_f; } - + // Range reduction for atan: if |x| >= 1, use atan(x) = pi/2 - atan(1/x) let mut res = if x_f >= FBig::::ONE.with_precision(work_precision).value() { - let pi = work_context.pi::(); + let pi = work_context.pi::().value(); let inv_x = FBig::ONE / x_f; (pi / 2) - work_context.atan_internal(&inv_x) } else { work_context.atan_internal(&x_f) }; - + if sign == Sign::Negative { res = -res; } - FpResult::Normal(res.with_precision(self.precision).value().repr) + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) } /// Internal series for arctangent using high-stability formula. @@ -399,12 +428,12 @@ impl Context { } /// Calculate the arctangent of y / x. - /// + /// /// Handles signed infinities according to IEEE 754 standards. /// Returns `NaN` if both arguments are zero. pub fn atan2(&self, y: &Repr, x: &Repr) -> FpResult { - if y.significand.is_zero() && x.significand.is_zero() && y.exponent == 0 && x.exponent == 0 { - return FpResult::NaN; + if y.is_zero() && x.is_zero() { + return FpResult::NaN; } let guard_digits = 50; @@ -413,35 +442,47 @@ impl Context { // Handle Infinities according to IEEE 754 if y.is_infinite() || x.is_infinite() { - let pi = work_context.pi::(); + let pi = work_context.pi::().value(); let half_pi: FBig = pi.clone() / 2; let res = if y.is_infinite() { if x.is_infinite() { // atan2(inf, inf) if y.sign() == Sign::Positive { - if x.sign() == Sign::Positive { pi.clone() / 4 } else { pi.clone() * 3 / 4 } + if x.sign() == Sign::Positive { + pi.clone() / 4 + } else { + pi.clone() * 3 / 4 + } } else { - if x.sign() == Sign::Positive { + if x.sign() == Sign::Positive { let pi4: FBig = pi.clone() / 4; -pi4 - } else { + } else { let pi34: FBig = pi.clone() * 3 / 4; -pi34 } } } else { // atan2(inf, finite) - if y.sign() == Sign::Positive { half_pi } else { -half_pi } + if y.sign() == Sign::Positive { + half_pi + } else { + -half_pi + } } } else { // atan2(finite, inf) - if x.sign() == Sign::Positive { - FBig::::ZERO.with_precision(work_precision).value() - } else { - if y.sign() == Sign::Positive { pi } else { -pi } + if x.sign() == Sign::Positive { + FBig::::ZERO.with_precision(work_precision).value() + } else { + if y.sign() == Sign::Positive { + pi + } else { + -pi + } } }; - return FpResult::Normal(res.with_precision(self.precision).value().repr); + 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); @@ -450,28 +491,28 @@ impl Context { if x_f > FBig::::ZERO { self.atan(&(y_f / x_f).repr) } else if x_f < FBig::::ZERO { - let pi = work_context.pi::(); + let pi = work_context.pi::().value(); let y_sign = y_f.sign(); - let atan_yx_res = self.atan(&(y_f.clone() / x_f).repr); - if let FpResult::Normal(atan_yx_repr) = atan_yx_res { - let atan_yx = FBig::::new(atan_yx_repr, work_context); + let atan_yx_res = self.atan(&(y_f / x_f).repr); + if let FpResult::Normal(atan_yx_rounded) = atan_yx_res { + let atan_yx = FBig::::new(atan_yx_rounded.value(), work_context); let res = if y_sign != Sign::Negative { atan_yx + pi } else { atan_yx - pi }; - return FpResult::Normal(res.with_precision(self.precision).value().repr); + return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)); } FpResult::NaN } else { // x == 0 case - let pi = work_context.pi::(); + 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).value().repr) + FpResult::Normal(half_pi.with_precision(self.precision).map(|v| v.repr)) } else { let res = -half_pi; - FpResult::Normal(res.with_precision(self.precision).value().repr) + FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr)) } } } @@ -482,48 +523,70 @@ impl FBig { /// /// # Panics /// Panics if the input is infinite or the result is not representable as a normal value. - #[inline] pub fn sin(&self) -> Self { self.context.sin(&self.repr).value(&self.context) } - + #[inline] + 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] pub fn cos(&self) -> Self { self.context.cos(&self.repr).value(&self.context) } - + #[inline] + 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] pub fn sin_cos(&self) -> (Self, Self) { + #[inline] + 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] pub fn tan(&self) -> FpResult { self.context.tan(&self.repr) } - + #[inline] + 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] pub fn asin(&self) -> FpResult { self.context.asin(&self.repr) } - + #[inline] + 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] pub fn acos(&self) -> FpResult { self.context.acos(&self.repr) } - + #[inline] + 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] pub fn atan(&self) -> Self { self.context.atan(&self.repr).value(&self.context) } - + #[inline] + 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] pub fn atan2(&self, x: &Self) -> FpResult { self.context.atan2(&self.repr, &x.repr) } + #[inline] + pub fn atan2(&self, x: &Self) -> FpResult { + self.context.atan2(&self.repr, &x.repr) + } } diff --git a/float/tests/trig.rs b/float/tests/trig.rs index d732ef9..ceea041 100644 --- a/float/tests/trig.rs +++ b/float/tests/trig.rs @@ -1,11 +1,11 @@ -use dashu_float::{DBig, Repr}; use dashu_float::ops::Abs; +use dashu_float::{DBig, 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"); } @@ -28,7 +28,7 @@ fn test_sin_cos() { 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()); @@ -39,11 +39,11 @@ fn test_tan() { 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 @@ -61,20 +61,28 @@ fn test_atan2() { 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)); - + 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)); - + 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)); + 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(); @@ -83,10 +91,11 @@ fn test_atan2() { #[test] fn test_atan2_infinities() { - let ctx = DBig::ZERO.with_precision(30).value().context().clone(); + 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>(); + let pi = ctx.pi::<10>().value(); let pi_4 = &pi / 4; let pi_3_4 = &pi * 3 / 4; diff --git a/float/tests/trig_random.rs b/float/tests/trig_random.rs index 338831c..a974b31 100644 --- a/float/tests/trig_random.rs +++ b/float/tests/trig_random.rs @@ -1,13 +1,13 @@ -use dashu_float::DBig; +use core::str::FromStr; +use dashu_float::math::FpResult; use dashu_float::ops::Abs; use dashu_float::round::mode::HalfEven; -use rug::Float; +use dashu_float::{DBig, FBig}; use rand_v08::prelude::*; -use core::str::FromStr; - +use rug::Float; mod helper_macros; -/// Reproduction case for a bug discovered during fuzzing where very small +/// Reproduction case for a bug discovered during fuzzing where very small /// numbers with many digits triggered an assertion failure in the rounding logic. #[test] fn test_reproduce_assertion_failure() { @@ -24,9 +24,13 @@ fn test_pi_fuzz() { let pi_dashu = DBig::pi(prec).with_rounding::(); let bits = (prec * 3322 + 999) / 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(), -(prec as isize)), - "Pi mismatch at prec={prec}: dashu={pi_dashu}, rug={pi_rug}"); + 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(), -(prec as isize)), + "Pi mismatch at prec={prec}: dashu={pi_dashu}, rug={pi_rug}" + ); } } @@ -35,11 +39,17 @@ fn random_dbig(rng: &mut R, large_exp: bool) -> DBig { let sign = if rng.gen_bool(0.5) { 1 } else { -1 }; let num_digits = rng.gen_range(1..100); let mut s = String::new(); - if sign == -1 { s.push('-'); } + if sign == -1 { + s.push('-'); + } for _ in 0..num_digits { s.push(char::from_digit(rng.gen_range(0..10), 10).unwrap()); } - let exponent = if large_exp { rng.gen_range(-500..500) } else { rng.gen_range(-10..10) }; + let exponent = if large_exp { + rng.gen_range(-2000..2000) + } else { + rng.gen_range(-10..10) + }; s.push_str(&format!("e{exponent}")); DBig::from_str(&s).unwrap_or(DBig::ZERO) } @@ -48,25 +58,24 @@ fn random_dbig(rng: &mut R, large_exp: bool) -> DBig { fn test_trig_fuzz_comprehensive() { let mut rng = StdRng::seed_from_u64(42); let precisions = [10, 20, 50, 100]; - - for i in 0..1000 { + + 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}"); - } - }; - + 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; @@ -74,11 +83,15 @@ fn test_trig_fuzz_comprehensive() { 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(), -(prec as isize)), - "Sin mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={sin_d}, rug={sin_r}"); + 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(), -(prec as isize)), + "Sin mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={sin_d}, rug={sin_r}" + ); } } } @@ -87,33 +100,244 @@ fn test_trig_fuzz_comprehensive() { 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 = match std::panic::catch_unwind(|| { - dashu_ctx.atan2(y_dashu.repr(), x_dashu.repr()).value(&dashu_ctx) + dashu_ctx + .atan2(y_dashu.repr(), x_dashu.repr()) + .value(&dashu_ctx) }) { Ok(v) => v, Err(_) => { panic!("PANIC at iteration {i}, prec {prec}, y = {y_str}, x = {x_str}"); } }; - + let bits = (prec as u32 * 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::(); + + 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(), -(prec as isize)), "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.gen_range(min..max); + DBig::from_str(&format!("{val:.15}")).unwrap() +} + +#[test] +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 = (prec as u32 * 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(), -(prec as isize)), + "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(), -(prec as isize)), + "Acos mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={acos_d}, rug={acos_r}" + ); + } + } +} + +#[test] +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.gen_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 = (prec as u32 * 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(), -(prec as isize)), + "Edge Asin mismatch: x={x_str}, prec={prec}" + ); + } + } +} + +#[test] +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 = (prec as u32 * 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(), -(prec as isize)), + "Large-exponent tan regression failed at prec={prec}: dashu={tan_d}, rug={tan_r}" + ); + } +} + +#[test] +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(), -(prec as isize)), + "sin²+cos²≠1 at iteration {i}, prec={prec}, x={x_dashu:e}, sum={sum}" + ); + } + } + } +} + +#[test] +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 = match std::panic::catch_unwind(|| { + dashu_ctx.cos(x_dashu.repr()).value(&dashu_ctx) + }) { + Ok(v) => v, + Err(_) => 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(), -(prec as isize)), + "Cos mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={cos_d}, rug={cos_r}" + ); + } + } +} + +#[test] +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 = match std::panic::catch_unwind(|| { + dashu_ctx.tan(x_dashu.repr()).value(&dashu_ctx) + }) { + Ok(v) => v, + Err(_) => 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(), -(prec as isize)), + "Tan mismatch at iteration {i}, x={x_str}, prec={prec}: dashu={tan_d}, rug={tan_r}" + ); + } + } +}