astro-float-num 0.3.6

//! Hyperbolic arccosine.

use crate::common::consts::ONE;
use crate::common::util::count_leading_zeroes_skip_first;
use crate::common::util::round_p;
use crate::defs::Error;
use crate::defs::RoundingMode;
use crate::num::BigFloatNumber;
use crate::Consts;
use crate::Sign;
use crate::EXPONENT_MAX;
use crate::WORD_BIT_SIZE;

impl BigFloatNumber {
    /// Computes the hyperbolic arccosine of a number with precision `p`. The result is rounded using the rounding mode `rm`.
    /// This function requires constants cache `cc` for computing the result.
    /// Precision is rounded upwards to the word size.
    ///
    /// ## Errors
    ///
    ///  - ExponentOverflow: the result is too large or too small number.
    ///  - MemoryAllocation: failed to allocate memory.
    ///  - InvalidArgument: when `self` < 1, or the precision is incorrect.
    pub fn acosh(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Result<Self, Error> {
        let p = round_p(p);

        let cmpone = self.cmp(&ONE);
        if cmpone == 0 {
            return Self::new2(p, Sign::Pos, self.inexact());
        } else if cmpone < 0 {
            return Err(Error::InvalidArgument);
        }

        let p_max = p.max(self.mantissa_max_bit_len());

        let mut x = self.clone()?;
        x.set_inexact(false);

        let mut p_inc = WORD_BIT_SIZE;
        let mut p_wrk = p_max + p_inc;

        loop {
            let mut p_x = p_wrk;

            let mut ret = if (self.exponent() as isize - 1) / 2 > p_x as isize + 2 {
                // acosh(x) = ln(2*x)
                if self.exponent() == EXPONENT_MAX {
                    // ln(2) + ln(x)
                    p_x += 2;

                    let lnx = x.ln(p_x, RoundingMode::None, cc)?;

                    let ln2 = cc.ln_2_num(p_x, RoundingMode::None)?;

                    ln2.add(&lnx, p_x, RoundingMode::None)?
                } else {
                    x.set_exponent(self.exponent() + 1);

                    let lnx = x.ln(p_x, RoundingMode::None, cc)?;

                    x.set_exponent(self.exponent());

                    lnx
                }
            } else {
                // ln(x + sqrt(x*x - 1))
                let mut additional_prec = 0;
                if self.exponent() == 1 {
                    additional_prec = count_leading_zeroes_skip_first(self.mantissa().digits());
                }

                p_x += 6 + additional_prec;

                let xx = x.mul(&x, p_x, RoundingMode::None)?;

                let d1 = xx.sub(&ONE, p_x, RoundingMode::None)?;

                let d2 = d1.sqrt(p_x, RoundingMode::None)?;

                let d3 = d2.add(&x, p_x, RoundingMode::FromZero)?;

                d3.ln(p_x, RoundingMode::None, cc)?
            };

            if ret.try_set_precision(p, rm, p_wrk)? {
                ret.set_inexact(ret.inexact() | self.inexact());
                return Ok(ret);
            }

            p_wrk += p_inc;
            p_inc = round_p(p_wrk / 5);
        }
    }
}

#[cfg(test)]
mod tests {

    use crate::Exponent;

    use super::*;

    #[test]
    fn test_acosh() {
        let mut cc = Consts::new().unwrap();
        let rm = RoundingMode::ToEven;
        /* let n1 = BigFloatNumber::from_words(&[144, 9223372036854775808], Sign::Pos, 1).unwrap();
        let n2 = n1.acosh(128, RoundingMode::Down, &mut cc).unwrap();
        println!("{:?}", n2);
        //println!("{:?}", n2.format(crate::Radix::Bin, rm).unwrap());
        return; */

        // near 1
        let p = 448;
        let n1 = BigFloatNumber::parse("1.0000000000000000000000000000000000000000000B56A0EBA6F7D47E21A7B2A7806A698BABAF2F05BC61E2F8FB50FE0B98F55B181AC9C8_e+0", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();
        let n2 = n1.acosh(p, rm, &mut cc).unwrap();
        let n3 = BigFloatNumber::parse("4.C31368910963B1A1BCFC0EDBD393FB7A5E876F9751D93A20E7E48EC0D16090ADA5F46DF2184D32A19C500088EA09CBD4F23DF713113D8A58_e-16", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();

        // println!("{:?}", n1.format(crate::Radix::Bin, rm).unwrap());
        // println!("{:?}", n2.format(crate::Radix::Hex, rm).unwrap());

        assert!(n2.cmp(&n3) == 0);

        // large exp
        let p = 320;
        let n1 = BigFloatNumber::parse("1.921FB54442D18469898CC51701B839A200000000000000004D3C337F7C8D419EBBFC39B4BEC14AF6_e+1000", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();
        let n2 = n1.acosh(p, rm, &mut cc).unwrap();
        let n3 = BigFloatNumber::parse("2.C5DAB0AF9025886C3364C7B6D6741EB19D4FB009D3F92CA21B77498D9F0666363C665F2F324EAEC8_e+3", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();

        //println!("{:?}", n2.format(crate::Radix::Hex, rm).unwrap());

        assert!(n2.cmp(&n3) == 0);

        let d1 = BigFloatNumber::max_value(p).unwrap();
        let d2 = BigFloatNumber::min_value(p).unwrap();
        let d3 = BigFloatNumber::min_positive(p).unwrap();

        // MAX
        let mut d4 = d1.acosh(p, RoundingMode::Down, &mut cc).unwrap();

        d4.set_exponent(d4.exponent() - 1);
        let d5 = d4.cosh(p, rm, &mut cc).unwrap();
        let mut d5 = d5.mul(&d5, p, rm).unwrap();
        d5.set_exponent(d5.exponent() + 1);

        let mut eps = ONE.clone().unwrap();
        eps.set_exponent(
            d1.exponent() - p as Exponent + core::mem::size_of::<Exponent>() as Exponent * 8,
        );

        assert!(
            d1.sub(&d5, p, RoundingMode::ToEven)
                .unwrap()
                .abs()
                .unwrap()
                .cmp(&eps)
                < 0
        );

        assert!(d2.acosh(p, rm, &mut cc).unwrap_err() == Error::InvalidArgument);
        assert!(d3.acosh(p, rm, &mut cc).unwrap_err() == Error::InvalidArgument);

        assert!(ONE.acosh(p, rm, &mut cc).unwrap().is_zero());
    }

    #[ignore]
    #[test]
    #[cfg(feature = "std")]
    fn acosh_perf() {
        let mut cc = Consts::new().unwrap();
        let mut n = vec![];
        let p = 160;
        for _ in 0..10000 {
            n.push(BigFloatNumber::random_normal(p, 0, 5).unwrap());
        }

        for _ in 0..5 {
            let start_time = std::time::Instant::now();
            for ni in n.iter() {
                let _f = ni.acosh(p, RoundingMode::ToEven, &mut cc).unwrap();
            }
            let time = start_time.elapsed();
            println!("{}", time.as_millis());
        }
    }
}