use crate::{classify::*, const_interval, interval::*};
use gmp_mpfr_sys::mpfr;
use rug::Float;
fn mpfr_fn(
f: unsafe extern "C" fn(*mut mpfr::mpfr_t, *const mpfr::mpfr_t, mpfr::rnd_t) -> i32,
x: f64,
rnd: mpfr::rnd_t,
) -> f64 {
let mut x = Float::with_val(f64::MANTISSA_DIGITS, x);
unsafe {
f(x.as_raw_mut(), x.as_raw(), rnd);
mpfr::get_d(x.as_raw(), rnd)
}
}
fn mpfr_fn2(
f: unsafe extern "C" fn(
*mut mpfr::mpfr_t,
*const mpfr::mpfr_t,
*const mpfr::mpfr_t,
mpfr::rnd_t,
) -> i32,
x: f64,
y: f64,
rnd: mpfr::rnd_t,
) -> f64 {
let mut x = Float::with_val(f64::MANTISSA_DIGITS, x);
let y = Float::with_val(f64::MANTISSA_DIGITS, y);
unsafe {
f(x.as_raw_mut(), x.as_raw(), y.as_raw(), rnd);
mpfr::get_d(x.as_raw(), rnd)
}
}
fn mpfr_fn_si(
f: unsafe extern "C" fn(*mut mpfr::mpfr_t, *const mpfr::mpfr_t, i64, mpfr::rnd_t) -> i32,
x: f64,
y: i64,
rnd: mpfr::rnd_t,
) -> f64 {
let mut x = Float::with_val(f64::MANTISSA_DIGITS, x);
unsafe {
f(x.as_raw_mut(), x.as_raw(), y, rnd);
mpfr::get_d(x.as_raw(), rnd)
}
}
macro_rules! mpfr_fn {
($mpfr_f:ident, $f_rd:ident, $f_ru:ident) => {
fn $f_rd(x: f64) -> f64 {
mpfr_fn(mpfr::$mpfr_f, x, mpfr::rnd_t::RNDD)
}
fn $f_ru(x: f64) -> f64 {
mpfr_fn(mpfr::$mpfr_f, x, mpfr::rnd_t::RNDU)
}
};
}
macro_rules! mpfr_fn2 {
($mpfr_f:ident, $f_rd:ident, $f_ru:ident) => {
fn $f_rd(x: f64, y: f64) -> f64 {
mpfr_fn2(mpfr::$mpfr_f, x, y, mpfr::rnd_t::RNDD)
}
fn $f_ru(x: f64, y: f64) -> f64 {
mpfr_fn2(mpfr::$mpfr_f, x, y, mpfr::rnd_t::RNDU)
}
};
}
macro_rules! mpfr_fn_si {
($mpfr_f:ident, $f_rd:ident, $f_ru:ident) => {
fn $f_rd(x: f64, y: i32) -> f64 {
mpfr_fn_si(mpfr::$mpfr_f, x, y as i64, mpfr::rnd_t::RNDD)
}
fn $f_ru(x: f64, y: i32) -> f64 {
mpfr_fn_si(mpfr::$mpfr_f, x, y as i64, mpfr::rnd_t::RNDU)
}
};
}
mpfr_fn!(acos, acos_rd, acos_ru);
mpfr_fn!(acosh, acosh_rd, acosh_ru);
mpfr_fn!(asin, asin_rd, asin_ru);
mpfr_fn!(asinh, asinh_rd, asinh_ru);
mpfr_fn!(atan, atan_rd, atan_ru);
mpfr_fn2!(atan2, atan2_rd, atan2_ru);
mpfr_fn!(atanh, atanh_rd, atanh_ru);
mpfr_fn!(cos, cos_rd, cos_ru);
mpfr_fn!(cosh, cosh_rd, cosh_ru);
mpfr_fn!(exp, exp_rd, exp_ru);
mpfr_fn!(exp10, exp10_rd, exp10_ru);
mpfr_fn!(exp2, exp2_rd, exp2_ru);
mpfr_fn!(log, ln_rd, ln_ru);
mpfr_fn!(log10, log10_rd, log10_ru);
mpfr_fn!(log2, log2_rd, log2_ru);
mpfr_fn2!(pow, pow_rd, pow_ru);
mpfr_fn_si!(pow_si, pown_rd, pown_ru);
mpfr_fn!(sin, sin_rd, sin_ru);
mpfr_fn!(sinh, sinh_rd, sinh_ru);
mpfr_fn!(tan, tan_rd, tan_ru);
mpfr_fn!(tanh, tanh_rd, tanh_ru);
fn rem_euclid_2(x: f64) -> f64 {
if 2.0 * (x / 2.0).floor() == x {
0.0
} else {
1.0
}
}
macro_rules! impl_log {
($(#[$meta:meta])* $f:ident, $f_impl:ident, $f_rd:ident, $f_ru:ident) => {
$(#[$meta])*
pub fn $f(self) -> Self {
self.$f_impl().0
}
#[allow(clippy::many_single_char_names)]
fn $f_impl(self) -> (Self, Decoration) {
const DOM: Interval = const_interval!(0.0, f64::INFINITY);
let x = self.intersection(DOM);
let a = x.inf_raw();
let b = x.sup_raw();
if x.is_empty() || b <= 0.0 {
return (Self::EMPTY, Decoration::Trv);
}
let y = Self::with_infsup_raw($f_rd(a), $f_ru(b));
let d = if self.interior(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
(y, d)
}
};
}
macro_rules! impl_mono_inc {
($(#[$meta:meta])* $f:ident, $f_rd:ident, $f_ru:ident) => {
$(#[$meta])*
pub fn $f(self) -> Self {
if self.is_empty() {
return self;
}
Self::with_infsup_raw($f_rd(self.inf_raw()), $f_ru(self.sup_raw()))
}
};
}
impl Interval {
pub fn acos(self) -> Self {
self.acos_impl().0
}
fn acos_impl(self) -> (Self, Decoration) {
const DOM: Interval = const_interval!(-1.0, 1.0);
let x = self.intersection(DOM);
if x.is_empty() {
return (x, Decoration::Trv);
}
let y = Self::with_infsup_raw(acos_rd(x.sup_raw()), acos_ru(x.inf_raw()));
let d = if self.subset(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
(y, d)
}
pub fn acosh(self) -> Self {
self.acosh_impl().0
}
fn acosh_impl(self) -> (Self, Decoration) {
const DOM: Interval = const_interval!(1.0, f64::INFINITY);
let x = self.intersection(DOM);
if x.is_empty() {
return (x, Decoration::Trv);
}
let y = Self::with_infsup_raw(acosh_rd(x.inf_raw()), acosh_ru(x.sup_raw()));
let d = if self.subset(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
(y, d)
}
pub fn asin(self) -> Self {
self.asin_impl().0
}
fn asin_impl(self) -> (Self, Decoration) {
const DOM: Interval = const_interval!(-1.0, 1.0);
let x = self.intersection(DOM);
if x.is_empty() {
return (x, Decoration::Trv);
}
let y = Self::with_infsup_raw(asin_rd(x.inf_raw()), asin_ru(x.sup_raw()));
let d = if self.subset(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
(y, d)
}
impl_mono_inc!(
asinh,
asinh_rd,
asinh_ru
);
impl_mono_inc!(
atan,
atan_rd,
atan_ru
);
pub fn atan2(self, rhs: Self) -> Self {
self.atan2_impl(rhs).0
}
#[allow(clippy::many_single_char_names)]
fn atan2_impl(self, rhs: Self) -> (Self, Decoration) {
let (x, y) = (rhs, self);
let a = x.inf_raw();
let b = x.sup_raw();
let c = y.inf_raw();
let d = y.sup_raw();
use IntervalClass2::*;
match x.classify2(y) {
E_E | E_M | E_N0 | E_N1 | E_P0 | E_P1 | E_Z | M_E | N0_E | N1_E | P0_E | P1_E | Z_E
| Z_Z => (Self::EMPTY, Decoration::Trv),
M_M | M_N0 | N0_M | N0_N0 => (
Self::with_infsup_raw(-Self::PI.sup_raw(), Self::PI.sup_raw()),
Decoration::Trv,
),
P0_P0 => (
Self::with_infsup_raw(0.0, Self::FRAC_PI_2.sup_raw()),
Decoration::Trv,
),
P0_P1 | P1_P0 | P1_P1 | P1_Z | Z_P1 => (
Self::with_infsup_raw(atan2_rd(c, b), atan2_ru(d, a)),
Decoration::Com,
),
M_P0 | M_Z => (
Self::with_infsup_raw(0.0, Self::PI.sup_raw()),
Decoration::Trv,
),
M_P1 => (
Self::with_infsup_raw(atan2_rd(c, b), atan2_ru(c, a)),
Decoration::Com,
),
N0_P0 => (
Self::with_infsup_raw(Self::FRAC_PI_2.inf_raw(), Self::PI.sup_raw()),
Decoration::Trv,
),
N0_P1 | N1_P1 => (
Self::with_infsup_raw(atan2_rd(d, b), atan2_ru(c, a)),
Decoration::Com,
),
N1_P0 => (
Self::with_infsup_raw(atan2_rd(d, b), Self::PI.sup_raw()),
Decoration::Dac,
),
N1_M | N1_N0 => (
Self::with_infsup_raw(-Self::PI.sup_raw(), Self::PI.sup_raw()),
Decoration::Def,
),
N0_N1 | N1_N1 => (
Self::with_infsup_raw(atan2_rd(d, a), atan2_ru(c, b)),
Decoration::Com,
),
M_N1 => (
Self::with_infsup_raw(atan2_rd(d, a), atan2_ru(d, b)),
Decoration::Com,
),
P0_N0 => (
Self::with_infsup_raw(-Self::FRAC_PI_2.sup_raw(), 0.0),
Decoration::Trv,
),
P0_N1 | P1_N0 | P1_N1 | Z_N1 => (
Self::with_infsup_raw(atan2_rd(c, a), atan2_ru(d, b)),
Decoration::Com,
),
P0_M | Z_M => (
Self::with_infsup_raw(-Self::FRAC_PI_2.sup_raw(), Self::FRAC_PI_2.sup_raw()),
Decoration::Trv,
),
P1_M => (
Self::with_infsup_raw(atan2_rd(c, a), atan2_ru(d, a)),
Decoration::Com,
),
N0_Z => (Self::PI, Decoration::Trv),
N1_Z => (Self::PI, Decoration::Dac),
P0_Z => (Self::zero(), Decoration::Trv),
Z_N0 => (-Self::FRAC_PI_2, Decoration::Trv),
Z_P0 => (Self::FRAC_PI_2, Decoration::Trv),
}
}
pub fn atanh(self) -> Self {
self.atanh_impl().0
}
#[allow(clippy::many_single_char_names)]
fn atanh_impl(self) -> (Self, Decoration) {
const DOM: Interval = const_interval!(-1.0, 1.0);
let x = self.intersection(DOM);
let a = x.inf_raw();
let b = x.sup_raw();
if x.is_empty() || b <= -1.0 || a >= 1.0 {
return (Self::EMPTY, Decoration::Trv);
}
let y = Self::with_infsup_raw(atanh_rd(a), atanh_ru(b));
let d = if self.interior(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
(y, d)
}
pub fn cos(self) -> Self {
if self.is_empty() {
return self;
}
let a = self.inf_raw();
let b = self.sup_raw();
let q_nowrap = (self / Self::PI).floor();
let qa = q_nowrap.inf_raw();
let qb = q_nowrap.sup_raw();
let n = if a == b {
0.0
} else {
qb - qa
};
let q = rem_euclid_2(qa);
if n == 0.0 {
if q == 0.0 {
Self::with_infsup_raw(cos_rd(b), cos_ru(a))
} else {
Self::with_infsup_raw(cos_rd(a), cos_ru(b))
}
} else if n <= 1.0 {
if q == 0.0 {
Self::with_infsup_raw(-1.0, cos_ru(a).max(cos_ru(b)))
} else {
Self::with_infsup_raw(cos_rd(a).min(cos_rd(b)), 1.0)
}
} else {
const_interval!(-1.0, 1.0)
}
}
pub fn cosh(self) -> Self {
if self.is_empty() {
return self;
}
let a = self.inf_raw();
let b = self.sup_raw();
if b < 0.0 {
Self::with_infsup_raw(cosh_rd(b), cosh_ru(a))
} else if a > 0.0 {
Self::with_infsup_raw(cosh_rd(a), cosh_ru(b))
} else {
Self::with_infsup_raw(1.0, cosh_ru((-a).max(b)))
}
}
impl_mono_inc!(
exp,
exp_rd,
exp_ru
);
impl_mono_inc!(
exp10,
exp10_rd,
exp10_ru
);
impl_mono_inc!(
exp2,
exp2_rd,
exp2_ru
);
impl_log!(
ln,
ln_impl,
ln_rd,
ln_ru
);
impl_log!(
log10,
log10_impl,
log10_rd,
log10_ru
);
impl_log!(
log2,
log2_impl,
log2_rd,
log2_ru
);
pub fn pow(self, rhs: Self) -> Self {
self.pow_impl(rhs).0
}
#[allow(clippy::many_single_char_names)]
fn pow_impl(self, rhs: Self) -> (Self, Decoration) {
const DOM: Interval = const_interval!(0.0, f64::INFINITY);
let x = self.intersection(DOM);
if x.either_empty(rhs) {
return (Self::EMPTY, Decoration::Trv);
}
let a = x.inf_raw();
let b = x.sup_raw();
let c = rhs.inf_raw();
let d = rhs.sup_raw();
if d <= 0.0 {
if b == 0.0 {
return (Self::EMPTY, Decoration::Trv);
}
let dec = if self.interior(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
if b < 1.0 {
(Self::with_infsup_raw(pow_rd(b, d), pow_ru(a, c)), dec)
} else if a > 1.0 {
(Self::with_infsup_raw(pow_rd(b, c), pow_ru(a, d)), dec)
} else {
(Self::with_infsup_raw(pow_rd(b, c), pow_ru(a, c)), dec)
}
} else if c > 0.0 {
let dec = if self.subset(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
if b < 1.0 {
(Self::with_infsup_raw(pow_rd(a, d), pow_ru(b, c)), dec)
} else if a > 1.0 {
(Self::with_infsup_raw(pow_rd(a, c), pow_ru(b, d)), dec)
} else {
(Self::with_infsup_raw(pow_rd(a, d), pow_ru(b, d)), dec)
}
} else {
if b == 0.0 {
return (Self::zero(), Decoration::Trv);
}
let z_ac = pow_ru(a, c);
let z_ad = pow_rd(a, d);
let z_bc = pow_rd(b, c);
let z_bd = pow_ru(b, d);
let dec = if self.interior(DOM) {
Decoration::Com
} else {
Decoration::Trv
};
(Self::with_infsup_raw(z_ad.min(z_bc), z_ac.max(z_bd)), dec)
}
}
pub fn pown(self, rhs: i32) -> Self {
self.pown_impl(rhs).0
}
fn pown_impl(self, rhs: i32) -> (Self, Decoration) {
if self.is_empty() {
return (self, Decoration::Trv);
}
let mut a = self.inf_raw();
let mut b = self.sup_raw();
#[allow(clippy::collapsible_else_if, clippy::collapsible_if)]
if rhs < 0 {
let d = if a <= 0.0 && b >= 0.0 {
Decoration::Trv
} else {
Decoration::Com
};
if a == 0.0 && b == 0.0 {
return (Self::EMPTY, d);
}
if rhs % 2 == 0 {
let abs = self.abs();
(
Self::with_infsup_raw(pown_rd(abs.sup_raw(), rhs), pown_ru(abs.inf_raw(), rhs)),
d,
)
} else {
if a < 0.0 && b > 0.0 {
(Self::ENTIRE, d)
} else {
if a == 0.0 {
a = 0.0; }
if b == 0.0 {
b = -0.0; }
(Self::with_infsup_raw(pown_rd(b, rhs), pown_ru(a, rhs)), d)
}
}
} else {
if rhs % 2 == 0 {
let abs = self.abs();
(
Self::with_infsup_raw(pown_rd(abs.inf_raw(), rhs), pown_ru(abs.sup_raw(), rhs)),
Decoration::Com,
)
} else {
(
Self::with_infsup_raw(pown_rd(a, rhs), pown_ru(b, rhs)),
Decoration::Com,
)
}
}
}
pub fn sin(self) -> Self {
if self.is_empty() {
return self;
}
let a = self.inf_raw();
let b = self.sup_raw();
let q_nowrap = (self / Self::FRAC_PI_2).floor();
let qa = q_nowrap.inf_raw();
let qb = q_nowrap.sup_raw();
let n = if a == b { 0.0 } else { qb - qa };
let q = qa.rem_euclid(4.0);
if q == 0.0 && n < 1.0 || q == 3.0 && n < 2.0 {
Self::with_infsup_raw(sin_rd(a), sin_ru(b))
} else if q == 1.0 && n < 2.0 || q == 2.0 && n < 1.0 {
Self::with_infsup_raw(sin_rd(b), sin_ru(a))
} else if q == 0.0 && n < 3.0 || q == 3.0 && n < 4.0 {
Self::with_infsup_raw(sin_rd(a).min(sin_rd(b)), 1.0)
} else if q == 1.0 && n < 4.0 || q == 2.0 && n < 3.0 {
Self::with_infsup_raw(-1.0, sin_ru(a).max(sin_ru(b)))
} else {
const_interval!(-1.0, 1.0)
}
}
impl_mono_inc!(
sinh,
sinh_rd,
sinh_ru
);
pub fn tan(self) -> Self {
self.tan_impl().0
}
fn tan_impl(self) -> (Self, Decoration) {
if self.is_empty() {
return (self, Decoration::Trv);
}
let a = self.inf_raw();
let b = self.sup_raw();
let q_nowrap = (self / Self::FRAC_PI_2).floor();
let qa = q_nowrap.inf_raw();
let qb = q_nowrap.sup_raw();
let n = if a == b { 0.0 } else { qb - qa };
let q = rem_euclid_2(qa);
let cont =
qb != f64::INFINITY && b <= (Self::with_infsup_raw(qb, qb) * Self::FRAC_PI_2).inf_raw();
if q == 0.0 && (n < 1.0 || n == 1.0 && cont) || q == 1.0 && (n < 2.0 || n == 2.0 && cont) {
(Self::with_infsup_raw(tan_rd(a), tan_ru(b)), Decoration::Com)
} else {
(Self::ENTIRE, Decoration::Trv)
}
}
impl_mono_inc!(
tanh,
tanh_rd,
tanh_ru
);
}
macro_rules! impl_dec {
($f:ident) => {
pub fn $f(self) -> Self {
Self::set_dec(self.x.$f(), self.d)
}
};
($f:ident, $f_impl:ident) => {
pub fn $f(self) -> Self {
let (y, d) = self.x.$f_impl();
Self::set_dec(y, self.d.min(d))
}
};
}
macro_rules! impl_dec2 {
($f:ident, $f_impl:ident) => {
pub fn $f(self, rhs: Self) -> Self {
let (z, d) = self.x.$f_impl(rhs.x);
Self::set_dec(z, self.d.min(rhs.d.min(d)))
}
};
}
impl DecInterval {
impl_dec!(acos, acos_impl);
impl_dec!(acosh, acosh_impl);
impl_dec!(asin, asin_impl);
impl_dec!(asinh);
impl_dec!(atan);
impl_dec2!(atan2, atan2_impl);
impl_dec!(atanh, atanh_impl);
impl_dec!(cos);
impl_dec!(cosh);
impl_dec!(exp);
impl_dec!(exp10);
impl_dec!(exp2);
impl_dec!(ln, ln_impl);
impl_dec!(log10, log10_impl);
impl_dec!(log2, log2_impl);
impl_dec2!(pow, pow_impl);
pub fn pown(self, rhs: i32) -> Self {
let (y, d) = self.x.pown_impl(rhs);
Self::set_dec(y, self.d.min(d))
}
impl_dec!(sin);
impl_dec!(sinh);
impl_dec!(tan, tan_impl);
impl_dec!(tanh);
}
#[cfg(test)]
mod tests {
use crate::*;
use DecInterval as DI;
use Interval as I;
#[test]
fn nai() {
assert!(DI::NAI.acos().is_nai());
assert!(DI::NAI.acosh().is_nai());
assert!(DI::NAI.asin().is_nai());
assert!(DI::NAI.asinh().is_nai());
assert!(DI::NAI.atan().is_nai());
assert!(DI::NAI.atan2(DI::EMPTY).is_nai());
assert!(DI::EMPTY.atan2(DI::NAI).is_nai());
assert!(DI::NAI.atanh().is_nai());
assert!(DI::NAI.cos().is_nai());
assert!(DI::NAI.cosh().is_nai());
assert!(DI::NAI.exp().is_nai());
assert!(DI::NAI.exp10().is_nai());
assert!(DI::NAI.exp2().is_nai());
assert!(DI::NAI.ln().is_nai());
assert!(DI::NAI.log10().is_nai());
assert!(DI::NAI.log2().is_nai());
assert!(DI::NAI.pow(DI::EMPTY).is_nai());
assert!(DI::EMPTY.pow(DI::NAI).is_nai());
assert!(DI::NAI.pown(1).is_nai());
assert!(DI::NAI.sin().is_nai());
assert!(DI::NAI.sinh().is_nai());
assert!(DI::NAI.tan().is_nai());
assert!(DI::NAI.tanh().is_nai());
}
#[test]
fn tan() {
assert!(interval!(std::f64::consts::FRAC_PI_4, I::FRAC_PI_2.inf())
.unwrap()
.tan()
.is_common_interval());
assert!(interval!(-std::f64::consts::FRAC_PI_4, I::FRAC_PI_2.inf())
.unwrap()
.tan()
.is_common_interval());
assert!(
interval!(-3.0 * std::f64::consts::FRAC_PI_4, -I::FRAC_PI_2.sup())
.unwrap()
.tan()
.is_common_interval()
);
assert!(
interval!(-5.0 * std::f64::consts::FRAC_PI_4, -I::FRAC_PI_2.sup())
.unwrap()
.tan()
.is_common_interval()
);
}
}