oxinum_complex/native/

transcendental.rs

//! Transcendental functions for native [`BigComplex`]: `abs`, `arg`, `exp`,
//! `ln`, `sqrt`, and `pow`.
//!
//! Every method works at a working precision of `prec + 10` bits internally
//! (the `guard`) and rounds each delivered component to `prec` bits with the
//! caller's [`RoundingMode`]. 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.

use oxinum_core::{OxiNumError, OxiNumResult};
use oxinum_float::native::{BigFloat, RoundingMode};

use super::BigComplex;

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

impl BigComplex {
    /// The magnitude `|z| = sqrt(a² + b²)` as a real [`BigFloat`] at `prec` bits.
    ///
    /// Returns the canonical zero for a zero input (avoiding a `sqrt(0)` round
    /// trip). The squared magnitude is taken from [`BigComplex::norm_sqr`].
    ///
    /// # Errors
    ///
    /// Propagates any error from [`BigFloat::sqrt`] (none expected for the
    /// non-negative `norm_sqr`).
    ///
    /// # Examples
    ///
    /// ```
    /// use oxinum_complex::native::{BigComplex, RoundingMode};
    /// let z = BigComplex::from_f64(3.0, 4.0, 80).expect("finite");
    /// let m = z.abs(80, RoundingMode::HalfEven).expect("abs");
    /// assert!((m.to_f64() - 5.0).abs() < 1e-12);
    /// ```
    pub fn abs(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigFloat> {
        if self.is_zero() {
            return Ok(BigFloat::zero(prec));
        }
        let guard = prec.saturating_add(GUARD);
        let nrm = self.norm_sqr().with_precision(guard, mode);
        Ok(nrm.sqrt(guard, mode)?.with_precision(prec, mode))
    }

    /// The argument `arg(z) = atan2(b, a)` as a real [`BigFloat`] at `prec`
    /// bits, the principal value in `(−π, π]`.
    ///
    /// # Errors
    ///
    /// Propagates any error from [`BigFloat::atan2`].
    ///
    /// # Examples
    ///
    /// ```
    /// use oxinum_complex::native::{BigComplex, RoundingMode};
    /// let i = BigComplex::i(80, RoundingMode::HalfEven);
    /// let a = i.arg(80, RoundingMode::HalfEven).expect("arg");
    /// assert!((a.to_f64() - std::f64::consts::FRAC_PI_2).abs() < 1e-12);
    /// ```
    pub fn arg(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigFloat> {
        let guard = prec.saturating_add(GUARD);
        let re = self.re.clone().with_precision(guard, mode);
        let im = self.im.clone().with_precision(guard, mode);
        Ok(im.atan2(&re, guard, mode)?.with_precision(prec, mode))
    }

    /// The complex exponential `exp(z) = eᵃ·(cos b + i·sin b)` at `prec` bits.
    ///
    /// # Errors
    ///
    /// Propagates errors from [`BigFloat::exp`] (e.g. overflow when `a` is huge)
    /// and from the real trig routines.
    ///
    /// # Examples
    ///
    /// ```
    /// use oxinum_complex::native::{BigComplex, RoundingMode};
    /// // exp(iπ) ≈ −1 + 0i  (Euler's identity).
    /// let z = BigComplex::from_f64(0.0, std::f64::consts::PI, 80).expect("finite");
    /// let e = z.exp(80, RoundingMode::HalfEven).expect("exp");
    /// assert!((e.re().to_f64() + 1.0).abs() < 1e-12);
    /// assert!(e.im().to_f64().abs() < 1e-12);
    /// ```
    pub fn exp(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
        let guard = prec.saturating_add(GUARD);
        let a = self.re.clone().with_precision(guard, mode);
        let b = self.im.clone().with_precision(guard, mode);

        let exp_a = a.exp(guard, mode)?;
        let cos_b = b.cos(guard, mode)?;
        let sin_b = b.sin(guard, mode)?;

        let re = (&exp_a * &cos_b).with_precision(prec, mode);
        let im = (&exp_a * &sin_b).with_precision(prec, mode);
        Ok(BigComplex { re, im })
    }

    /// The principal complex logarithm
    /// `ln(z) = ½·ln(a² + b²) + i·atan2(b, a)` at `prec` bits.
    ///
    /// # Errors
    ///
    /// - [`OxiNumError::Domain`] if `z` is zero (`ln(0)` is undefined).
    /// - Propagates errors from [`BigFloat::ln`] / [`BigFloat::atan2`].
    ///
    /// # Examples
    ///
    /// ```
    /// use oxinum_complex::native::{BigComplex, RoundingMode};
    /// // ln(−1) ≈ 0 + iπ.
    /// let z = BigComplex::from_f64(-1.0, 0.0, 80).expect("finite");
    /// let l = z.ln(80, RoundingMode::HalfEven).expect("ln");
    /// assert!(l.re().to_f64().abs() < 1e-12);
    /// assert!((l.im().to_f64() - std::f64::consts::PI).abs() < 1e-12);
    /// ```
    pub fn ln(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
        if self.is_zero() {
            return Err(OxiNumError::Domain("ln(0) is undefined".into()));
        }
        let guard = prec.saturating_add(GUARD);
        let a = self.re.clone().with_precision(guard, mode);
        let b = self.im.clone().with_precision(guard, mode);

        // re = ½·ln(a² + b²)
        let nrm = self.norm_sqr().with_precision(guard, mode);
        let ln_nrm = nrm.ln(guard, mode)?;
        let half = BigFloat::from_f64(0.5, guard)?;
        let re = (&ln_nrm * &half).with_precision(prec, mode);

        // im = atan2(b, a)
        let im = b.atan2(&a, guard, mode)?.with_precision(prec, mode);

        Ok(BigComplex { re, im })
    }

    /// The principal square root `sqrt(z)` at `prec` bits.
    ///
    /// 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.
    ///
    /// # Errors
    ///
    /// Propagates errors from [`BigFloat::sqrt`] (none expected: radicands are
    /// clamped non-negative).
    ///
    /// # Examples
    ///
    /// ```
    /// use oxinum_complex::native::{BigComplex, RoundingMode};
    /// // sqrt(2i) = 1 + i.
    /// let z = BigComplex::from_f64(0.0, 2.0, 80).expect("finite");
    /// let r = z.sqrt(80, RoundingMode::HalfEven).expect("sqrt");
    /// assert!((r.re().to_f64() - 1.0).abs() < 1e-12);
    /// assert!((r.im().to_f64() - 1.0).abs() < 1e-12);
    /// ```
    pub fn sqrt(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
        if self.is_zero() {
            return Ok(BigComplex::zero(prec));
        }
        let guard = prec.saturating_add(GUARD);

        // Real-axis fast path: b == 0.
        if self.im.is_zero() {
            let a = self.re.clone().with_precision(guard, mode);
            if !a.is_sign_negative() {
                // a >= 0: re = √a, im = 0.
                let re = a.sqrt(guard, mode)?.with_precision(prec, mode);
                return Ok(BigComplex {
                    re,
                    im: BigFloat::zero(prec),
                });
            } else {
                // a < 0: re = 0, im = √(−a).
                let neg_a = (-&a).with_precision(guard, mode);
                let im = neg_a.sqrt(guard, mode)?.with_precision(prec, mode);
                return Ok(BigComplex {
                    re: BigFloat::zero(prec),
                    im,
                });
            }
        }

        // General case.
        let a = self.re.clone().with_precision(guard, mode);
        let two = BigFloat::from_i64(2, guard, mode);

        // m = |z| at guard precision.
        let m = {
            let nrm = self.norm_sqr().with_precision(guard, mode);
            nrm.sqrt(guard, mode)?
        };

        let zero = BigFloat::zero(guard);

        // re = sqrt((m + a) / 2), clamping a tiny-negative radicand up to 0.
        let re_radicand = {
            let s = &m + &a;
            let r = s.div_ref_with_mode(&two, mode)?;
            if r < zero {
                zero.clone()
            } else {
                r
            }
        };
        let re = re_radicand.sqrt(guard, mode)?.with_precision(prec, mode);

        // im_mag = sqrt((m − a) / 2), same clamp.
        let im_radicand = {
            let s = &m - &a;
            let r = s.div_ref_with_mode(&two, mode)?;
            if r < zero {
                zero.clone()
            } else {
                r
            }
        };
        let im_mag = im_radicand.sqrt(guard, mode)?;

        // Apply sign(b): for b < 0 the imaginary part is negative.
        let im = if self.im.is_sign_negative() {
            (-&im_mag).with_precision(prec, mode)
        } else {
            im_mag.with_precision(prec, mode)
        };

        Ok(BigComplex { re, im })
    }

    /// The complex power `z^w = exp(w · ln z)` at `prec` bits.
    ///
    /// 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 [`BigComplex::ln`] / [`BigComplex::exp`].
    ///
    /// # Examples
    ///
    /// ```
    /// use oxinum_complex::native::{BigComplex, RoundingMode};
    /// // i^2 = −1.
    /// let i = BigComplex::i(80, RoundingMode::HalfEven);
    /// let two = BigComplex::from_f64(2.0, 0.0, 80).expect("finite");
    /// let r = i.pow(&two, 80, RoundingMode::HalfEven).expect("pow");
    /// assert!((r.re().to_f64() + 1.0).abs() < 1e-12);
    /// assert!(r.im().to_f64().abs() < 1e-12);
    /// ```
    pub fn pow(&self, w: &BigComplex, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
        if self.is_zero() {
            return if w.is_zero() {
                Ok(BigComplex::one(prec, mode))
            } else {
                Ok(BigComplex::zero(prec))
            };
        }
        let guard = prec.saturating_add(GUARD);
        let ln_z = self.ln(guard, mode)?;
        let prod = w.mul_core(&ln_z);
        prod.exp(prec, mode)
    }
}

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

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

    const PREC: u32 = 80;
    const MODE: RoundingMode = RoundingMode::HalfEven;

    fn c(re: f64, im: f64) -> BigComplex {
        BigComplex::from_f64(re, im, PREC).expect("finite parts")
    }

    #[test]
    fn abs_three_four_is_five() {
        let m = c(3.0, 4.0).abs(PREC, MODE).expect("abs");
        assert!((m.to_f64() - 5.0).abs() < 1e-9, "|3+4i| = {}", m.to_f64());
    }

    #[test]
    fn abs_zero_is_zero() {
        let m = BigComplex::zero(PREC).abs(PREC, MODE).expect("abs");
        assert!(m.is_zero());
    }

    #[test]
    fn arg_of_i_is_half_pi() {
        let a = BigComplex::i(PREC, MODE).arg(PREC, MODE).expect("arg");
        assert!(
            (a.to_f64() - std::f64::consts::FRAC_PI_2).abs() < 1e-9,
            "arg(i) = {}",
            a.to_f64()
        );
    }

    #[test]
    fn exp_i_pi_is_minus_one() {
        // exp(iπ) ≈ −1.
        let z = c(0.0, std::f64::consts::PI);
        let e = z.exp(PREC, MODE).expect("exp");
        assert!(
            (e.re().to_f64() + 1.0).abs() < 1e-9,
            "re = {}",
            e.re().to_f64()
        );
        assert!(e.im().to_f64().abs() < 1e-9, "im = {}", e.im().to_f64());
    }

    #[test]
    fn ln_minus_one_is_i_pi() {
        let l = c(-1.0, 0.0).ln(PREC, MODE).expect("ln");
        assert!(l.re().to_f64().abs() < 1e-9, "re = {}", l.re().to_f64());
        assert!(
            (l.im().to_f64() - std::f64::consts::PI).abs() < 1e-9,
            "im = {}",
            l.im().to_f64()
        );
    }

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

    #[test]
    fn sqrt_minus_one_is_i() {
        let r = c(-1.0, 0.0).sqrt(PREC, MODE).expect("sqrt");
        assert!(r.re().to_f64().abs() < 1e-9, "re = {}", r.re().to_f64());
        assert!(
            (r.im().to_f64() - 1.0).abs() < 1e-9,
            "im = {}",
            r.im().to_f64()
        );
    }

    #[test]
    fn sqrt_two_i_is_one_plus_i() {
        let r = c(0.0, 2.0).sqrt(PREC, MODE).expect("sqrt");
        assert!(
            (r.re().to_f64() - 1.0).abs() < 1e-9,
            "re = {}",
            r.re().to_f64()
        );
        assert!(
            (r.im().to_f64() - 1.0).abs() < 1e-9,
            "im = {}",
            r.im().to_f64()
        );
    }

    #[test]
    fn sqrt_positive_real() {
        // sqrt(4) = 2.
        let r = c(4.0, 0.0).sqrt(PREC, MODE).expect("sqrt");
        assert!((r.re().to_f64() - 2.0).abs() < 1e-9);
        assert!(r.im().to_f64().abs() < 1e-12);
    }

    #[test]
    fn sqrt_squared_roundtrip() {
        // (sqrt(z))² ≈ z for a general z.
        let z = c(2.0, -3.0);
        let r = z.sqrt(PREC, MODE).expect("sqrt");
        let sq = r.mul_core(&r);
        assert!(
            (sq.re().to_f64() - 2.0).abs() < 1e-9,
            "re = {}",
            sq.re().to_f64()
        );
        assert!(
            (sq.im().to_f64() + 3.0).abs() < 1e-9,
            "im = {}",
            sq.im().to_f64()
        );
    }

    #[test]
    fn pow_i_squared_is_minus_one() {
        let r = BigComplex::i(PREC, MODE)
            .pow(&c(2.0, 0.0), PREC, MODE)
            .expect("pow");
        assert!(
            (r.re().to_f64() + 1.0).abs() < 1e-9,
            "re = {}",
            r.re().to_f64()
        );
        assert!(r.im().to_f64().abs() < 1e-9, "im = {}", r.im().to_f64());
    }

    #[test]
    fn pow_zero_zero_is_one() {
        let r = BigComplex::zero(PREC)
            .pow(&BigComplex::zero(PREC), PREC, MODE)
            .expect("pow");
        assert!((r.re().to_f64() - 1.0).abs() < 1e-12);
        assert!(r.im().to_f64().abs() < 1e-12);
    }

    #[test]
    fn pow_zero_base_nonzero_exp_is_zero() {
        let r = BigComplex::zero(PREC)
            .pow(&c(2.0, 1.0), PREC, MODE)
            .expect("pow");
        assert!(r.is_zero());
    }
}