use crate::InnerFloat::{Finite, Infinity, NaN, Zero};
use crate::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,
};
use core::cmp::Ordering::{self, *};
use core::mem::swap;
use malachite_base::num::arithmetic::traits::{
Agm, CeilingLogBase2, IsPowerOf2, Ln, LnAssign, Sign,
};
use malachite_base::num::basic::floats::PrimitiveFloat;
use malachite_base::num::basic::integers::PrimitiveInt;
use malachite_base::num::basic::traits::{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::natural::arithmetic::float_extras::float_can_round;
use malachite_nz::platform::Limb;
use malachite_q::Rational;
fn ln_prec_round_normal_ref(x: &Float, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
if *x == 1u32 {
return (Float::ZERO, Equal);
}
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);
}
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_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 mut x_hi = x_lo.clone();
x_hi.increment();
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 mut x_hi = x_lo.clone();
x_hi.increment();
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;
}
}
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 x_exp = i32::saturating_from(x.floor_log_base_2_abs()).saturating_add(1);
if x_exp >= Self::MAX_EXPONENT - 1 || x_exp <= 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)
}
}
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)
}