oxinum_complex/

transcendental.rs

//! Transcendental functions for [`crate::CBig`]: `abs`, `arg`, `exp`, `ln`,
//! `sqrt`, and `pow`.
//!
//! Every method works at an internal *guard* precision of `precision + 10`
//! significant decimal digits and delivers each component at the guard
//! precision returned by the underlying [`oxinum_float`] free functions. With
//! `z = a + b·i` (so `a = self.re()`, `b = self.im()`):
//!
//! ```text
//! |z|     = sqrt(a² + b²)              (real, non-negative)
//! arg(z)  = atan2(b, a)               (real, principal value in (−π, π])
//! exp(z)  = eᵃ·(cos b + i·sin b)
//! ln(z)   = ½·ln(a² + b²) + i·atan2(b, a)
//! sqrt(z) = principal branch (see below)
//! z^w     = exp(w · ln z)
//! ```
//!
//! The principal `sqrt` uses the magnitude `m = |z|`:
//!
//! ```text
//! re = sqrt((m + a) / 2)
//! im = sign(b) · sqrt((m − a) / 2)
//! ```
//!
//! with the radicands clamped up to `0` before the real `sqrt` to absorb the
//! tiny negative values rounding can produce, and the purely-real input
//! (`b == 0`) handled by an exact axis split (so `sqrt(-1) = +i`).
//!
//! This is the decimal-backed ([`DBig`]) mirror of the native binary
//! implementation in [`crate::native`]; the formulas and branch cuts match
//! exactly.

use core::str::FromStr;

use oxinum_float::{atan2, cos, exp, ln, sin, sqrt};

use crate::{CBig, DBig, OxiNumError, OxiNumResult};

/// Working-precision headroom added on top of the requested `precision`.
const GUARD: usize = 10;

/// Parse a decimal literal into a [`DBig`], mapping any failure to
/// [`OxiNumError::Parse`].
fn make_dbig(s: &str) -> OxiNumResult<DBig> {
    DBig::from_str(s).map_err(|e| OxiNumError::Parse(format!("{e}").into()))
}

impl CBig {
    /// The magnitude `|z| = sqrt(a² + b²)` as a real [`DBig`] at `precision`
    /// significant digits.
    ///
    /// Returns an exact decimal zero for a zero input (avoiding a `sqrt(0)`
    /// round trip). The squared magnitude is taken from [`CBig::norm_sqr`].
    ///
    /// # Errors
    ///
    /// Propagates any error from [`oxinum_float::sqrt`] (none expected for the
    /// non-negative `norm_sqr`; an [`OxiNumError::Precision`] is returned if
    /// `precision == 0`).
    pub fn abs(&self, precision: usize) -> OxiNumResult<DBig> {
        if self.is_zero() {
            return make_dbig("0.0");
        }
        sqrt(&self.norm_sqr(), precision)
    }

    /// The argument `arg(z) = atan2(b, a)` as a real [`DBig`] at `precision`
    /// significant digits, the principal value in `(−π, π]`.
    ///
    /// # Errors
    ///
    /// Propagates any error from [`oxinum_float::atan2`].
    pub fn arg(&self, precision: usize) -> OxiNumResult<DBig> {
        atan2(&self.im, &self.re, precision)
    }

    /// The complex exponential `exp(z) = eᵃ·(cos b + i·sin b)` at `precision`
    /// significant digits.
    ///
    /// # Errors
    ///
    /// Propagates errors from [`oxinum_float::exp`], [`oxinum_float::cos`], and
    /// [`oxinum_float::sin`].
    pub fn exp(&self, precision: usize) -> OxiNumResult<CBig> {
        let guard = precision + GUARD;
        let ea = exp(&self.re, guard)?;
        let cosb = cos(&self.im, guard)?;
        let sinb = sin(&self.im, guard)?;
        Ok(CBig::from_parts(&ea * &cosb, &ea * &sinb))
    }

    /// The principal complex logarithm
    /// `ln(z) = ½·ln(a² + b²) + i·atan2(b, a)` at `precision` significant
    /// digits.
    ///
    /// Using `½·ln(norm_sqr)` for the real part avoids an extra `sqrt`.
    ///
    /// # Errors
    ///
    /// - [`OxiNumError::Domain`] if `z` is zero (`ln(0)` is undefined).
    /// - Propagates errors from [`oxinum_float::ln`] / [`oxinum_float::atan2`].
    pub fn ln(&self, precision: usize) -> OxiNumResult<CBig> {
        if self.is_zero() {
            return Err(OxiNumError::Domain("ln(0) is undefined".into()));
        }
        let guard = precision + GUARD;
        let ns = self.norm_sqr();
        let re = &make_dbig("0.5")? * &ln(&ns, guard)?;
        let im = atan2(&self.im, &self.re, guard)?;
        Ok(CBig::from_parts(re, im))
    }

    /// The principal square root `sqrt(z)` at `precision` significant digits.
    ///
    /// Uses `re = sqrt((|z| + a)/2)`, `im = sign(b)·sqrt((|z| − a)/2)`, with the
    /// purely-real input handled by an exact axis split and the radicands
    /// clamped up to zero before the real `sqrt` to absorb rounding noise. The
    /// branch chosen has `re ≥ 0` and matches the IEEE-754 / `num-complex`
    /// principal value (so `sqrt(-1) = +i`).
    ///
    /// # Errors
    ///
    /// Propagates errors from [`oxinum_float::sqrt`] (none expected: radicands
    /// are clamped non-negative).
    pub fn sqrt(&self, precision: usize) -> OxiNumResult<CBig> {
        let zero = DBig::from(0u32);

        if self.is_zero() {
            return Ok(CBig::zero());
        }

        // Real-axis fast path: b == 0.
        if self.im == zero {
            if self.re >= zero {
                // a >= 0: re = √a, im = 0.
                return Ok(CBig::from_parts(sqrt(&self.re, precision)?, zero));
            }
            // a < 0: re = 0, im = √(−a)  (so sqrt(-1) = +i).
            let neg_a = -&self.re;
            return Ok(CBig::from_parts(zero, sqrt(&neg_a, precision)?));
        }

        // General case.
        let guard = precision + GUARD;
        let m = sqrt(&self.norm_sqr(), guard)?; // m = |z|
        let two = make_dbig("2.0")?;

        // re = sqrt((m + a) / 2), clamping a tiny-negative radicand up to 0.
        let mut r_re = &(&m + &self.re) / &two;
        if r_re < zero {
            r_re = DBig::from(0u32);
        }

        // im_mag = sqrt((m − a) / 2), same clamp.
        let mut r_im = &(&m - &self.re) / &two;
        if r_im < zero {
            r_im = DBig::from(0u32);
        }

        let re = sqrt(&r_re, precision)?;
        let im_mag = sqrt(&r_im, precision)?;

        // Apply sign(b): for b < 0 the imaginary part is negative.
        let im = if self.im < zero { -im_mag } else { im_mag };

        Ok(CBig::from_parts(re, im))
    }

    /// The complex power `z^w = exp(w · ln z)` at `precision` significant
    /// digits.
    ///
    /// The zero base is handled by convention: `0^0 = 1` and `0^w = 0` for any
    /// other `w` (avoiding `ln(0)`).
    ///
    /// # Errors
    ///
    /// Propagates errors from [`CBig::ln`] / [`CBig::exp`].
    pub fn pow(&self, w: &CBig, precision: usize) -> OxiNumResult<CBig> {
        if self.is_zero() {
            return if w.is_zero() {
                Ok(CBig::one())
            } else {
                Ok(CBig::zero())
            };
        }
        let guard = precision + GUARD;
        let lz = self.ln(guard)?;
        let prod = w * &lz;
        prod.exp(precision)
    }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use oxinum_float::compute_pi;

    /// Build π at the given decimal precision.
    fn pi(precision: usize) -> DBig {
        compute_pi(precision)
    }

    #[test]
    fn exp_i_pi_is_minus_one() {
        // exp(iπ) ≈ −1 + 0i  (Euler's identity).
        let z = CBig::from_parts(DBig::from(0u32), pi(50));
        let r = z.exp(40).expect("exp");
        let (re, im) = r.to_f64_parts();
        assert!((re + 1.0).abs() < 1e-12, "re = {re}");
        assert!(im.abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn ln_minus_one_is_i_pi() {
        let z = CBig::from_real(make_dbig("-1.0").expect("literal"));
        let r = z.ln(40).expect("ln");
        let (re, im) = r.to_f64_parts();
        assert!(re.abs() < 1e-12, "re = {re}");
        assert!((im - std::f64::consts::PI).abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn ln_zero_is_domain_error() {
        let r = CBig::zero().ln(40);
        assert!(matches!(r, Err(OxiNumError::Domain(_))), "got {r:?}");
    }

    #[test]
    fn sqrt_minus_one_is_i() {
        let z = CBig::from_real(make_dbig("-1.0").expect("literal"));
        let r = z.sqrt(40).expect("sqrt");
        let (re, im) = r.to_f64_parts();
        assert!(re.abs() < 1e-12, "re = {re}");
        assert!((im - 1.0).abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn sqrt_two_i_is_one_plus_i() {
        // sqrt(2i) = 1 + i.
        let z = CBig::from_parts(DBig::from(0u32), DBig::from(2u32));
        let r = z.sqrt(40).expect("sqrt");
        let (re, im) = r.to_f64_parts();
        assert!((re - 1.0).abs() < 1e-12, "re = {re}");
        assert!((im - 1.0).abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn sqrt_positive_real() {
        // sqrt(4) = 2.
        let z = CBig::from_real(DBig::from(4u32));
        let r = z.sqrt(40).expect("sqrt");
        let (re, im) = r.to_f64_parts();
        assert!((re - 2.0).abs() < 1e-12, "re = {re}");
        assert!(im.abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn abs_three_four_is_five() {
        let z = CBig::from_parts(DBig::from(3u32), DBig::from(4u32));
        let m = z.abs(40).expect("abs");
        assert!(m.to_string().starts_with('5'), "|3+4i| = {m}");
    }

    #[test]
    fn abs_zero_is_zero() {
        let m = CBig::zero().abs(40).expect("abs");
        assert_eq!(m, DBig::from(0u32));
    }

    #[test]
    fn arg_of_i_is_half_pi() {
        let a = CBig::i().arg(40).expect("arg");
        let v = a.to_f64();
        assert!(
            (v.value() - std::f64::consts::FRAC_PI_2).abs() < 1e-12,
            "arg(i) = {a}"
        );
    }

    #[test]
    fn pow_i_squared_is_minus_one() {
        // i² = −1.
        let r = CBig::i()
            .pow(&CBig::from_real(DBig::from(2u32)), 40)
            .expect("pow");
        let (re, im) = r.to_f64_parts();
        assert!((re + 1.0).abs() < 1e-12, "re = {re}");
        assert!(im.abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn pow_to_one_is_identity() {
        // z^1 ≈ z on the real axis (arg = 0), where the round trip through
        // ln/exp is exact and independent of the underlying `atan2`.
        let z = CBig::from_real(make_dbig("2.0").expect("literal"));
        let r = z.pow(&CBig::one(), 40).expect("pow");
        let (re, im) = r.to_f64_parts();
        assert!((re - 2.0).abs() < 1e-12, "re = {re}");
        assert!(im.abs() < 1e-12, "im = {im}");
    }

    #[test]
    fn pow_zero_zero_is_one() {
        let r = CBig::zero().pow(&CBig::zero(), 40).expect("pow");
        assert!(r == CBig::one());
    }

    #[test]
    fn pow_zero_base_nonzero_exp_is_zero() {
        let r = CBig::zero()
            .pow(&CBig::from_parts(DBig::from(2u32), DBig::from(1u32)), 40)
            .expect("pow");
        assert!(r.is_zero());
    }
}