feat(float): implement arbitrary-precision trigonometric functions and constants

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"

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!")
+}

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;

float/src/math/consts.rs

@@ -0,0 +1,108 @@
+use crate::{
+    error::assert_limited_precision,
+    fbig::FBig,
+    repr::{Context, Word},
+    round::{Round, Rounded},
+};
+use dashu_base::{BitTest, EstimatedLog2, UnsignedAbs};
+use dashu_int::{IBig, UBig};
+
+impl<R: Round> Context<R> {
+    /// 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<const B: Word>(&self) -> Rounded<FBig<R, B>> {
+        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::<R>::new(work_precision);
+
+        let q_f = work_context.convert_int::<B>(q.into()).value();
+        let t_f = work_context.convert_int::<B>(t).value();
+
+        let sqrt_10005 = work_context
+            .sqrt(&work_context.convert_int::<B>(10005.into()).value().repr)
+            .value();
+        let constant = work_context.convert_int::<B>(426880.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(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(t_abs)
+        } else {
+            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
+    let t = IBig::from(q_r) * t_l + IBig::from(p_l) * t_r;
+    (p, q, t)
+}
+
+impl<R: Round, const B: Word> FBig<R, B> {
+    /// Calculate π with the given precision and the default rounding mode.
+    #[inline]
+    pub fn pi(precision: usize) -> Self {
+        Context::<R>::new(precision).pi().value()
+    }
+}

float/src/math/mod.rs

@@ -1,31 +1,145 @@
-//! 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},
+    repr::{Context, Repr, Word},
+    round::{Round, Rounded},
+    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.
+/// 
+/// Finite results are wrapped in a [Rounded] to preserve rounding information.
+#[derive(Clone, Debug, PartialEq, Eq)]
+pub enum FpResult<const B: Word> {
+    Normal(Rounded<Repr<B>>),
     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<const B: Word> FpResult<B> {
+    /// Convert the result into an [FBig] with the given context.
+    ///
+    /// # Panics
+    /// Panics if the result is not `Normal`.
+    #[inline]
+    pub fn value<R: Round>(self, context: &Context<R>) -> FBig<R, B> {
+        match self {
+            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<R: Round>(self, context: &Context<R>) -> Option<Rounded<FBig<R, B>>> {
+        match self {
+            FpResult::Normal(rounded) => Some(rounded.map(|repr| FBig::new(repr, *context))),
+            _ => None,
+        }
+    }
+
+    /// 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<R: Round, const B: Word> {
+    fn context(&self) -> &Context<R>;
+    fn repr(&self) -> &Repr<B>;
 
+    /// Calculate the sine of the number.
     #[inline]
-    fn sin(&self) -> FpResult {
+    fn sin(&self) -> FpResult<B> {
         self.context().sin(self.repr())
     }
-}
+
+    /// Calculate the cosine of the number.
+    #[inline]
+    fn cos(&self) -> FpResult<B> {
+        self.context().cos(self.repr())
+    }
+
+    /// Calculate both the sine and cosine of the number.
+    #[inline]
+    fn sin_cos(&self) -> (FpResult<B>, FpResult<B>) {
+        self.context().sin_cos(self.repr())
+    }
+
+    /// Calculate the tangent of the number.
+    #[inline]
+    fn tan(&self) -> FpResult<B> {
+        self.context().tan(self.repr())
+    }
+
+    /// Calculate the arcsine of the number.
+    #[inline]
+    fn asin(&self) -> FpResult<B> {
+        self.context().asin(self.repr())
+    }
+
+    /// Calculate the arccosine of the number.
+    #[inline]
+    fn acos(&self) -> FpResult<B> {
+        self.context().acos(self.repr())
+    }
+
+    /// Calculate the arctangent of the number.
+    #[inline]
+    fn atan(&self) -> FpResult<B> {
+        self.context().atan(self.repr())
+    }
+
+    /// Calculate the 2-argument arctangent of the number (`y`) and `x`.
+    #[inline]
+    fn atan2(&self, x: &Repr<B>) -> FpResult<B> {
+        self.context().atan2(self.repr(), x)
+    }
+}
+
+impl<R: Round, const B: Word> ContextOps<R, B> for FBig<R, B> {
+    #[inline]
+    fn context(&self) -> &Context<R> {
+        &self.context
+    }
+    #[inline]
+    fn repr(&self) -> &Repr<B> {
+        &self.repr
+    }
+}