astro_float_num/ops/
tan.rs

Help
1//! Tangent.
2
3use crate::common::consts::ONE;
4use crate::common::consts::TRIG_EXP_THRES;
5use crate::common::consts::TWO;
6use crate::common::util::calc_add_cost;
7use crate::common::util::calc_mul_cost;
8use crate::common::util::round_p;
9use crate::defs::Error;
10use crate::defs::RoundingMode;
11use crate::num::BigFloatNumber;
12use crate::ops::consts::Consts;
13use crate::ops::series::series_cost_optimize;
14use crate::ops::series::ArgReductionEstimator;
15use crate::ops::series::PolycoeffGen;
16use crate::ops::util::compute_small_exp;
17use crate::Exponent;
18use crate::EXPONENT_MIN;
19use crate::WORD_BIT_SIZE;
20
21// Polynomial coefficient generator (for tan it only used for cost estmation).
22struct TanPolycoeffGen {
23    iter_cost: usize,
24}
25
26impl TanPolycoeffGen {
27    fn new(p: usize) -> Result<Self, Error> {
28        let iter_cost =
29            9 * calc_mul_cost(p) + 2 * (calc_add_cost(p) + calc_add_cost(WORD_BIT_SIZE));
30
31        Ok(TanPolycoeffGen { iter_cost })
32    }
33}
34
35impl PolycoeffGen for TanPolycoeffGen {
36    fn next(&mut self, _rm: RoundingMode) -> Result<&BigFloatNumber, Error> {
37        Ok(&ONE)
38    }
39
40    #[inline]
41    fn iter_cost(&self) -> usize {
42        self.iter_cost
43    }
44}
45
46struct TanArgReductionEstimator {}
47
48impl ArgReductionEstimator for TanArgReductionEstimator {
49    /// Estimates cost of reduction n times for number with precision p.
50    fn reduction_cost(n: usize, p: usize) -> u64 {
51        let cost_mul = calc_mul_cost(p);
52        let cost_add = calc_add_cost(p);
53        n as u64 * (2 * cost_mul + cost_add) as u64
54    }
55
56    /// Given m, the negative power of 2 of a number, returns the negative power of 2 if reduction is applied n times.
57    #[inline]
58    fn reduction_effect(n: usize, m: isize) -> usize {
59        (n as isize + m) as usize
60    }
61}
62
63impl BigFloatNumber {
64    /// Computes the tangent of a number with precision `p`. The result is rounded using the rounding mode `rm`.
65    /// This function requires constants cache `cc` for computing the result.
66    /// Precision is rounded upwards to the word size.
67    ///
68    /// ## Errors
69    ///
70    ///  - ExponentOverflow: the result is too large or too small number.
71    ///  - MemoryAllocation: failed to allocate memory.
72    ///  - InvalidArgument: the precision is incorrect.
73    pub fn tan(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Result<Self, Error> {
74        let p = round_p(p);
75
76        if self.is_zero() {
77            return Self::new2(p, self.sign(), self.inexact());
78        }
79
80        let mut p_inc = WORD_BIT_SIZE;
81        let mut p_wrk = p.max(self.mantissa_max_bit_len());
82
83        compute_small_exp!(self, self.exponent() as isize * 2 - 1, false, p_wrk, p, rm);
84
85        p_wrk += p_inc;
86
87        let mut add_p = (3 - TRIG_EXP_THRES) as usize;
88        loop {
89            let mut x = self.clone()?;
90
91            let p_x = p_wrk + add_p;
92            x.set_precision(p_x, RoundingMode::None)?;
93
94            x = x.reduce_trig_arg(cc, RoundingMode::None)?;
95
96            let (t, _) = x.trig_arg_pi_proximity(cc, RoundingMode::None)?;
97            if add_p < t {
98                add_p = t;
99            } else {
100                let mut ret = x.tan_series(RoundingMode::None)?;
101
102                if ret.try_set_precision(p, rm, p_wrk)? {
103                    ret.set_inexact(ret.inexact() | self.inexact());
104                    break Ok(ret);
105                }
106
107                p_wrk += p_inc;
108                p_inc = round_p(p_wrk / 5);
109            }
110        }
111    }
112
113    fn tan_series(mut self, rm: RoundingMode) -> Result<Self, Error> {
114        let p = self.mantissa_max_bit_len();
115
116        let polycoeff_gen = TanPolycoeffGen::new(p)?;
117        let (reduction_times, _niter, e_eff) = series_cost_optimize::<TanArgReductionEstimator>(
118            p,
119            &polycoeff_gen,
120            -(self.exponent() as isize),
121            1,
122            true,
123        );
124
125        let add_prec = reduction_times as isize * 4 + 9 - e_eff as isize;
126        let p_arg = p + if add_prec > 0 { add_prec as usize } else { 0 };
127        self.set_precision(p_arg, rm)?;
128
129        let arg_holder;
130        let arg = if reduction_times > 0 {
131            arg_holder = self.tan_arg_reduce(reduction_times)?;
132            &arg_holder
133        } else {
134            &self
135        };
136
137        let ret = Self::tan_series_run(arg, rm)?;
138
139        if reduction_times > 0 {
140            ret.tan_arg_restore(reduction_times, rm)
141        } else {
142            Ok(ret)
143        }
144    }
145
146    /// Tangent series
147    fn tan_series_run(&self, rm: RoundingMode) -> Result<Self, Error> {
148        //  sin + cos series combined
149        // tan(x) = x * (((3! - x^2 * 1!) * 5! + x^4 * 3!) * 7! - ...) / (((2! - x^2 * 1!) * 4! + x^4 * 2!) * 6! - ...) / (3*5*7*...)
150
151        let p = self.mantissa_max_bit_len();
152        let mut xx = self.mul(self, p, rm)?;
153        xx.inv_sign();
154        let mut xxacc = BigFloatNumber::from_word(1, 1)?;
155        let mut fct = BigFloatNumber::from_word(2, 1)?;
156        let mut inc = BigFloatNumber::from_word(2, 1)?;
157        let mut q1 = BigFloatNumber::from_word(1, 1)?;
158        let mut p1 = BigFloatNumber::from_word(1, 1)?;
159        let mut q2 = BigFloatNumber::from_word(1, 1)?;
160        let mut p2 = BigFloatNumber::from_word(1, 1)?;
161
162        while fct.exponent() as isize - (xxacc.exponent() as isize) <= p as isize {
163            // -x^2, +x^4, -x^6, ...
164            xxacc = xxacc.mul(&xx, p, rm)?;
165
166            // cos
167            p1 = p1.mul(&fct, p, rm)?;
168            let n1 = xxacc.mul(&q1, p, rm)?;
169            p1 = p1.add(&n1, p, rm)?;
170
171            q1 = q1.mul(&fct, p, rm)?;
172
173            inc = inc.add(&ONE, inc.mantissa_max_bit_len(), rm)?;
174            if fct.mantissa_max_bit_len() < p {
175                fct = fct.mul_full_prec(&inc)?;
176            } else {
177                fct = fct.mul(&inc, p, rm)?;
178            }
179
180            // sin
181            p2 = p2.mul(&fct, p, rm)?;
182            let n1 = xxacc.mul(&q2, p, rm)?;
183            p2 = p2.add(&n1, p, rm)?;
184
185            q2 = q2.mul(&fct, p, rm)?;
186
187            inc = inc.add(&ONE, inc.mantissa_max_bit_len(), rm)?;
188            if fct.mantissa_max_bit_len() < p {
189                fct = fct.mul_full_prec(&inc)?;
190            } else {
191                fct = fct.mul(&inc, p, rm)?;
192            }
193        }
194
195        let n0 = p2.mul(&q1, p, rm)?;
196        let n1 = n0.mul(self, p, rm)?;
197        let n2 = p1.mul(&q2, p, rm)?;
198
199        let mut ret = n1.div(&n2, p, rm)?;
200        ret.set_inexact(true);
201
202        Ok(ret)
203    }
204
205    // reduce argument n times.
206    fn tan_arg_reduce(&self, n: usize) -> Result<Self, Error> {
207        // tan(3*x) = 3*tan(x) - tan(x)^3 / (1 - 3*tan(x)^2)
208        let mut ret = self.clone()?;
209        let p = ret.mantissa_max_bit_len();
210        if ret.exponent() < EXPONENT_MIN + n as Exponent {
211            ret.set_exponent(EXPONENT_MIN);
212            for _ in 0..n - (ret.exponent() - EXPONENT_MIN) as usize {
213                ret = ret.div(&TWO, p, RoundingMode::FromZero)?;
214            }
215        } else {
216            ret.set_exponent(ret.exponent() - n as Exponent);
217        }
218        Ok(ret)
219    }
220
221    // restore value for the argument reduced n times.
222    fn tan_arg_restore(&self, n: usize, rm: RoundingMode) -> Result<Self, Error> {
223        // tan(2*x) = 2*tan(x) / (1 - tan(x)^2)
224
225        let mut val = self.clone()?;
226        let p = val.mantissa_max_bit_len();
227
228        for _ in 0..n {
229            let val_sq = val.mul(&val, p, rm)?;
230            let q = ONE.sub(&val_sq, p, rm)?;
231            val.set_exponent(val.exponent() + 1);
232            val = val.div(&q, p, rm)?;
233        }
234
235        Ok(val)
236    }
237}
238
239#[cfg(test)]
240mod tests {
241
242    use crate::common::util::random_subnormal;
243
244    use super::*;
245
246    #[test]
247    fn test_tangent() {
248        let p = 320;
249        let mut cc = Consts::new().unwrap();
250        let rm = RoundingMode::ToEven;
251        let mut n1 = BigFloatNumber::from_word(2, p).unwrap();
252        n1.set_exponent(0);
253        let _n2 = n1.tan(p, rm, &mut cc).unwrap();
254        //println!("{:?}", n2.format(crate::Radix::Dec, rm).unwrap());
255
256        // asymptotic & extrema testing
257        let mut half_pi = cc.pi_num(128, RoundingMode::None).unwrap();
258        half_pi.set_exponent(1);
259        half_pi.set_precision(p, RoundingMode::None).unwrap();
260
261        let n2 = half_pi.tan(p, rm, &mut cc).unwrap();
262        let n3 = BigFloatNumber::parse(
263            "3.1F0B46DCBD63D29899ECF829DA54DE0EE0852B2569B572B793E50817CEF4C77D959712B45E2B7E4C_e+20",
264            crate::Radix::Hex,
265            p,
266            RoundingMode::None, &mut cc,
267        )
268        .unwrap();
269
270        assert!(n2.cmp(&n3) == 0);
271
272        // large exponent
273        half_pi.set_exponent(256);
274        let n2 = half_pi.tan(p, rm, &mut cc).unwrap();
275        let n3 = BigFloatNumber::parse("4.ECDEC5EF3A1EA5339A46BC0C490F52A86A033C56BCDD413E36C657EB7757F073500B013B9A7B43C0_e+0", crate::Radix::Hex, p, RoundingMode::None, &mut cc).unwrap();
276
277        assert!(n2.cmp(&n3) == 0);
278
279        let d3 = BigFloatNumber::min_positive(p).unwrap();
280        let zero = BigFloatNumber::new(1).unwrap();
281        let n1 = random_subnormal(p);
282
283        assert!(d3.tan(p, rm, &mut cc).unwrap().cmp(&d3) == 0);
284        assert!(zero.tan(p, rm, &mut cc).unwrap().is_zero());
285        assert!(n1.tan(p, rm, &mut cc).unwrap().cmp(&n1) == 0);
286    }
287
288    #[ignore]
289    #[test]
290    #[cfg(feature = "std")]
291    fn tan_perf() {
292        let p = 160;
293        let mut cc = Consts::new().unwrap();
294        let mut n = vec![];
295        for _ in 0..10000 {
296            n.push(BigFloatNumber::random_normal(p, 0, 5).unwrap());
297        }
298
299        for _ in 0..5 {
300            let start_time = std::time::Instant::now();
301            for ni in n.iter() {
302                let _f = ni.tan(p, RoundingMode::ToEven, &mut cc).unwrap();
303            }
304            let time = start_time.elapsed();
305            println!("{}", time.as_millis());
306        }
307    }
308}
astro_float_num/ops/tan.rs

astro_float_num/ops/
tan.rs
