use crate::InnerFloat::{Infinity, NaN, Zero};
use crate::TWICE_WIDTH;
use crate::float::arithmetic::exp::{exp_overflow, one_neighbor};
use crate::float::arithmetic::round_near_x::float_round_near_x;
use crate::{
Float, emulate_float_to_float_fn, emulate_rational_to_float_fn, float_infinity, float_nan,
float_zero, floor_and_ceiling,
};
use core::cmp::Ordering::{self, *};
use core::cmp::max;
use malachite_base::num::arithmetic::traits::{
CeilingLogBase2, ExpXMinus1, ExpXMinus1Assign, PowerOf2, Sign,
};
use malachite_base::num::basic::floats::PrimitiveFloat;
use malachite_base::num::basic::integers::PrimitiveInt;
use malachite_base::num::basic::traits::{NegativeOne, One, Zero as ZeroTrait};
use malachite_base::num::conversion::traits::{ExactFrom, RoundingFrom, SaturatingFrom};
use malachite_base::num::logic::traits::SignificantBits;
use malachite_base::rounding_modes::RoundingMode::{self, *};
use malachite_nz::integer::Integer;
use malachite_nz::natural::arithmetic::float_extras::float_can_round;
use malachite_nz::platform::{Limb, SignedLimb};
use malachite_q::Rational;
fn exp_x_minus_1_prec_round_normal(x: &Float, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
let ex = i64::from(x.get_exponent().unwrap());
if ex < 0 {
let (err, dir) = if *x > 0u32 {
(-ex, true)
} else {
(-ex + 1, false)
};
let err = u64::exact_from(err);
if err > prec + 1
&& let Some(result) = float_round_near_x(x, err, dir, prec, rm)
{
return result;
}
}
if x.is_sign_negative() && ex == i64::from(Float::MIN_EXPONENT) {
return exp_x_minus_1_rational_near_zero(&Rational::exact_from(x), prec, rm);
}
assert_ne!(rm, Exact, "Inexact exp_x_minus_1");
const BP: u64 = 64;
if x.is_sign_negative() && ex > 5 {
let log2_up = Float::ln_2_prec_round(BP, Up).0;
let t = x.div_prec_round_ref_val(log2_up, BP, Ceiling).0; let exp_t = i64::from(t.get_exponent().unwrap());
let (clamped, err, s_est) = if exp_t > 31 {
(
true,
u64::exact_from(Float::MAX_EXPONENT),
if exp_t > 64 {
u64::MAX
} else {
u64::power_of_2(u64::exact_from(exp_t - 1))
},
)
} else {
let neg_ceil = -Integer::rounding_from(&t, Ceiling).0;
const MAX_EXP: Integer = Integer::const_from_signed(Float::MAX_EXPONENT as SignedLimb);
let clamped = neg_ceil >= MAX_EXP;
let err = if clamped {
u64::exact_from(Float::MAX_EXPONENT)
} else {
u64::exact_from(&neg_ceil)
};
(clamped, err, u64::saturating_from(&neg_ceil))
};
if let Some(result) = float_round_near_x(&Float::NEGATIVE_ONE, err, false, prec, rm) {
return result;
}
if clamped {
return exp_x_minus_1_deep_negative(x, prec, rm, s_est);
}
}
let mut working_prec = prec + prec.ceiling_log_base_2() + 6;
if ex < 0 {
working_prec += u64::exact_from(-ex);
}
let mut increment = Limb::WIDTH;
loop {
let mut t = x.exp_prec_ref(working_prec).0;
if t.is_infinite() {
return exp_overflow(prec, rm);
}
let exp_te = i64::from(t.get_exponent().unwrap());
t.sub_prec_assign(Float::ONE, working_prec); let t_exp = i64::from(t.get_exponent().unwrap());
let err = working_prec - u64::exact_from(max(exp_te - t_exp, 0) + 1);
if float_can_round(t.significand_ref().unwrap(), err, prec, rm) {
return Float::from_float_prec_round(t, prec, rm);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
fn exp_x_minus_1_deep_negative(
x: &Float,
prec: u64,
rm: RoundingMode,
s_est: u64,
) -> (Float, Ordering) {
if s_est >= prec.saturating_add(2) {
return match rm {
Ceiling | Down => (-one_neighbor(prec, false), Greater), _ => (-Float::one_prec(prec), Less), };
}
let xr = Rational::exact_from(x);
let mut working_prec = prec.saturating_add(2).saturating_sub(s_est) + TWICE_WIDTH;
let mut increment = Limb::WIDTH;
loop {
let (ln_2_lo, ln_2_hi) = floor_and_ceiling(Float::ln_2_prec_round(working_prec, Floor));
let y_lo = &xr / Rational::exact_from(&ln_2_lo);
let y_hi = &xr / Rational::exact_from(&ln_2_hi);
let y_lo = Float::from_rational_prec_round(y_lo, working_prec, Floor).0;
let y_hi = Float::from_rational_prec_round(y_hi, working_prec, Ceiling).0;
let (e_lo, mut o_lo) = y_lo.power_of_2_x_minus_1_prec_round(prec, rm);
let (e_hi, mut o_hi) = y_hi.power_of_2_x_minus_1_prec_round(prec, rm);
if o_lo == Equal {
o_lo = o_hi;
}
if o_hi == Equal {
o_hi = o_lo;
}
if o_lo == o_hi && e_lo == e_hi {
return (e_lo, o_lo);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
pub(crate) fn exp_x_minus_1_rational_near_zero(
x: &Rational,
prec: u64,
rm: RoundingMode,
) -> (Float, Ordering) {
let negative = *x < 0u32;
let mut s = Rational::ZERO; let mut term = Rational::ONE; let mut k = 1u64;
loop {
term *= x;
term /= Rational::from(k); let s_next = &s + &term; let (lo, hi) = if negative {
if s < s_next {
(s.clone(), s_next.clone())
} else {
(s_next.clone(), s.clone())
}
} else {
let next = (&term * x) / Rational::from(k + 1); (s_next.clone(), &s_next + next / (Rational::ONE - x))
};
s = s_next;
k += 1;
let (f_lo, mut o_lo) = Float::from_rational_prec_round_ref(&lo, prec, rm);
let (f_hi, mut o_hi) = Float::from_rational_prec_round_ref(&hi, prec, rm);
if o_lo == Equal {
o_lo = o_hi;
}
if o_hi == Equal {
o_hi = o_lo;
}
if o_lo == o_hi && f_lo == f_hi {
return (f_lo, o_lo);
}
}
}
fn exp_x_minus_1_rational_helper(x: &Rational, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
assert_ne!(rm, Exact, "Inexact exp_x_minus_1");
let positive = x.sign() == Greater;
let exp_x = x.floor_log_base_2_abs() + 1; if exp_x <= Float::MIN_EXPONENT_I64 {
return exp_x_minus_1_rational_near_zero(x, prec, rm);
}
if exp_x >= Float::MAX_EXPONENT_I64 {
if positive {
return exp_overflow(prec, rm);
}
let err = u64::exact_from(Float::MAX_EXPONENT);
if let Some(result) = float_round_near_x(&Float::NEGATIVE_ONE, err, false, prec, rm) {
return result;
}
return match rm {
Ceiling | Down => (-one_neighbor(prec, false), Greater), _ => (-Float::one_prec(prec), Less), };
}
let mut working_prec = prec + 10;
let mut increment = Limb::WIDTH;
loop {
let (x_lo, x_o) = Float::from_rational_prec_round_ref(x, working_prec, Floor);
if x_o == Equal {
return x_lo.exp_x_minus_1_prec_round(prec, rm);
}
let (x_lo, x_hi) = floor_and_ceiling((x_lo, x_o));
let (e_lo, o_lo) = x_lo.exp_x_minus_1_prec_round_ref(prec, rm);
let (e_hi, o_hi) = x_hi.exp_x_minus_1_prec_round_ref(prec, rm);
if o_lo == o_hi && e_lo == e_hi {
return (e_lo, o_lo);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
impl Float {
#[inline]
pub fn exp_x_minus_1_prec_round(self, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
self.exp_x_minus_1_prec_round_ref(prec, rm)
}
#[inline]
pub fn exp_x_minus_1_prec_round_ref(&self, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
assert_ne!(prec, 0);
match self {
Self(NaN) => (float_nan!(), Equal),
float_infinity!() => (float_infinity!(), Equal),
Self(Infinity { sign: false }) => (Self::from_signed_prec(-1i32, prec).0, Equal),
Self(Zero { sign }) => (Self(Zero { sign: *sign }), Equal),
_ => exp_x_minus_1_prec_round_normal(self, prec, rm),
}
}
#[inline]
pub fn exp_x_minus_1_prec(self, prec: u64) -> (Self, Ordering) {
self.exp_x_minus_1_prec_round(prec, Nearest)
}
#[inline]
pub fn exp_x_minus_1_prec_ref(&self, prec: u64) -> (Self, Ordering) {
self.exp_x_minus_1_prec_round_ref(prec, Nearest)
}
#[allow(clippy::needless_pass_by_value)]
#[inline]
pub fn exp_x_minus_1_rational_prec_round(
x: Rational,
prec: u64,
rm: RoundingMode,
) -> (Self, Ordering) {
Self::exp_x_minus_1_rational_prec_round_ref(&x, prec, rm)
}
pub fn exp_x_minus_1_rational_prec_round_ref(
x: &Rational,
prec: u64,
rm: RoundingMode,
) -> (Self, Ordering) {
assert_ne!(prec, 0);
if *x == 0u32 {
return (float_zero!(), Equal);
}
exp_x_minus_1_rational_helper(x, prec, rm)
}
#[allow(clippy::needless_pass_by_value)]
#[inline]
pub fn exp_x_minus_1_rational_prec(x: Rational, prec: u64) -> (Self, Ordering) {
Self::exp_x_minus_1_rational_prec_round_ref(&x, prec, Nearest)
}
#[inline]
pub fn exp_x_minus_1_rational_prec_ref(x: &Rational, prec: u64) -> (Self, Ordering) {
Self::exp_x_minus_1_rational_prec_round_ref(x, prec, Nearest)
}
#[inline]
pub fn exp_x_minus_1_round(self, rm: RoundingMode) -> (Self, Ordering) {
let prec = self.significant_bits();
self.exp_x_minus_1_prec_round(prec, rm)
}
#[inline]
pub fn exp_x_minus_1_round_ref(&self, rm: RoundingMode) -> (Self, Ordering) {
self.exp_x_minus_1_prec_round_ref(self.significant_bits(), rm)
}
#[inline]
pub fn exp_x_minus_1_prec_round_assign(&mut self, prec: u64, rm: RoundingMode) -> Ordering {
let (result, o) = core::mem::take(self).exp_x_minus_1_prec_round(prec, rm);
*self = result;
o
}
#[inline]
pub fn exp_x_minus_1_prec_assign(&mut self, prec: u64) -> Ordering {
self.exp_x_minus_1_prec_round_assign(prec, Nearest)
}
#[inline]
pub fn exp_x_minus_1_round_assign(&mut self, rm: RoundingMode) -> Ordering {
let prec = self.significant_bits();
self.exp_x_minus_1_prec_round_assign(prec, rm)
}
}
impl ExpXMinus1 for Float {
type Output = Self;
#[inline]
fn exp_x_minus_1(self) -> Self {
let prec = self.significant_bits();
self.exp_x_minus_1_prec(prec).0
}
}
impl ExpXMinus1 for &Float {
type Output = Float;
#[inline]
fn exp_x_minus_1(self) -> Float {
self.exp_x_minus_1_prec_round_ref(self.significant_bits(), Nearest)
.0
}
}
impl ExpXMinus1Assign for Float {
#[inline]
fn exp_x_minus_1_assign(&mut self) {
let prec = self.significant_bits();
self.exp_x_minus_1_prec_round_assign(prec, Nearest);
}
}
#[inline]
#[allow(clippy::type_repetition_in_bounds)]
pub fn primitive_float_exp_x_minus_1<T: PrimitiveFloat>(x: T) -> T
where
Float: From<T> + PartialOrd<T>,
for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
{
emulate_float_to_float_fn(Float::exp_x_minus_1_prec, x)
}
#[inline]
#[allow(clippy::type_repetition_in_bounds)]
pub fn primitive_float_exp_x_minus_1_rational<T: PrimitiveFloat>(x: &Rational) -> T
where
Float: PartialOrd<T>,
for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
{
emulate_rational_to_float_fn(Float::exp_x_minus_1_rational_prec_ref, x)
}