use crate::InnerFloat::{Finite, Infinity, NaN, Zero};
use crate::float::basic::extended::ExtendedFloat;
use crate::{
Float, emulate_float_to_float_fn, emulate_rational_to_float_fn, float_either_zero,
float_infinity, float_nan, float_negative_infinity, float_zero, floor_and_ceiling,
significand_bits,
};
use alloc::vec;
use core::cmp::Ordering::{self, *};
use core::mem::{swap, take};
use malachite_base::num::arithmetic::traits::{
Abs, Agm, CeilingLogBase2, IsPowerOf2, Ln, LnAssign, Parity, PowerOf2, Sign,
};
use malachite_base::num::basic::floats::PrimitiveFloat;
use malachite_base::num::basic::integers::PrimitiveInt;
use malachite_base::num::basic::traits::{NegativeInfinity, 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;
use malachite_q::Rational;
pub(crate) fn ln_1_plus_rational_brackets(e: &Rational, wprec: u64) -> (Rational, Rational) {
let negative = *e < 0u32;
let mut pow = e.clone(); let mut s = e.clone(); let mut k = 1u64;
loop {
pow *= e; k += 1;
let mut term = &pow / Rational::from(k); if k.even() {
term = -term;
}
let s_next = &s + &term; let (lo, hi) = if negative {
let bound = -((&pow * e) / (Rational::from(k + 1) * (Rational::ONE + e))).abs();
(&s_next + &bound, s_next.clone())
} else if s < s_next {
(s.clone(), s_next.clone())
} else {
(s_next.clone(), s.clone())
};
s = s_next;
let width = &hi - &lo;
if width == 0u32
|| width.floor_log_base_2_abs() < lo.floor_log_base_2_abs() - i64::exact_from(wprec) - 2
{
return (lo, hi);
}
}
}
pub(crate) enum SliverOfOne {
No,
Representable(Float),
Underflow,
}
pub(crate) fn sliver_of_one(x: &Float) -> SliverOfOne {
let e = i64::from(x.get_exponent().unwrap());
if (e == 0 || e == 1)
&& x.get_prec().unwrap() >= u64::exact_from(-i64::from(Float::MIN_EXPONENT) - 8)
{
let (mut d, o) = x.sub_prec_round_ref_val(Float::ONE, x.get_prec().unwrap() + 1, Floor);
if o != Equal {
return SliverOfOne::Underflow;
}
if i64::from(d.get_exponent().unwrap()) <= i64::from(Float::MIN_EXPONENT) + 4 {
let sig = d.significand_ref().unwrap();
let min_prec = significand_bits(sig) - sig.trailing_zeros().unwrap();
let o = d.set_prec_round(min_prec, Floor);
debug_assert_eq!(o, Equal);
return SliverOfOne::Representable(d);
}
}
SliverOfOne::No
}
fn ln_prec_round_normal_ref(x: &Float, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
if *x == 1u32 {
return (Float::ZERO, Equal);
}
match sliver_of_one(x) {
SliverOfOne::Representable(d) => return d.ln_1_plus_x_prec_round(prec, rm),
SliverOfOne::Underflow => {
return Float::ln_rational_prec_round(Rational::exact_from(x), prec, rm);
}
SliverOfOne::No => {}
}
assert_ne!(rm, Exact, "Inexact ln");
let x_exp = i64::from(x.get_exponent().unwrap());
let mut working_prec = prec + (prec.ceiling_log_base_2() << 1) + 10;
let mut increment = Limb::WIDTH;
let mut previous_m = 0;
let mut x = x.clone();
loop {
let m = i64::exact_from((working_prec + 3) >> 1)
.checked_sub(x_exp)
.unwrap();
x <<= m - previous_m;
previous_m = m;
assert!(x.is_normal());
let tmp2 = Float::pi_prec(working_prec).0
/ (Float::ONE.agm(
const { Float::const_from_unsigned(4) }
.div_prec_round_val_ref(&x, working_prec, Floor)
.0,
) << 1u32);
let exp2 = tmp2.get_exponent();
let tmp1 = tmp2
- Float::ln_2_prec(working_prec)
.0
.mul_prec(Float::from(m), working_prec)
.0;
if let (Some(exp1), Some(exp2)) = (tmp1.get_exponent(), exp2) {
let cancel = u64::saturating_from(exp2 - exp1);
if float_can_round(
tmp1.significand_ref().unwrap(),
working_prec.saturating_sub(cancel).saturating_sub(4),
prec,
rm,
) {
return Float::from_float_prec_round(tmp1, prec, rm);
}
working_prec += cancel + working_prec.ceiling_log_base_2();
} else {
working_prec += working_prec.ceiling_log_base_2();
}
working_prec += increment;
increment = working_prec >> 1;
}
}
fn ln_prec_round_normal(mut x: Float, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
if x == 1u32 {
return (Float::ZERO, Equal);
}
match sliver_of_one(&x) {
SliverOfOne::Representable(d) => return d.ln_1_plus_x_prec_round(prec, rm),
SliverOfOne::Underflow => {
return Float::ln_rational_prec_round(Rational::exact_from(&x), prec, rm);
}
SliverOfOne::No => {}
}
assert_ne!(rm, Exact, "Inexact ln");
let x_exp = i64::from(x.get_exponent().unwrap());
let mut working_prec = prec + (prec.ceiling_log_base_2() << 1) + 10;
let mut increment = Limb::WIDTH;
let mut previous_m = 0;
loop {
let m = i64::exact_from((working_prec + 3) >> 1)
.checked_sub(x_exp)
.unwrap();
x <<= m - previous_m;
previous_m = m;
assert!(x.is_normal());
let tmp2 = Float::pi_prec(working_prec).0
/ (Float::ONE.agm(
const { Float::const_from_unsigned(4) }
.div_prec_round_val_ref(&x, working_prec, Floor)
.0,
) << 1u32);
let exp2 = tmp2.get_exponent();
let tmp1 = tmp2
- Float::ln_2_prec(working_prec)
.0
.mul_prec(Float::from(m), working_prec)
.0;
if let (Some(exp1), Some(exp2)) = (tmp1.get_exponent(), exp2) {
let cancel = u64::saturating_from(exp2 - exp1);
if float_can_round(
tmp1.significand_ref().unwrap(),
working_prec.saturating_sub(cancel).saturating_sub(4),
prec,
rm,
) {
return Float::from_float_prec_round(tmp1, prec, rm);
}
working_prec += cancel + working_prec.ceiling_log_base_2();
} else {
working_prec += working_prec.ceiling_log_base_2();
}
working_prec += increment;
increment = working_prec >> 1;
}
}
pub(crate) fn ln_prec_round_normal_extended(
x: ExtendedFloat,
prec: u64,
rm: RoundingMode,
) -> (Float, Ordering) {
if x.exp == 1 && x.x.is_power_of_2() {
return (Float::ZERO, Equal);
}
assert_ne!(rm, Exact, "Inexact ln");
let x_exp = x.exp;
let mut working_prec = prec + (prec.ceiling_log_base_2() << 1) + 10;
let mut increment = Limb::WIDTH;
let mut m = i64::exact_from((working_prec + 3) >> 1)
.checked_sub(x.exp)
.unwrap();
let mut previous_m = m;
let mut x = Float::exact_from(x << m);
let mut first = true;
loop {
if first {
first = false;
} else {
m = i64::exact_from((working_prec + 3) >> 1)
.checked_sub(x_exp)
.unwrap();
x <<= m - previous_m;
previous_m = m;
}
assert!(x.is_normal());
let tmp2 = Float::pi_prec(working_prec).0
/ (Float::ONE.agm(
const { Float::const_from_unsigned(4) }
.div_prec_round_val_ref(&x, working_prec, Floor)
.0,
) << 1u32);
let exp2 = tmp2.get_exponent();
let tmp1 = tmp2
- Float::ln_2_prec(working_prec)
.0
.mul_prec(Float::from(m), working_prec)
.0;
if let (Some(exp1), Some(exp2)) = (tmp1.get_exponent(), exp2) {
let cancel = u64::saturating_from(exp2 - exp1);
if float_can_round(
tmp1.significand_ref().unwrap(),
working_prec.saturating_sub(cancel).saturating_sub(4),
prec,
rm,
) {
return Float::from_float_prec_round(tmp1, prec, rm);
}
working_prec += cancel + working_prec.ceiling_log_base_2();
} else {
working_prec += working_prec.ceiling_log_base_2();
}
working_prec += increment;
increment = working_prec >> 1;
}
}
fn ln_rational_near_one(eps: &Rational, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
let negative = *eps < 0u32;
let mut pow = eps.clone(); let mut s = eps.clone(); let mut k = 1u64;
loop {
pow *= eps;
k += 1;
let mut term = &pow / Rational::from(k); if k.even() {
term = -term;
}
let s_next = &s + &term; let (lo, hi) = if negative {
let bound = (&pow * eps) / (Rational::from(k + 1) * (Rational::ONE + eps.clone()));
let bound = -bound.abs();
(&s_next + bound, s_next.clone())
} else {
if s < s_next {
(s.clone(), s_next.clone())
} else {
(s_next.clone(), s.clone())
}
};
s = s_next;
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 ln_rational_helper(x: &Rational, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
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 ln_prec_round_normal(x_lo, prec, rm);
}
let (x_lo, x_hi) = floor_and_ceiling((x_lo, x_o));
let (ln_lo, mut o_lo) = ln_prec_round_normal(x_lo, prec, rm);
let (ln_hi, mut o_hi) = ln_prec_round_normal(x_hi, prec, rm);
if o_lo == Equal {
o_lo = o_hi;
}
if o_hi == Equal {
o_hi = o_lo;
}
if o_lo == o_hi && ln_lo == ln_hi {
return (ln_lo, o_lo);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
fn ln_rational_helper_extended(x: &Rational, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
let mut working_prec = prec + 10;
let mut increment = Limb::WIDTH;
loop {
let (x_lo, x_o) = ExtendedFloat::from_rational_prec_round_ref(x, working_prec, Floor);
if x_o == Equal {
return ln_prec_round_normal_extended(x_lo, prec, rm);
}
let (x_lo, x_hi) = crate::float::basic::extended::floor_and_ceiling((x_lo, x_o));
let (ln_lo, mut o_lo) = ln_prec_round_normal_extended(x_lo, prec, rm);
let (ln_hi, mut o_hi) = ln_prec_round_normal_extended(x_hi, prec, rm);
if o_lo == Equal {
o_lo = o_hi;
}
if o_hi == Equal {
o_hi = o_lo;
}
if o_lo == o_hi && ln_lo == ln_hi {
return (ln_lo, o_lo);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
#[allow(clippy::too_many_arguments)]
fn log_ui_s(
p: &mut [Integer],
b: &mut [Integer],
t: &mut [Integer],
q: &mut u64,
n1: u64,
n2: u64,
p_val: i64,
k: u64,
need_p: bool,
) {
if n2 == n1 + 1 {
p[0] = Integer::from(if n1 == 1 { p_val } else { -p_val });
*q = k;
b[0] = Integer::from(n1);
t[0] = p[0].clone();
} else {
let m = (n1 >> 1) + (n2 >> 1) + (n1 & n2 & 1);
log_ui_s(p, b, t, q, n1, m, p_val, k, true);
let mut q1 = 0;
let (p_head, p_tail) = p.split_first_mut().unwrap();
let (b_head, b_tail) = b.split_first_mut().unwrap();
let (t_head, t_tail) = t.split_first_mut().unwrap();
log_ui_s(p_tail, b_tail, t_tail, &mut q1, m, n2, p_val, k, need_p);
t_tail[0] *= &*p_head * &*b_head;
*t_head = ((&*t_head * &b_tail[0]) << q1) + &t_tail[0];
if need_p {
*p_head *= &p_tail[0];
}
*q += q1;
*b_head *= &b_tail[0];
}
}
fn ln_unsigned_prec_round_normal(n: u64, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
assert_ne!(rm, Exact, "Inexact ln");
let three_n = 3u128 * u128::from(n);
let k = u64::from(128 - three_n.leading_zeros() - 2);
let mut p = i64::exact_from(i128::from(n) - i128::power_of_2(k));
let mut kk = k;
if p != 0 {
let zeros = p.trailing_zeros();
p >>= zeros;
kk -= u64::from(zeros);
}
let mut w = prec + prec.ceiling_log_base_2() + 10;
loop {
let abs_p = p.unsigned_abs();
let n_terms = if abs_p == 0 {
2
} else {
let log2_abs_p = if abs_p == 1 {
0
} else {
abs_p.ceiling_log_base_2()
};
w.div_ceil(kk - log2_abs_p).max(2)
};
let integer_bits = n_terms.saturating_mul(n_terms.ceiling_log_base_2().saturating_add(kk));
if integer_bits.saturating_add(64) >= Float::MAX_EXPONENT_U64 {
return Float::from(n).ln_prec_round(prec, rm);
}
let lg_n = usize::exact_from(n_terms.ceiling_log_base_2() + 1);
let mut scratch = vec![Integer::ZERO; lg_n * 3];
split_into_chunks_mut!(scratch, lg_n, [p_arr, b_arr], t_arr);
let mut q0 = 0;
log_ui_s(p_arr, b_arr, t_arr, &mut q0, 1, n_terms, p, kk, false);
let t_num = Float::from_integer_prec(take(&mut t_arr[0]), w).0;
let t_den = Float::from_integer_prec(take(&mut b_arr[0]), w).0 << q0;
let t = t_num / t_den + Float::ln_2_prec(w).0 * Float::from_unsigned_prec(k, w).0;
let err = (k + 6).ceiling_log_base_2() + 1;
if float_can_round(
t.significand_ref().unwrap(),
w.saturating_sub(err),
prec,
rm,
) {
return Float::from_float_prec_round(t, prec, rm);
}
w += w >> 1;
}
}
impl Float {
#[inline]
pub fn ln_prec_round(self, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
assert_ne!(prec, 0);
match self {
Self(NaN | Infinity { sign: false } | Finite { sign: false, .. }) => {
(float_nan!(), Equal)
}
float_either_zero!() => (float_negative_infinity!(), Equal),
float_infinity!() => (float_infinity!(), Equal),
_ => ln_prec_round_normal(self, prec, rm),
}
}
#[inline]
pub fn ln_prec_round_ref(&self, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
assert_ne!(prec, 0);
match self {
Self(NaN | Infinity { sign: false } | Finite { sign: false, .. }) => {
(float_nan!(), Equal)
}
float_either_zero!() => (float_negative_infinity!(), Equal),
float_infinity!() => (float_infinity!(), Equal),
_ => ln_prec_round_normal_ref(self, prec, rm),
}
}
#[inline]
pub fn ln_prec(self, prec: u64) -> (Self, Ordering) {
self.ln_prec_round(prec, Nearest)
}
#[inline]
pub fn ln_prec_ref(&self, prec: u64) -> (Self, Ordering) {
self.ln_prec_round_ref(prec, Nearest)
}
#[inline]
pub fn ln_round(self, rm: RoundingMode) -> (Self, Ordering) {
let prec = self.significant_bits();
self.ln_prec_round(prec, rm)
}
#[inline]
pub fn ln_round_ref(&self, rm: RoundingMode) -> (Self, Ordering) {
let prec = self.significant_bits();
self.ln_prec_round_ref(prec, rm)
}
#[inline]
pub fn ln_prec_round_assign(&mut self, prec: u64, rm: RoundingMode) -> Ordering {
let mut x = Self::ZERO;
swap(self, &mut x);
let o;
(*self, o) = x.ln_prec_round(prec, rm);
o
}
#[inline]
pub fn ln_prec_assign(&mut self, prec: u64) -> Ordering {
self.ln_prec_round_assign(prec, Nearest)
}
#[inline]
pub fn ln_round_assign(&mut self, rm: RoundingMode) -> Ordering {
let prec = self.significant_bits();
self.ln_prec_round_assign(prec, rm)
}
#[allow(clippy::needless_pass_by_value)]
#[inline]
pub fn ln_rational_prec_round(x: Rational, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
Self::ln_rational_prec_round_ref(&x, prec, rm)
}
pub fn ln_rational_prec_round_ref(
x: &Rational,
prec: u64,
rm: RoundingMode,
) -> (Self, Ordering) {
assert_ne!(prec, 0);
match x.sign() {
Equal => return (float_negative_infinity!(), Equal),
Less => return (float_nan!(), Equal),
Greater => {}
}
if *x == 1u32 {
return (float_zero!(), Equal);
}
assert_ne!(rm, Exact, "Inexact ln");
let eps = x - Rational::ONE;
if eps.floor_log_base_2_abs() <= i64::from(Self::MIN_EXPONENT) + 4 {
return ln_rational_near_one(&eps, prec, rm);
}
let x_exp = i32::saturating_from(x.floor_log_base_2_abs()).saturating_add(1);
if x_exp >= const { Self::MAX_EXPONENT - 1 } || x_exp <= const { Self::MIN_EXPONENT + 1 } {
ln_rational_helper_extended(x, prec, rm)
} else {
ln_rational_helper(x, prec, rm)
}
}
#[inline]
pub fn ln_rational_prec(x: Rational, prec: u64) -> (Self, Ordering) {
Self::ln_rational_prec_round(x, prec, Nearest)
}
#[inline]
pub fn ln_rational_prec_ref(x: &Rational, prec: u64) -> (Self, Ordering) {
Self::ln_rational_prec_round_ref(x, prec, Nearest)
}
pub fn ln_unsigned_prec_round(n: u64, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
assert_ne!(prec, 0);
match n {
0 => (Self::NEGATIVE_INFINITY, Equal),
1 => (Self::ZERO, Equal),
2 => Self::ln_2_prec_round(prec, rm),
_ => ln_unsigned_prec_round_normal(n, prec, rm),
}
}
#[inline]
pub fn ln_unsigned_prec(n: u64, prec: u64) -> (Self, Ordering) {
Self::ln_unsigned_prec_round(n, prec, Nearest)
}
}
impl Ln for Float {
type Output = Self;
#[inline]
fn ln(self) -> Self {
let prec = self.significant_bits();
self.ln_prec_round(prec, Nearest).0
}
}
impl Ln for &Float {
type Output = Float;
#[inline]
fn ln(self) -> Float {
let prec = self.significant_bits();
self.ln_prec_round_ref(prec, Nearest).0
}
}
impl LnAssign for Float {
#[inline]
fn ln_assign(&mut self) {
let prec = self.significant_bits();
self.ln_prec_round_assign(prec, Nearest);
}
}
#[inline]
#[allow(clippy::type_repetition_in_bounds)]
pub fn primitive_float_ln<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::ln_prec, x)
}
#[inline]
#[allow(clippy::type_repetition_in_bounds)]
pub fn primitive_float_ln_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::ln_rational_prec_ref, x)
}