g_math 0.7.0

//! Decimal natural logarithm: `ln(x)` at compute dp.
//!
//! # Algorithm
//!
//! 1. **Domain check**: `x > 0` or `DomainError`.
//! 2. **Range reduction** (halving/doubling): reduce `x` to `m ∈ [1, 2)` by extracting
//!    `k` such that `x = 2^k × m`. Then `ln(x) = k × ln(2) + ln(m)`.
//! 3. **Core series** for `m ∈ [1, 2)`: use `atanh` series with `s = (m-1)/(m+1) ∈ [0, 1/3)`:
//!    `ln(m) = 2 × (s + s³/3 + s⁵/5 + s⁷/7 + ...)`
//!    For `|s| ≤ 1/3`, each term shrinks by factor `1/9`, giving very fast convergence.
//! 4. **`ln(2)` constant**: computed once via the same `atanh` series applied to `1/3`:
//!    `ln(2) = 2 × atanh(1/3)`. Cached thread-local.

#[allow(unused_imports)]
use super::word::DecimalWord;
#[allow(unused_imports)]
use crate::fixed_point::width::WideInt;
#[allow(unused_imports)]
use super::decimal_compute::{decimal_compute_one, decimal_compute_halve, decimal_compute_add, decimal_compute_sub, decimal_compute_mul, decimal_compute_div, decimal_compute_cmp, decimal_compute_neg};
use crate::fixed_point::domains::symbolic::rational::rational_number::OverflowDetected;


/// Compute `ln(2)` at decimal compute dp, cached.
///
/// Formula: `ln(2) = 2 × atanh(1/3) = 2 × (1/3 + (1/3)³/3 + (1/3)⁵/5 + ...)`
fn ln2_at_compute<C: DecimalWord>() -> Result<C, OverflowDetected> {
    struct Ln2<C>(C);
    Ok(super::cache::typed_cached(|| Ok(Ln2(compute_ln2::<C>()?)))?.0)
}

fn compute_ln2<C: DecimalWord>() -> Result<C, OverflowDetected> {

    let one = decimal_compute_one::<C>();
    let three = C::decimal_compute_from_int(3);
    let one_third = decimal_compute_div(one, three)?;

    // atanh(1/3) via direct series (|s| = 1/3, fast convergence)
    let s = one_third;
    let s_sq = decimal_compute_mul(s, s);
    let mut term = s;
    let mut sum = s;

    // atanh(s) = s + s³/3 + s⁵/5 + ...
    // All positive terms.
    let max_terms = (C::DP as u32) + 20;
    for k in 1..=max_terms {
        term = decimal_compute_mul(term, s_sq);
        if C::decimal_compute_is_zero(&term) {
            break;
        }
        let divisor = (2 * k as u64) + 1;
        let contribution = C::decimal_compute_div_int(term, divisor);
        if C::decimal_compute_is_zero(&contribution) {
            break;
        }
        sum = decimal_compute_add(sum, contribution);
    }

    // ln(2) = 2 × atanh(1/3)
    let ln2 = decimal_compute_add(sum, sum);

    Ok(ln2)
}

/// Compute `ln(x)` for `x > 0` at compute dp.
pub fn decimal_ln<C: DecimalWord>(x: C) -> Result<C, OverflowDetected> {
    if C::decimal_compute_is_negative(&x) || C::decimal_compute_is_zero(&x) {
        return Err(OverflowDetected::DomainError);
    }

    let one = decimal_compute_one::<C>();
    let two = C::decimal_compute_from_int(2);

    // Range reduction: find k such that x / 2^k ∈ [1, 2)
    // Halve while ≥ 2, double while < 1
    let mut m = x;
    let mut k: i64 = 0;

    // First: if x is much larger than 2, halve repeatedly
    while decimal_compute_cmp(&m, &two) != std::cmp::Ordering::Less {
        m = decimal_compute_halve(m);
        k += 1;
    }

    // If x < 1, double repeatedly
    while decimal_compute_cmp(&m, &one) == std::cmp::Ordering::Less {
        // Double: m = m + m
        m = decimal_compute_add(m, m);
        k -= 1;
    }

    // Now m ∈ [1, 2). Compute ln(m) via atanh series.
    // s = (m - 1) / (m + 1), |s| ∈ [0, 1/3)
    let m_minus_one = decimal_compute_sub(m, one);
    let m_plus_one = decimal_compute_add(m, one);
    let s = decimal_compute_div(m_minus_one, m_plus_one)?;

    // ln(m) = 2 × atanh(s) = 2 × (s + s³/3 + s⁵/5 + ...)
    let s_sq = decimal_compute_mul(s, s);
    let mut term = s;
    let mut atanh_sum = s;

    let max_terms = (C::DP as u32) + 20;
    for j in 1..=max_terms {
        term = decimal_compute_mul(term, s_sq);
        if C::decimal_compute_is_zero(&term) {
            break;
        }
        let divisor = (2 * j as u64) + 1;
        let contribution = C::decimal_compute_div_int(term, divisor);
        if C::decimal_compute_is_zero(&contribution) {
            break;
        }
        atanh_sum = decimal_compute_add(atanh_sum, contribution);
    }

    let ln_m = decimal_compute_add(atanh_sum, atanh_sum); // 2 × atanh(s)

    // Final: ln(x) = k × ln(2) + ln(m)
    if k == 0 {
        return Ok(ln_m);
    }

    let ln2 = ln2_at_compute::<C>()?;
    let k_ln2 = if k > 0 {
        C::mul_compute_by_int(ln2, k as u64)
    } else {
        decimal_compute_neg(C::mul_compute_by_int(ln2, (-k) as u64))
    };

    Ok(decimal_compute_add(k_ln2, ln_m))
}

#[cfg(test)]
mod tests {
    //! On the 38-digit compute word (I256), in every build.
    use super::*;
    type ComputeStorage = crate::fixed_point::I256;
    #[allow(unused_imports)]
    use super::super::decimal_compute::*;
    use super::super::decimal_compute::pow10_compute_ct;
    use crate::fixed_point::i256::I256;

    fn parse_decimal_str(s: &str) -> I256 {
        let mut result = I256::from_i128(0);
        let ten = I256::from_i128(10);
        for ch in s.chars() {
            let digit = ch.to_digit(10).expect("non-digit");
            result = result * ten + I256::from_i128(digit as i128);
        }
        result
    }

    #[test]
    fn ln_one_is_zero() {
        let result = decimal_ln(decimal_compute_one::<ComputeStorage>()).unwrap();
        assert_eq!(result, decimal_compute_zero::<ComputeStorage>());
    }

    #[test]
    fn ln_negative_errors() {
        let neg_one = decimal_compute_neg(decimal_compute_one::<ComputeStorage>());
        assert!(decimal_ln(neg_one).is_err());
    }

    #[test]
    fn ln_zero_errors() {
        assert!(decimal_ln(decimal_compute_zero::<ComputeStorage>()).is_err());
    }

    /// mpmath: ln(2) = 0.69314718055994530941723212145817656807550013436025...
    #[test]
    fn ln_two_mpmath() {
        let two = decimal_compute_from_int::<ComputeStorage>(2);
        let result = decimal_ln(two).unwrap();
        let expected = parse_decimal_str("69314718055994530941723212145817656808");
        let diff = if result > expected { result - expected } else { expected - result };
        let tolerance = I256::from_i128(1_000_000);
        assert!(
            diff < tolerance,
            "ln(2) precision: got={:?} expected={:?} diff={:?}",
            result, expected, diff
        );
    }

    /// mpmath: ln(10) = 2.30258509299404568401799145468436420760110148862877...
    #[test]
    fn ln_ten_mpmath() {
        let ten = decimal_compute_from_int::<ComputeStorage>(10);
        let result = decimal_ln(ten).unwrap();
        let expected = parse_decimal_str("230258509299404568401799145468436420760");
        let diff = if result > expected { result - expected } else { expected - result };
        let tolerance = I256::from_i128(10_000_000);
        assert!(
            diff < tolerance,
            "ln(10) precision: got={:?} expected={:?} diff={:?}",
            result, expected, diff
        );
    }

    /// mpmath: ln(0.5) = -0.69314718055994530941723212145817656807550013436026...
    #[test]
    fn ln_half_mpmath() {
        let half = pow10_compute_ct::<ComputeStorage>(37) * I256::from_i128(5);
        let result = decimal_ln(half).unwrap();
        let expected_abs = parse_decimal_str("69314718055994530941723212145817656808");
        // result should be -expected_abs
        let result_abs = if decimal_compute_is_negative(&result) {
            decimal_compute_neg(result)
        } else {
            result
        };
        let diff = if result_abs > expected_abs { result_abs - expected_abs } else { expected_abs - result_abs };
        let tolerance = I256::from_i128(1_000_000);
        assert!(
            diff < tolerance,
            "ln(0.5) precision: got={:?} expected_abs={:?}",
            result, expected_abs
        );
        assert!(decimal_compute_is_negative(&result), "ln(0.5) should be negative");
    }
}