rival/interval/
trig.rs

1use super::value::{Endpoint, Ival, IvalClass, classify};
2use crate::{
3    interval::core::endpoint_unary,
4    mpfr::{
5        mpfr_cos, mpfr_cosu, mpfr_div, mpfr_floor_inplace, mpfr_get_exp, mpfr_pi,
6        mpfr_round_inplace, mpfr_sin, mpfr_sinu, mpfr_tan, mpfr_tanu, zero,
7    },
8};
9use rug::{Assign, Float, float::Round};
10
11const RANGE_REDUCE_PRECISION_CAP: u32 = 1 << 20;
12
13#[derive(Clone, Copy, Debug, PartialEq, Eq)]
14enum PeriodClass {
15    TooWide,
16    NearZero,
17    RangeReduce,
18}
19
20fn classify_ival_periodic(x: &Ival, period_quarter_bitlen: i64) -> PeriodClass {
21    let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
22
23    if xlo.is_infinite() || xhi.is_infinite() {
24        return PeriodClass::TooWide;
25    }
26
27    let (lo_exp, hi_exp) = (mpfr_get_exp(xlo), mpfr_get_exp(xhi));
28    if lo_exp < period_quarter_bitlen && hi_exp < period_quarter_bitlen {
29        return PeriodClass::NearZero;
30    }
31
32    let lo_ulp = lo_exp.saturating_sub(xlo.prec() as i64);
33    let hi_ulp = hi_exp.saturating_sub(xhi.prec() as i64);
34
35    if lo_ulp > 0 || hi_ulp > 0 {
36        if xlo == xhi {
37            PeriodClass::RangeReduce
38        } else {
39            PeriodClass::TooWide
40        }
41    } else {
42        PeriodClass::RangeReduce
43    }
44}
45
46fn range_reduce_precision(xlo: &Float, xhi: &Float, curr_prec: u32) -> u32 {
47    let lo = (mpfr_get_exp(xlo) + 2 * (xlo.prec() as i64)).max(curr_prec as i64) as u32;
48    let hi = (mpfr_get_exp(xhi) + 2 * (xhi.prec() as i64)).max(curr_prec as i64) as u32;
49    lo.max(hi).min(RANGE_REDUCE_PRECISION_CAP).max(curr_prec)
50}
51
52fn ival_div_pi(x: &Ival, prec: u32, round_fn: fn(&mut Float) -> bool) -> (Float, Float) {
53    let (mut pi_lo, mut pi_hi) = (zero(prec), zero(prec));
54    mpfr_pi(&mut pi_lo, Round::Down);
55    mpfr_pi(&mut pi_hi, Round::Up);
56
57    let (mut q_lo, mut q_hi) = (zero(prec), zero(prec));
58    mpfr_div(x.lo.as_float(), &pi_hi, &mut q_lo, Round::Down);
59    mpfr_div(x.hi.as_float(), &pi_lo, &mut q_hi, Round::Up);
60    round_fn(&mut q_lo);
61    round_fn(&mut q_hi);
62    (q_lo, q_hi)
63}
64
65fn ival_div_n_half(
66    x: &Ival,
67    n: u64,
68    prec: u32,
69    round_fn: fn(&mut Float) -> bool,
70) -> (Float, Float) {
71    let n_half = Float::with_val(prec, n) / 2u32;
72
73    let (mut q_lo, mut q_hi) = (zero(prec), zero(prec));
74    mpfr_div(x.lo.as_float(), &n_half, &mut q_lo, Round::Down);
75    mpfr_div(x.hi.as_float(), &n_half, &mut q_hi, Round::Up);
76    round_fn(&mut q_lo);
77    round_fn(&mut q_hi);
78    (q_lo, q_hi)
79}
80
81fn bfeven(x: &Float) -> bool {
82    let prec = x.prec();
83    let mut t = Float::with_val(prec, x);
84    t /= 2;
85    mpfr_floor_inplace(&mut t);
86    t *= 2;
87    t == *x
88}
89
90#[inline]
91fn bfodd(x: &Float) -> bool {
92    !bfeven(x)
93}
94
95fn bfsub_is_one(a: &Float, b: &Float) -> bool {
96    let prec = a.prec().max(b.prec());
97    let mut d = Float::with_val(prec, b);
98    d -= a;
99    d == 1
100}
101
102fn period_quarter_bitlen(n: u64, divisor: u64) -> i64 {
103    let quarter = n / divisor;
104    if quarter > 0 {
105        (quarter.ilog2() + 1) as i64
106    } else {
107        0
108    }
109}
110
111fn endpoint_min<F>(f: &F, lo: &Endpoint, hi: &Endpoint, out: &mut Float) -> bool
112where
113    F: Fn(&Float, &mut Float, Round) -> bool,
114{
115    let prec = out.prec();
116    let mut tmp = zero(prec);
117
118    let imm_lo = endpoint_unary(f, lo, out, Round::Down);
119    let imm_hi = endpoint_unary(f, hi, &mut tmp, Round::Down);
120
121    if tmp < *out {
122        out.assign(&tmp);
123        imm_hi
124    } else if tmp == *out {
125        imm_lo || imm_hi
126    } else {
127        imm_lo
128    }
129}
130
131fn endpoint_max<F>(f: &F, lo: &Endpoint, hi: &Endpoint, out: &mut Float) -> bool
132where
133    F: Fn(&Float, &mut Float, Round) -> bool,
134{
135    let prec = out.prec();
136    let mut tmp = zero(prec);
137
138    let imm_lo = endpoint_unary(f, lo, out, Round::Up);
139    let imm_hi = endpoint_unary(f, hi, &mut tmp, Round::Up);
140
141    if tmp > *out {
142        out.assign(&tmp);
143        imm_hi
144    } else if tmp == *out {
145        imm_lo || imm_hi
146    } else {
147        imm_lo
148    }
149}
150
151impl Ival {
152    /// Compute the interval cosine of `x`.
153    pub fn cos_assign(&mut self, x: &Ival) {
154        self.err = x.err;
155        let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
156
157        match classify_ival_periodic(x, 1) {
158            PeriodClass::TooWide => {
159                self.lo.as_float_mut().assign(-1);
160                self.hi.as_float_mut().assign(1);
161                self.lo.immovable = false;
162                self.hi.immovable = false;
163            }
164            PeriodClass::NearZero => match classify(x, false) {
165                IvalClass::Neg => self.monotonic_assign(&mpfr_cos, x),
166                IvalClass::Pos => self.comonotonic_assign(&mpfr_cos, x),
167                IvalClass::Mix => {
168                    self.set_prec(x.prec());
169                    self.lo.immovable =
170                        endpoint_min(&mpfr_cos, &x.lo, &x.hi, self.lo.as_float_mut());
171                    self.hi.as_float_mut().assign(1);
172                    self.hi.immovable = false;
173                }
174            },
175            PeriodClass::RangeReduce => {
176                let prec = range_reduce_precision(xlo, xhi, x.prec());
177                let (a, b) = ival_div_pi(x, prec, mpfr_floor_inplace);
178
179                if a == b && bfeven(&a) {
180                    self.comonotonic_assign(&mpfr_cos, x);
181                } else if a == b && bfodd(&a) {
182                    self.monotonic_assign(&mpfr_cos, x);
183                } else if bfsub_is_one(&a, &b) && bfeven(&a) {
184                    self.set_prec(x.prec());
185                    self.lo.as_float_mut().assign(-1);
186                    self.lo.immovable = false;
187                    self.hi.immovable =
188                        endpoint_max(&mpfr_cos, &x.lo, &x.hi, self.hi.as_float_mut());
189                } else if bfsub_is_one(&a, &b) && bfodd(&a) {
190                    self.set_prec(x.prec());
191                    self.lo.immovable =
192                        endpoint_min(&mpfr_cos, &x.lo, &x.hi, self.lo.as_float_mut());
193                    self.hi.as_float_mut().assign(1);
194                    self.hi.immovable = false;
195                } else {
196                    self.lo.as_float_mut().assign(-1);
197                    self.hi.as_float_mut().assign(1);
198                    self.lo.immovable = false;
199                    self.hi.immovable = false;
200                }
201            }
202        }
203    }
204
205    /// Compute the interval sine of `x`.
206    pub fn sin_assign(&mut self, x: &Ival) {
207        self.err = x.err;
208        let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
209
210        match classify_ival_periodic(x, 1) {
211            PeriodClass::TooWide => {
212                self.lo.as_float_mut().assign(-1);
213                self.hi.as_float_mut().assign(1);
214                self.lo.immovable = false;
215                self.hi.immovable = false;
216            }
217            PeriodClass::NearZero => {
218                self.monotonic_assign(&mpfr_sin, x);
219            }
220            PeriodClass::RangeReduce => {
221                let prec = range_reduce_precision(xlo, xhi, x.prec());
222                let (a, b) = ival_div_pi(x, prec, mpfr_round_inplace);
223
224                if a == b && bfeven(&a) {
225                    self.monotonic_assign(&mpfr_sin, x);
226                } else if a == b && bfodd(&a) {
227                    self.comonotonic_assign(&mpfr_sin, x);
228                } else if bfsub_is_one(&a, &b) && bfodd(&a) {
229                    self.set_prec(x.prec());
230                    self.lo.as_float_mut().assign(-1);
231                    self.lo.immovable = false;
232                    self.hi.immovable =
233                        endpoint_max(&mpfr_sin, &x.lo, &x.hi, self.hi.as_float_mut());
234                } else if bfsub_is_one(&a, &b) && bfeven(&a) {
235                    self.set_prec(x.prec());
236                    self.lo.immovable =
237                        endpoint_min(&mpfr_sin, &x.lo, &x.hi, self.lo.as_float_mut());
238                    self.hi.as_float_mut().assign(1);
239                    self.hi.immovable = false;
240                } else {
241                    self.lo.as_float_mut().assign(-1);
242                    self.hi.as_float_mut().assign(1);
243                    self.lo.immovable = false;
244                    self.hi.immovable = false;
245                }
246            }
247        }
248    }
249
250    /// Compute the interval tangent of `x`.
251    pub fn tan_assign(&mut self, x: &Ival) {
252        let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
253        let immovable = x.lo.immovable && x.hi.immovable;
254
255        match classify_ival_periodic(x, 0) {
256            PeriodClass::TooWide => {
257                self.lo.as_float_mut().assign(f64::NEG_INFINITY);
258                self.hi.as_float_mut().assign(f64::INFINITY);
259                self.lo.immovable = immovable;
260                self.hi.immovable = immovable;
261                self.err.partial = true;
262                self.err.total = x.err.total;
263            }
264            PeriodClass::NearZero => {
265                self.monotonic_assign(&mpfr_tan, x);
266                self.err = x.err;
267            }
268            PeriodClass::RangeReduce => {
269                let prec = range_reduce_precision(xlo, xhi, x.prec());
270                let (a, b) = ival_div_pi(x, prec, mpfr_round_inplace);
271
272                if a == b {
273                    self.monotonic_assign(&mpfr_tan, x);
274                    self.err = x.err;
275                } else {
276                    self.lo.as_float_mut().assign(f64::NEG_INFINITY);
277                    self.hi.as_float_mut().assign(f64::INFINITY);
278                    self.lo.immovable = immovable;
279                    self.hi.immovable = immovable;
280                    self.err.partial = true;
281                    self.err.total = x.err.total;
282                }
283            }
284        }
285    }
286
287    pub(crate) fn cosu_assign(&mut self, x: &Ival, n: u64) {
288        self.err = x.err;
289        let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
290        let period_qtr = period_quarter_bitlen(n, 4);
291        let cosu = |x: &Float, out: &mut Float, rnd: Round| mpfr_cosu(x, n, out, rnd);
292
293        match classify_ival_periodic(x, period_qtr) {
294            PeriodClass::TooWide => {
295                self.lo.as_float_mut().assign(-1);
296                self.hi.as_float_mut().assign(1);
297                self.lo.immovable = false;
298                self.hi.immovable = false;
299            }
300            PeriodClass::NearZero => match classify(x, false) {
301                IvalClass::Neg => self.monotonic_assign(&cosu, x),
302                IvalClass::Pos => self.comonotonic_assign(&cosu, x),
303                IvalClass::Mix => {
304                    self.set_prec(x.prec());
305                    self.lo.immovable = endpoint_min(&cosu, &x.lo, &x.hi, self.lo.as_float_mut());
306                    self.hi.as_float_mut().assign(1);
307                    self.hi.immovable = false;
308                }
309            },
310            PeriodClass::RangeReduce => {
311                let prec = range_reduce_precision(xlo, xhi, x.prec());
312                let (a, b) = ival_div_n_half(x, n, prec, mpfr_floor_inplace);
313
314                if a == b && bfeven(&a) {
315                    self.comonotonic_assign(&cosu, x);
316                } else if a == b && bfodd(&a) {
317                    self.monotonic_assign(&cosu, x);
318                } else if bfsub_is_one(&a, &b) && bfeven(&a) {
319                    self.set_prec(x.prec());
320                    self.lo.as_float_mut().assign(-1);
321                    self.lo.immovable = false;
322                    self.hi.immovable = endpoint_max(&cosu, &x.lo, &x.hi, self.hi.as_float_mut());
323                } else if bfsub_is_one(&a, &b) && bfodd(&a) {
324                    self.set_prec(x.prec());
325                    self.lo.immovable = endpoint_min(&cosu, &x.lo, &x.hi, self.lo.as_float_mut());
326                    self.hi.as_float_mut().assign(1);
327                    self.hi.immovable = false;
328                } else {
329                    self.lo.as_float_mut().assign(-1);
330                    self.hi.as_float_mut().assign(1);
331                    self.lo.immovable = false;
332                    self.hi.immovable = false;
333                }
334            }
335        }
336    }
337
338    pub(crate) fn sinu_assign(&mut self, x: &Ival, n: u64) {
339        self.err = x.err;
340        let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
341        let period_qtr = period_quarter_bitlen(n, 4);
342        let sinu = |x: &Float, out: &mut Float, rnd: Round| mpfr_sinu(x, n, out, rnd);
343
344        match classify_ival_periodic(x, period_qtr) {
345            PeriodClass::TooWide => {
346                self.lo.as_float_mut().assign(-1);
347                self.hi.as_float_mut().assign(1);
348                self.lo.immovable = false;
349                self.hi.immovable = false;
350            }
351            PeriodClass::NearZero => {
352                self.monotonic_assign(&sinu, x);
353            }
354            PeriodClass::RangeReduce => {
355                let prec = range_reduce_precision(xlo, xhi, x.prec());
356                let (a, b) = ival_div_n_half(x, n, prec, mpfr_round_inplace);
357
358                if a == b && bfeven(&a) {
359                    self.monotonic_assign(&sinu, x);
360                } else if a == b && bfodd(&a) {
361                    self.comonotonic_assign(&sinu, x);
362                } else if bfsub_is_one(&a, &b) && bfodd(&a) {
363                    self.set_prec(x.prec());
364                    self.lo.as_float_mut().assign(-1);
365                    self.lo.immovable = false;
366                    self.hi.immovable = endpoint_max(&sinu, &x.lo, &x.hi, self.hi.as_float_mut());
367                } else if bfsub_is_one(&a, &b) && bfeven(&a) {
368                    self.set_prec(x.prec());
369                    self.lo.immovable = endpoint_min(&sinu, &x.lo, &x.hi, self.lo.as_float_mut());
370                    self.hi.as_float_mut().assign(1);
371                    self.hi.immovable = false;
372                } else {
373                    self.lo.as_float_mut().assign(-1);
374                    self.hi.as_float_mut().assign(1);
375                    self.lo.immovable = false;
376                    self.hi.immovable = false;
377                }
378            }
379        }
380    }
381
382    pub(crate) fn tanu_assign(&mut self, x: &Ival, n: u64) {
383        let (xlo, xhi) = (x.lo.as_float(), x.hi.as_float());
384        let period_qtr = period_quarter_bitlen(n, 8);
385        let immovable = x.lo.immovable && x.hi.immovable;
386        let tanu = |x: &Float, out: &mut Float, rnd: Round| mpfr_tanu(x, n, out, rnd);
387
388        match classify_ival_periodic(x, period_qtr) {
389            PeriodClass::TooWide => {
390                self.lo.as_float_mut().assign(f64::NEG_INFINITY);
391                self.hi.as_float_mut().assign(f64::INFINITY);
392                self.lo.immovable = immovable;
393                self.hi.immovable = immovable;
394                self.err.partial = true;
395                self.err.total = x.err.total;
396            }
397            PeriodClass::NearZero => {
398                self.monotonic_assign(&tanu, x);
399                self.err = x.err;
400            }
401            PeriodClass::RangeReduce => {
402                let prec = range_reduce_precision(xlo, xhi, x.prec());
403                let (a, b) = ival_div_n_half(x, n, prec, mpfr_round_inplace);
404
405                if a == b {
406                    self.monotonic_assign(&tanu, x);
407                    self.err = x.err;
408                } else {
409                    self.lo.as_float_mut().assign(f64::NEG_INFINITY);
410                    self.hi.as_float_mut().assign(f64::INFINITY);
411                    self.lo.immovable = immovable;
412                    self.hi.immovable = immovable;
413                    self.err.partial = true;
414                    self.err.total = x.err.total;
415                }
416            }
417        }
418    }
419}
rival/interval/trig.rs

rival/interval/
trig.rs
