num_valid/backends/rug/

raw.rs

#![deny(rustdoc::broken_intra_doc_links)]

//! # Arbitrary-Precision Floating-Point Raw Implementations
//!
//! This module provides raw trait implementations for the [`rug`](https://docs.rs/rug/) backend,
//! enabling arbitrary-precision floating-point arithmetic with configurable precision.
//!
//! ## Purpose and Role
//!
//! The primary role of this module is to implement the core raw traits
//! ([`RawScalarTrait`](crate::core::traits::raw::RawScalarTrait),
//! [`RawRealTrait`](crate::core::traits::raw::RawRealTrait), etc.) for [`rug::Float`]
//! and [`rug::Complex`]. These implementations provide the computational foundations
//! that validated wrappers build upon.
//!
//! ## Precision Model
//!
//! Unlike the native `f64` backend with fixed 53-bit precision, the rug backend uses
//! **const generic precision**:
//!
//! - Precision is specified in bits (e.g., 100, 200, 500)
//! - All operations preserve precision of the operands
//! - Precision mismatches are detected during validation
//!
//! ## MPFR Constant Optimization
//!
//! For mathematical constants, this module uses MPFR's precomputed constants
//! when available, providing significant performance improvements:
//!
//! - [`rug::float::Constant::Pi`] - ~10x faster than computing `acos(-1)`
//! - [`rug::float::Constant::Log2`] - ~10x faster than computing `ln(2)`
//!
//! ## Memory and Performance Characteristics
//!
//! - `rug::Float` and `rug::Complex` are heap-allocated, non-`Copy` types
//! - All operations support reference-based variants to minimize cloning
//! - Memory usage scales with precision (approximately `precision/8` bytes per value)

use crate::{
    core::{
        errors::{
            ErrorsRawRealToInteger, ErrorsValidationRawComplex, ErrorsValidationRawReal,
            capture_backtrace,
        },
        traits::{
            raw::{
                RawComplexTrait, RawRealTrait, RawScalarHyperbolic, RawScalarPow, RawScalarTrait,
                RawScalarTrigonometric, RawScalarValidate,
            },
            validation::FpChecks,
        },
    },
    functions::{Conjugate, NegAssign, Rounding, Sign},
    kernels::PrecisionBits,
};
use duplicate::duplicate_item;
use rug::{
    float::Constant as MpfrConstant,
    ops::{CompleteRound, Pow},
};
use std::{
    cmp::Ordering,
    hash::{Hash, Hasher},
    num::FpCategory,
    ops::Neg,
};

#[duplicate_item(
    T;
    [rug::Float];
    [rug::Complex];
)]
impl RawScalarTrigonometric for T {
    #[duplicate_item(
        unchecked_method method;
        [unchecked_sin]  [sin];
        [unchecked_asin] [asin];
        [unchecked_cos]  [cos];
        [unchecked_acos] [acos];
        [unchecked_tan]  [tan];
        [unchecked_atan] [atan];)]
    #[inline(always)]
    fn unchecked_method(self) -> Self {
        T::method(self)
    }
}

#[duplicate_item(
    T;
    [rug::Float];
    [rug::Complex];
)]
impl RawScalarHyperbolic for T {
    #[duplicate_item(
        unchecked_method  method;
        [unchecked_sinh]  [sinh];
        [unchecked_asinh] [asinh];
        [unchecked_cosh]  [cosh];
        [unchecked_acosh] [acosh];
        [unchecked_tanh]  [tanh];
        [unchecked_atanh] [atanh];)]
    #[inline(always)]
    fn unchecked_method(self) -> Self {
        T::method(self)
    }
}

#[duplicate_item(
    T;
    [rug::Float];
    [rug::Complex];
)]
impl RawScalarPow for T {
    #[duplicate_item(
        unchecked_method               exponent_type;
        [unchecked_pow_exponent_i8]    [i8];
        [unchecked_pow_exponent_i16]   [i16];
        [unchecked_pow_exponent_i32]   [i32];
        [unchecked_pow_exponent_i64]   [i64];
        [unchecked_pow_exponent_i128]  [i128];
        [unchecked_pow_exponent_isize] [isize];
        [unchecked_pow_exponent_u8]    [u8];
        [unchecked_pow_exponent_u16]   [u16];
        [unchecked_pow_exponent_u32]   [u32];
        [unchecked_pow_exponent_u64]   [u64];
        [unchecked_pow_exponent_u128]  [u128];
        [unchecked_pow_exponent_usize] [usize];
    )]
    #[inline(always)]
    fn unchecked_method(self, exponent: &exponent_type) -> Self {
        T::pow(self, exponent)
    }
}

impl RawScalarValidate for rug::Float {
    type ValidationError = ErrorsValidationRawReal<rug::Float>;

    #[inline(always)]
    fn validate_raw_scalar_ref(&self) -> Result<(), Self::ValidationError> {
        match self.raw_classify() {
            FpCategory::Nan => Err(ErrorsValidationRawReal::<rug::Float>::IsNaN {
                value: self.clone(),
                backtrace: capture_backtrace(),
            }),
            FpCategory::Infinite => {
                if self.is_sign_positive() {
                    Err(ErrorsValidationRawReal::<rug::Float>::IsPosInfinity {
                        backtrace: capture_backtrace(),
                    })
                } else {
                    Err(ErrorsValidationRawReal::<rug::Float>::IsNegInfinity {
                        backtrace: capture_backtrace(),
                    })
                }
            }
            FpCategory::Subnormal => {
                unreachable!("rug::Float should not produce subnormal values");
            }
            FpCategory::Zero | FpCategory::Normal => {
                // No action needed
                Ok(())
            }
        }
    }
}

impl RawScalarTrait for rug::Float {
    fn raw_zero(precision: u32) -> Self {
        rug::Float::with_val(precision, 0.)
    }

    fn is_zero(&self) -> bool {
        rug::Float::is_zero(self)
    }

    fn raw_one(precision: u32) -> Self {
        rug::Float::with_val(precision, 1.)
    }

    #[duplicate_item(
        unchecked_method       method;
        [unchecked_reciprocal] [recip];
        [unchecked_exp]        [exp];
        [unchecked_sqrt]       [sqrt];
        [unchecked_ln]         [ln];
        [unchecked_log2]       [log2];
        [unchecked_log10]      [log10];
    )]
    #[inline(always)]
    fn unchecked_method(self) -> Self {
        rug::Float::method(self)
    }

    /// Multiplies and adds in one fused operation, rounding to the nearest with only one rounding error.
    ///
    /// `a.mul_add(b, c)` produces a result like `a * &b + &c`.
    #[inline(always)]
    fn unchecked_mul_add(self, b: &Self, c: &Self) -> Self {
        rug::Float::mul_add(self, b, c)
    }

    #[inline(always)]
    fn compute_hash_raw_scalar<H: Hasher>(&self, state: &mut H) {
        debug_assert!(
            self.is_finite(),
            "Hashing a non-finite rug::Float value (i.e., NaN or Infinity) may lead to inconsistent results."
        );
        // Always include precision in the hash
        self.prec().hash(state);

        // Handle signed zeros by normalizing to positive zero
        if self.is_zero() {
            // Ensure that -0.0 and 0.0 hash to the same value
            rug::Float::with_val(self.prec(), 0.0)
                .to_string_radix(10, None)
                .hash(state);
        } else {
            // For non-zero values, create a deterministic string representation
            // Use enough decimal places to capture the full precision
            self.to_string_radix(10, None).hash(state)
        }
    }
}

impl PrecisionBits for rug::Float {
    #[inline(always)]
    fn precision_bits(&self) -> u32 {
        self.prec()
    }
}

impl RawRealTrait for rug::Float {
    type RawComplex = rug::Complex;

    #[inline(always)]
    fn unchecked_abs(self) -> rug::Float {
        rug::Float::abs(self)
    }

    #[inline(always)]
    fn unchecked_atan2(self, denominator: &Self) -> Self {
        rug::Float::atan2(self, denominator)
    }

    #[inline(always)]
    fn unchecked_pow_exponent_real(self, exponent: &Self) -> Self {
        rug::Float::pow(self, exponent)
    }

    #[inline(always)]
    fn unchecked_hypot(self, other: &Self) -> Self {
        rug::Float::hypot(self, other)
    }

    #[inline(always)]
    fn unchecked_ln_1p(self) -> Self {
        rug::Float::ln_1p(self)
    }

    #[inline(always)]
    fn unchecked_exp_m1(self) -> Self {
        rug::Float::exp_m1(self)
    }

    /// Multiplies two pairs and adds them in one fused operation, rounding to the nearest with only one rounding error.
    /// `a.unchecked_mul_add_mul_mut(&b, &c, &d)` produces a result like `&a * &b + &c * &d`, but stores the result in `a` using its precision.
    #[inline(always)]
    fn unchecked_mul_add_mul_mut(&mut self, mul: &Self, add_mul1: &Self, add_mul2: &Self) {
        self.mul_add_mul_mut(mul, add_mul1, add_mul2);
    }

    /// Multiplies two pairs and subtracts them in one fused operation, rounding to the nearest with only one rounding error.
    /// `a.unchecked_mul_sub_mul_mut(&b, &c, &d)` produces a result like `&a * &b - &c * &d`, but stores the result in `a` using its precision.
    #[inline(always)]
    fn unchecked_mul_sub_mul_mut(&mut self, mul: &Self, sub_mul1: &Self, sub_mul2: &Self) {
        self.mul_sub_mul_mut(mul, sub_mul1, sub_mul2);
    }

    #[inline(always)]
    fn raw_total_cmp(&self, other: &Self) -> Ordering {
        rug::Float::total_cmp(self, other)
    }

    /// Clamps the value within the specified bounds.
    #[inline(always)]
    fn raw_clamp(self, min: &Self, max: &Self) -> Self {
        rug::Float::clamp(self, min, max)
    }

    #[inline(always)]
    fn raw_classify(&self) -> FpCategory {
        rug::Float::classify(self)
    }

    #[inline(always)]
    fn raw_two(precision: u32) -> Self {
        rug::Float::with_val(precision, 2.)
    }

    #[inline(always)]
    fn raw_one_div_2(precision: u32) -> Self {
        rug::Float::with_val(precision, 0.5)
    }

    #[inline(always)]
    fn raw_pi(precision: u32) -> Self {
        // Use MPFR's optimized precomputed π constant (10x faster than acos(-1))
        rug::Float::with_val(precision, MpfrConstant::Pi)
    }

    #[inline(always)]
    fn raw_two_pi(precision: u32) -> Self {
        rug::Float::with_val(precision, MpfrConstant::Pi) * 2
    }

    #[inline(always)]
    fn raw_pi_div_2(precision: u32) -> Self {
        rug::Float::with_val(precision, MpfrConstant::Pi) / 2
    }

    #[inline(always)]
    fn raw_max_finite(precision: u32) -> Self {
        // Create a Float with the desired precision
        // let f = rug::Float::new(PRECISION);

        // Fill the significand with all 1s: 1 - 2^(-precision)
        // This gives you the largest significand before it rolls over
        let one = rug::Float::with_val(precision, 1);
        let eps = rug::Float::with_val(precision, rug::Float::u_pow_u(2, precision)).recip();
        let significand = one - &eps;

        // Scale it to the maximum exponent (subtract 1 to avoid overflow to infinity)
        let max_exp = rug::float::exp_max() - 1;
        //        println!("max_exp = {max_exp}");
        significand * rug::Float::with_val(precision, rug::Float::u_pow_u(2, max_exp as u32))
    }

    #[inline(always)]
    fn raw_min_finite(precision: u32) -> Self {
        Self::raw_max_finite(precision).neg()
    }

    #[inline(always)]
    fn raw_epsilon(precision: u32) -> Self {
        rug::Float::u_pow_u(2, precision - 1)
            .complete(precision)
            .recip()
    }

    #[inline(always)]
    fn raw_ln_2(precision: u32) -> Self {
        // Use MPFR's optimized precomputed ln(2) constant
        rug::Float::with_val(precision, MpfrConstant::Log2)
    }

    #[inline(always)]
    fn raw_ln_10(precision: u32) -> Self {
        // ln(10) = ln(2) * log₂(10)
        // More efficient than computing ln(10) directly
        let ln2 = rug::Float::with_val(precision, MpfrConstant::Log2);
        let log2_10 = rug::Float::with_val(precision, 10).log2();
        ln2 * log2_10
    }

    #[inline(always)]
    fn raw_log10_2(precision: u32) -> Self {
        rug::Float::with_val(precision, 2.).log10()
    }

    #[inline(always)]
    fn raw_log2_10(precision: u32) -> Self {
        rug::Float::with_val(precision, 10.).log2()
    }

    #[inline(always)]
    fn raw_log2_e(precision: u32) -> Self {
        // log₂(e) = 1/ln(2)
        // More efficient than computing e then taking log₂
        rug::Float::with_val(precision, MpfrConstant::Log2).recip()
    }

    #[inline(always)]
    fn raw_log10_e(precision: u32) -> Self {
        // log₁₀(e) = 1/ln(10) = 1/(ln(2) * log₂(10))
        let ln2 = rug::Float::with_val(precision, MpfrConstant::Log2);
        let log2_10 = rug::Float::with_val(precision, 10).log2();
        (ln2 * log2_10).recip()
    }

    #[inline(always)]
    fn raw_e(precision: u32) -> Self {
        rug::Float::with_val(precision, 1.).exp()
    }

    #[inline(always)]
    fn truncate_to_usize(self) -> Result<usize, ErrorsRawRealToInteger<rug::Float, usize>> {
        if !self.is_finite() {
            return Err(ErrorsRawRealToInteger::NotFinite {
                value: self,
                backtrace: capture_backtrace(),
            });
        }

        let truncation = self.clone().trunc();

        let out_of_range = || ErrorsRawRealToInteger::OutOfRange {
            value: self,
            min: usize::MIN,
            max: usize::MAX,
            backtrace: capture_backtrace(),
        };

        let truncation_as_rug_integer = truncation.to_integer().ok_or_else(out_of_range.clone())?;

        truncation_as_rug_integer
            .to_usize()
            .ok_or_else(out_of_range)
    }
}

impl RawScalarValidate for rug::Complex {
    type ValidationError = ErrorsValidationRawComplex<ErrorsValidationRawReal<rug::Float>>;

    fn validate_raw_scalar_ref(&self) -> Result<(), Self::ValidationError> {
        let real_validation = self.real().validate_raw_scalar_ref();
        let imag_validation = self.imag().validate_raw_scalar_ref();
        match (real_validation, imag_validation) {
            (Ok(()), Ok(())) => Ok(()),
            (Ok(()), Err(imag_err)) => Err(ErrorsValidationRawComplex::InvalidImaginaryPart {
                source: Box::new(imag_err),
            }),
            (Err(real_err), Ok(())) => Err(ErrorsValidationRawComplex::InvalidRealPart {
                source: Box::new(real_err),
            }),
            (Err(real_err), Err(imag_err)) => Err(ErrorsValidationRawComplex::InvalidBothParts {
                real_error: Box::new(real_err),
                imag_error: Box::new(imag_err),
            }),
        }
    }
}

impl RawScalarTrait for rug::Complex {
    fn raw_zero(precision: u32) -> Self {
        rug::Complex::with_val(precision, (0., 0.))
    }

    fn is_zero(&self) -> bool {
        rug::Complex::is_zero(self)
    }

    fn raw_one(precision: u32) -> Self {
        rug::Complex::with_val(precision, (1., 0.))
    }

    #[duplicate_item(
        unchecked_method       method;
        [unchecked_reciprocal] [recip];
        [unchecked_exp]        [exp];
        [unchecked_sqrt]       [sqrt];
        [unchecked_ln]         [ln];
        [unchecked_log10]      [log10];
    )]
    #[inline(always)]
    fn unchecked_method(self) -> Self {
        rug::Complex::method(self)
    }

    #[inline(always)]
    fn unchecked_log2(self) -> Self {
        let ln_2 = rug::Float::with_val(self.real().prec(), 2.).ln();
        rug::Complex::ln(self) / ln_2
    }

    /// Multiplies and adds in one fused operation, rounding to the nearest with only one rounding error.
    ///
    /// `a.mul_add(b, c)` produces a result like `a * &b + &c`.
    #[inline(always)]
    fn unchecked_mul_add(self, b: &Self, c: &Self) -> Self {
        rug::Complex::mul_add(self, b, c)
    }

    fn compute_hash_raw_scalar<H: Hasher>(&self, state: &mut H) {
        self.raw_real_part().compute_hash_raw_scalar(state);
        self.raw_imag_part().compute_hash_raw_scalar(state);
    }
}

impl Conjugate for rug::Complex {
    #[inline(always)]
    fn conjugate(self) -> Self {
        rug::Complex::conj(self)
    }
}

impl RawComplexTrait for rug::Complex {
    type RawReal = rug::Float;

    fn new_unchecked_raw_complex(real: rug::Float, imag: rug::Float) -> Self {
        debug_assert_eq!(
            real.prec(),
            imag.prec(),
            "Different precision between real and imaginary part!"
        );
        rug::Complex::with_val(real.prec(), (real, imag))
    }

    /// Returns a mutable reference to the real part of the complex number.
    fn mut_raw_real_part(&mut self) -> &mut rug::Float {
        self.mut_real()
    }

    /// Returns a mutable reference to the imaginary part of the complex number.
    fn mut_raw_imag_part(&mut self) -> &mut rug::Float {
        self.mut_imag()
    }

    #[inline(always)]
    fn unchecked_abs(self) -> rug::Float {
        rug::Complex::abs(self).into_real_imag().0
    }

    #[inline(always)]
    fn raw_real_part(&self) -> &rug::Float {
        self.real()
    }

    #[inline(always)]
    fn raw_imag_part(&self) -> &rug::Float {
        self.imag()
    }

    #[inline(always)]
    fn unchecked_arg(self) -> rug::Float {
        rug::Complex::arg(self).into_real_imag().0
    }

    #[inline(always)]
    fn unchecked_pow_exponent_real(self, exponent: &rug::Float) -> Self {
        rug::Complex::pow(self, exponent)
    }
}

impl FpChecks for rug::Float {
    fn is_finite(&self) -> bool {
        rug::Float::is_finite(self)
    }

    fn is_infinite(&self) -> bool {
        rug::Float::is_infinite(self)
    }

    fn is_nan(&self) -> bool {
        rug::Float::is_nan(self)
    }
    fn is_normal(&self) -> bool {
        rug::Float::is_normal(self)
    }
}

impl FpChecks for rug::Complex {
    /// Returns `true` if the real and imaginary parts of `self` are neither infinite nor NaN.
    #[inline(always)]
    fn is_finite(&self) -> bool {
        self.real().is_finite() && self.imag().is_finite()
    }

    /// Returns `true` if the real or the imaginary part of `self` are positive infinity or negative infinity.
    #[inline(always)]
    fn is_infinite(&self) -> bool {
        !self.is_nan() && (self.real().is_infinite() || self.imag().is_infinite())
    }

    /// Returns `true` if the real or the imaginary part of `self` are NaN.
    #[inline(always)]
    fn is_nan(&self) -> bool {
        self.real().is_nan() || self.imag().is_nan()
    }

    /// Returns `true` if the real and the imaginary part of `self` are *normal* (i.e. neither zero, infinite, subnormal, or NaN).
    #[inline(always)]
    fn is_normal(&self) -> bool {
        self.real().is_normal() && self.imag().is_normal()
    }
}

//------------------------------------------------------------------------------------------------
#[duplicate_item(
    T;
    [rug::Float];
    [rug::Complex];
)]
/// Compound negation and assignment.
impl NegAssign for T {
    /// Performs the negation of `self`.
    fn neg_assign(&mut self) {
        <T as rug::ops::NegAssign>::neg_assign(self);
    }
}
//------------------------------------------------------------------------------------------------

impl Sign for rug::Float {
    /// Returns a number with the magnitude of `self` and the sign of `sign`.
    #[inline(always)]
    fn kernel_copysign(self, sign: &Self) -> Self {
        self.copysign(sign)
    }

    /// Returns `true` if the value is negative, −0 or NaN with a negative sign.
    #[inline(always)]
    fn kernel_is_sign_negative(&self) -> bool {
        self.is_sign_negative()
    }

    /// Returns `true` if the value is positive, +0 or NaN with a positive sign.
    #[inline(always)]
    fn kernel_is_sign_positive(&self) -> bool {
        self.is_sign_positive()
    }

    /// Returns the signum of the number.
    ///
    /// This method relies on `rug::Float::signum()`.
    /// The behavior for special values is as follows:
    /// - If the number is positive and not zero, returns `1.0`.
    /// - If the number is negative and not zero, returns `-1.0`.
    /// - If the number is zero, `rug::Float::signum()` returns `1.0` (representing positive zero).
    #[inline(always)]
    fn kernel_signum(self) -> Self {
        self.signum()
    }
}

impl Rounding for rug::Float {
    /// Returns the smallest integer greater than or equal to `self`.
    #[inline(always)]
    fn kernel_ceil(self) -> Self {
        self.ceil()
    }

    /// Returns the largest integer smaller than or equal to `self`.
    #[inline(always)]
    fn kernel_floor(self) -> Self {
        self.floor()
    }

    /// Returns the fractional part of `self`.
    #[inline(always)]
    fn kernel_fract(self) -> Self {
        self.fract()
    }

    /// Rounds `self` to the nearest integer, rounding half-way cases away from zero.
    #[inline(always)]
    fn kernel_round(self) -> Self {
        self.round()
    }

    /// Returns the nearest integer to a number. Rounds half-way cases to the number with an even least significant digit.
    ///
    /// This function always returns the precise result.
    ///
    /// # Examples
    /// ```
    /// use num_valid::{RealScalar, RealRugStrictFinite, functions::Rounding};
    ///
    /// const PRECISION: u32 = 100;
    ///
    /// let f = RealRugStrictFinite::<PRECISION>::try_from_f64(3.3).unwrap();
    /// let g = RealRugStrictFinite::<PRECISION>::try_from_f64(-3.3).unwrap();
    /// let h = RealRugStrictFinite::<PRECISION>::try_from_f64(3.5).unwrap();
    /// let i = RealRugStrictFinite::<PRECISION>::try_from_f64(-4.5).unwrap();
    ///
    /// assert_eq!(f.kernel_round_ties_even(), RealRugStrictFinite::<PRECISION>::try_from_f64(3.).unwrap());
    /// assert_eq!(g.kernel_round_ties_even(), RealRugStrictFinite::<PRECISION>::try_from_f64(-3.).unwrap());
    /// assert_eq!(h.kernel_round_ties_even(), RealRugStrictFinite::<PRECISION>::try_from_f64(4.).unwrap());
    /// assert_eq!(i.kernel_round_ties_even(), RealRugStrictFinite::<PRECISION>::try_from_f64(-4.).unwrap());
    /// ```
    #[inline(always)]
    fn kernel_round_ties_even(self) -> Self {
        self.round_even()
    }

    /// Returns the integer part of `self`. This means that non-integer numbers are always truncated towards zero.
    ///    
    /// # Examples
    /// ```
    /// use num_valid::{RealScalar, RealRugStrictFinite, functions::Rounding};
    ///
    /// const PRECISION: u32 = 100;
    ///
    /// let f = RealRugStrictFinite::<PRECISION>::try_from_f64(3.7).unwrap();
    /// let g = RealRugStrictFinite::<PRECISION>::try_from_f64(3.).unwrap();
    /// let h = RealRugStrictFinite::<PRECISION>::try_from_f64(-3.7).unwrap();
    ///
    /// assert_eq!(f.kernel_trunc(), RealRugStrictFinite::<PRECISION>::try_from_f64(3.).unwrap());
    /// assert_eq!(g.kernel_trunc(), RealRugStrictFinite::<PRECISION>::try_from_f64(3.).unwrap());
    /// assert_eq!(h.kernel_trunc(), RealRugStrictFinite::<PRECISION>::try_from_f64(-3.).unwrap());
    /// ```
    #[inline(always)]
    fn kernel_trunc(self) -> Self {
        self.trunc()
    }
}