astro_float_num/ops/
acos.rs

1//! Arccosine.
2
3use crate::common::consts::ONE;
4use crate::common::util::round_p;
5use crate::defs::Error;
6use crate::defs::RoundingMode;
7use crate::num::BigFloatNumber;
8use crate::ops::consts::Consts;
9use crate::Exponent;
10use crate::Sign;
11use crate::WORD_BIT_SIZE;
12
13const ACOS_EXP_THRES: Exponent = -32;
14
15impl BigFloatNumber {
16    /// Computes the arccosine of a number with precision `p`. The result is rounded using the rounding mode `rm`.
17    /// This function requires constants cache `cc` for computing the result.
18    /// Precision is rounded upwards to the word size.
19    ///
20    /// ## Errors
21    ///
22    ///  - InvalidArgument: argument is greater than 1 or smaller than -1, or the precision is incorrect.
23    ///  - MemoryAllocation: failed to allocate memory.
24    pub fn acos(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Result<Self, Error> {
25        let p = round_p(p);
26
27        let cmpone = self.abs_cmp(&ONE);
28        if cmpone == 0 && self.is_positive() {
29            return Self::new2(p, Sign::Pos, self.inexact());
30        } else if cmpone > 0 {
31            return Err(Error::InvalidArgument);
32        }
33
34        let mut p_inc = WORD_BIT_SIZE;
35        let mut p_wrk = p.max(self.mantissa_max_bit_len()) + p_inc;
36
37        let mut add_p = (1 - ACOS_EXP_THRES) as usize;
38
39        let mut x = self.clone()?;
40        x.set_inexact(false);
41
42        loop {
43            let p_x = p_wrk + add_p;
44            x.set_precision(p_x, RoundingMode::None)?;
45
46            let mut ret = x.asin(p_x, RoundingMode::None, cc)?;
47
48            let mut pi = cc.pi_num(p_x, RoundingMode::None)?;
49            pi.set_exponent(pi.exponent() - 1);
50
51            let ret2 = pi.add(&ret, p_x, RoundingMode::None)?;
52            ret = pi.sub(&ret, p_x, RoundingMode::None)?;
53
54            let t = ret
55                .exponent()
56                .unsigned_abs()
57                .max(ret2.exponent().unsigned_abs()) as usize
58                + 1; // ret near pi / 2 gives cancellation
59            if add_p < t {
60                add_p = t;
61            } else {
62                if ret.try_set_precision(p, rm, p_wrk)? {
63                    ret.set_inexact(ret.inexact() | self.inexact());
64                    return Ok(ret);
65                }
66
67                p_wrk += p_inc;
68                p_inc = round_p(p_wrk / 5);
69            }
70        }
71    }
72}
73
74#[cfg(test)]
75mod tests {
76
77    use crate::common::{consts::ONE, util::random_subnormal};
78
79    use super::*;
80
81    #[test]
82    fn test_arccosine() {
83        let mut cc = Consts::new().unwrap();
84
85        let rm = RoundingMode::ToEven;
86        let p = 64;
87        let mut n1 = BigFloatNumber::from_word(4294967295, p).unwrap();
88        n1.set_exponent(0);
89        //println!("{}", n1.format(crate::Radix::Dec, RoundingMode::None).unwrap());
90        let _n2 = n1.acos(p, rm, &mut cc).unwrap();
91        //println!("{:?}", n2.format(crate::Radix::Dec, rm).unwrap());
92
93        // near 1
94        let p = 320;
95        let n1 = BigFloatNumber::parse(
96            "F.FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF2DC85F7E77EC487_e-1",
97            crate::Radix::Hex,
98            p,
99            RoundingMode::None,
100            &mut cc,
101        )
102        .unwrap();
103        let n2 = n1.acos(p, rm, &mut cc).unwrap();
104        let n3 = BigFloatNumber::parse("5.2049C1114CF98E7B6DB49CCF999F4A5E697D73E5DA6BEC6578098357460BAFFB0C25779F1C63E8D8_e-21", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();
105
106        // println!("{:?}", n2.format(crate::Radix::Hex, rm).unwrap());
107
108        assert!(n2.cmp(&n3) == 0);
109
110        // near zero
111        let n1 = BigFloatNumber::parse("1.921FB54442D18469898CC51701B839A200000000000000004D3C337F7C8D419EBBFC39B4BEC14AF6_e-10", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();
112        let n2 = n1.acos(p, rm, &mut cc).unwrap();
113        let n3 = BigFloatNumber::parse("1.921FB54442D18467F76D0FD2BEE6B538C877D6FA13175F455EB9961B4834A04EF790AE0274E87FF6_e+0", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();
114
115        // println!("{:?}", n2.format(crate::Radix::Hex, rm).unwrap());
116
117        assert!(n2.cmp(&n3) == 0);
118
119        let d1 = BigFloatNumber::max_value(p).unwrap();
120        let d2 = BigFloatNumber::min_value(p).unwrap();
121        let d3 = BigFloatNumber::min_positive(p).unwrap();
122        let zero = BigFloatNumber::new(1).unwrap();
123
124        let mut half_pi = cc.pi_num(p, RoundingMode::ToEven).unwrap();
125        half_pi.set_exponent(1);
126
127        assert!(d1.acos(p, rm, &mut cc).is_err());
128        assert!(d2.acos(p, rm, &mut cc).is_err());
129        assert!(d3.acos(p, rm, &mut cc).unwrap().cmp(&half_pi) == 0);
130        assert!(zero.acos(p, rm, &mut cc).unwrap().cmp(&half_pi) == 0);
131        assert!(ONE.acos(p, rm, &mut cc).unwrap().is_zero());
132
133        // subnormal arg
134        let n1 = random_subnormal(p);
135        assert!(n1.acos(p, rm, &mut cc).unwrap().cmp(&half_pi) == 0);
136    }
137
138    #[ignore]
139    #[test]
140    #[cfg(feature = "std")]
141    fn arccosine_perf() {
142        let mut cc = Consts::new().unwrap();
143        let p = 133;
144
145        let mut n = vec![];
146        for _ in 0..10000 {
147            n.push(BigFloatNumber::random_normal(p, -5, 5).unwrap());
148        }
149
150        for _ in 0..5 {
151            let start_time = std::time::Instant::now();
152            for ni in n.iter() {
153                let _f = ni.acos(p, RoundingMode::ToEven, &mut cc).unwrap();
154            }
155            let time = start_time.elapsed();
156            println!("{}", time.as_millis());
157        }
158    }
159}
astro_float_num/ops/acos.rs

astro_float_num/ops/
acos.rs
