use crate::InnerFloat::{Finite, Infinity, NaN, Zero};
use crate::float::arithmetic::ln::{SliverOfOne, sliver_of_one};
use crate::float::arithmetic::log_base_2::extended_log_base_2_of_rational;
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,
};
use alloc::vec::Vec;
use core::cmp::Ordering::{self, *};
use malachite_base::num::arithmetic::traits::{
CeilingLogBase2, CheckedLogBase, DivisibleBy, Gcd, IsPowerOf2, LogBase, LogBaseAssign, Mod,
ModAdd, ModMul, ModPow, Pow, Sign,
};
use malachite_base::num::basic::floats::PrimitiveFloat;
use malachite_base::num::basic::integers::PrimitiveInt;
use malachite_base::num::basic::traits::{One, Two, Zero as ZeroTrait};
use malachite_base::num::conversion::traits::{ExactFrom, RoundingFrom};
use malachite_base::num::factorization::traits::ExpressAsPower;
use malachite_base::num::logic::traits::SignificantBits;
use malachite_base::rounding_modes::RoundingMode::{self, *};
use malachite_nz::integer::Integer;
use malachite_nz::natural::Natural;
use malachite_nz::natural::arithmetic::float_extras::float_can_round;
use malachite_nz::platform::Limb;
use malachite_q::Rational;
pub(crate) fn odd_significand_and_exponent(x: &Float) -> (Natural, i64) {
let sig = x.significand_ref().unwrap();
let strip = sig.trailing_zeros().unwrap();
let s = sig >> strip;
let t = i64::from(x.get_exponent().unwrap()) - i64::exact_from(sig.significant_bits())
+ i64::exact_from(strip);
(s, t)
}
pub(crate) fn dyadic_log_of_root(s: &Natural, t: i64, z: i64, h: &Natural) -> Option<i64> {
if *h == 1u32 {
return if *s == 1u32 && z != 0 && t.divisible_by(z) {
Some(t / z)
} else {
None
};
}
if *s == 1u32 {
return None;
}
let m = i64::exact_from(s.checked_log_base(h)?);
if z.checked_mul(m)? == t {
Some(m)
} else {
None
}
}
pub(crate) fn dyadic_primitive_root(s_b: &Natural, t_b: i64) -> (i64, Natural, u64) {
if *s_b == 1u32 {
return if t_b > 0 {
(1, Natural::ONE, u64::exact_from(t_b))
} else {
(-1, Natural::ONE, u64::exact_from(-t_b))
};
}
let (h0, k) = s_b.express_as_power().unwrap_or_else(|| (s_b.clone(), 1));
let e = if t_b == 0 {
k
} else {
k.gcd(t_b.unsigned_abs())
};
(t_b / i64::exact_from(e), h0.pow(k / e), e)
}
pub(crate) fn rational_root_parts(g: &Rational) -> (i64, Natural, Natural) {
let num = g.numerator_ref();
let den = g.denominator_ref();
let num_z = num.trailing_zeros().unwrap();
let den_z = den.trailing_zeros().unwrap();
(
i64::exact_from(num_z) - i64::exact_from(den_z),
num >> num_z,
den >> den_z,
)
}
pub(crate) fn dyadic_log_of_rational_root(s: &Natural, t: i64, g: &Rational) -> Option<i64> {
let (z, hn, hd) = rational_root_parts(g);
if hd == 1u32 {
dyadic_log_of_root(s, t, z, &hn)
} else if hn == 1u32 {
dyadic_log_of_root(s, t, -z, &hd).map(|m| -m)
} else {
None
}
}
pub(crate) fn rational_value_log_of_dyadic_root(x: &Rational, z: i64, h: &Natural) -> Option<i64> {
let num = x.numerator_ref();
let den = x.denominator_ref();
let num_z = num.trailing_zeros().unwrap();
let den_z = den.trailing_zeros().unwrap();
let v2 = i64::exact_from(num_z) - i64::exact_from(den_z);
let on = num >> num_z;
let od = den >> den_z;
if *h == 1u32 {
return if on == 1u32 && od == 1u32 && z != 0 && v2.divisible_by(z) {
Some(v2 / z)
} else {
None
};
}
if on != 1u32 && od == 1u32 {
let m = i64::exact_from(on.checked_log_base(h)?);
if z.checked_mul(m)? == v2 {
Some(m)
} else {
None
}
} else if on == 1u32 && od != 1u32 {
let m = i64::exact_from(od.checked_log_base(h)?);
if z.checked_mul(m)? == -v2 {
Some(-m)
} else {
None
}
} else {
None
}
}
const FILTER_PRIME: u64 = 0xFFFFFFFFFFFFFFC5;
fn implicit_sum_is_pow(high: &Natural, shift: u64, low: &Integer, h: &Natural, m: u64) -> bool {
let p = Natural::from(FILTER_PRIME);
let lhs = (high % &p)
.mod_mul(Natural::TWO.mod_pow(Natural::from(shift), &p), &p)
.mod_add(Natural::exact_from(low.mod_op(Integer::from(&p))), &p);
if (h % &p).mod_pow(Natural::from(m), &p) != lhs {
return false;
}
Integer::from(high << shift) + low == h.pow(m)
}
pub(crate) fn dyadic_1p_log_of_root(x: &Float, z: i64, h: &Natural) -> Option<i64> {
let (s, t) = odd_significand_and_exponent(x);
let neg = *x < 0u32;
if t >= 0 {
debug_assert!(!neg);
if t == 0 {
let n = &s + Natural::ONE;
let r = n.trailing_zeros().unwrap();
return dyadic_log_of_root(&(n >> r), i64::exact_from(r), z, h);
}
if z != 0 || *h == 1u32 {
return None;
}
let l = u64::exact_from(t) + s.significant_bits();
for m in pow_bit_length_candidates(h, l) {
if implicit_sum_is_pow(&s, u64::exact_from(t), &Integer::ONE, h, m) {
return Some(i64::exact_from(m));
}
}
return None;
}
let t_abs = u64::exact_from(-t);
if *h == 1u32 {
return if neg
&& z != 0
&& t.divisible_by(z)
&& s.significant_bits() == t_abs
&& (&s + Natural::ONE).is_power_of_2()
{
Some(t / z)
} else {
None
};
}
if z == 0 || !t.divisible_by(z) {
return None;
}
let m = t / z;
if m <= 0 {
return None;
}
let low = if neg {
-Integer::from(&s)
} else {
Integer::from(&s)
};
if implicit_sum_is_pow(&Natural::ONE, t_abs, &low, h, u64::exact_from(m)) {
Some(m)
} else {
None
}
}
fn pow_bit_length_candidates(h: &Natural, l: u64) -> Vec<u64> {
let h_float = Float::exact_from(h.clone());
let lo = Rational::exact_from(h_float.log_base_2_prec_round_ref(80, Floor).0);
let hi = Rational::exact_from(h_float.log_base_2_prec_round_ref(80, Ceiling).0);
let m_min = Integer::rounding_from(Rational::from(l - 1) / hi, Ceiling).0;
let m_max = Integer::rounding_from(Rational::from(l) / lo, Floor).0;
let mut candidates = Vec::new();
let mut m = m_min;
while m <= m_max && candidates.len() < 4 {
if m > 0u32 {
candidates.push(u64::exact_from(&m));
}
m += Integer::ONE;
}
candidates
}
pub(crate) fn rational_log_base(x: &Float, base: u64) -> Option<Rational> {
let (g, e_base) = base.express_as_power().unwrap_or((base, 1));
let z = i64::exact_from(g.trailing_zeros());
let h = Natural::from(g >> g.trailing_zeros());
let (s, t) = odd_significand_and_exponent(x);
let m = dyadic_log_of_root(&s, t, z, &h)?;
Some(Rational::from_signeds(m, i64::exact_from(e_base)))
}
pub(crate) fn rational_log_base_of_rational(x: &Rational, base: u64) -> Option<Rational> {
let (g, e_base) = base.express_as_power().unwrap_or((base, 1));
x.checked_log_base(g)
.map(|m| Rational::from_signeds(m, i64::exact_from(e_base)))
}
fn log_base_prec_round_normal(
x: &Float,
base: u64,
prec: u64,
rm: RoundingMode,
) -> (Float, Ordering) {
if *x == 1u32 {
return (Float::ZERO, Equal);
}
if let Some(q) = rational_log_base(x, base) {
return Float::from_rational_prec_round(q, prec, rm);
}
match sliver_of_one(x) {
SliverOfOne::Representable(d) => return d.log_base_1_plus_x_prec_round(base, prec, rm),
SliverOfOne::Underflow => {
return Float::log_base_rational_prec_round(Rational::exact_from(x), base, prec, rm);
}
SliverOfOne::No => {}
}
assert_ne!(rm, Exact, "Inexact log_base");
let base_float = Float::from(base);
let mut working_prec = prec + 4 + prec.ceiling_log_base_2();
let mut increment = Limb::WIDTH;
loop {
let t = x
.ln_prec_ref(working_prec)
.0
.div_prec(base_float.ln_prec_ref(working_prec).0, working_prec)
.0;
if float_can_round(t.significand_ref().unwrap(), working_prec - 4, prec, rm) {
return Float::from_float_prec_round(t, prec, rm);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
fn log_base_rational_prec_round_helper(
x: &Rational,
base: u64,
prec: u64,
rm: RoundingMode,
) -> (Float, Ordering) {
let base_float = Float::from(base);
let mut working_prec = prec + 6 + prec.ceiling_log_base_2();
let mut increment = Limb::WIDTH;
loop {
let num = extended_log_base_2_of_rational(x, working_prec);
let den = ExtendedFloat::from(base_float.log_base_2_prec_ref(working_prec).0);
let quotient = num.div_prec_val_ref(&den, working_prec).0;
if float_can_round(
quotient.x.significand_ref().unwrap(),
working_prec - 6,
prec,
rm,
) {
let (rounded, o) = Float::from_float_prec_round(quotient.x, prec, rm);
let mut result = ExtendedFloat::from(rounded);
result.exp = result.exp.checked_add(quotient.exp).unwrap();
return result.into_float_helper(prec, rm, o);
}
working_prec += increment;
increment = working_prec >> 1;
}
}
impl Float {
#[inline]
pub fn log_base_prec_round(self, base: u64, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
assert_ne!(prec, 0);
assert!(base > 1, "Logarithm base must be greater than 1");
if base.is_power_of_2() {
return self.log_base_power_of_2_prec_round(i64::from(base.trailing_zeros()), prec, rm);
}
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),
_ => log_base_prec_round_normal(&self, base, prec, rm),
}
}
#[inline]
pub fn log_base_prec_round_ref(
&self,
base: u64,
prec: u64,
rm: RoundingMode,
) -> (Self, Ordering) {
assert_ne!(prec, 0);
assert!(base > 1, "Logarithm base must be greater than 1");
if base.is_power_of_2() {
return self.log_base_power_of_2_prec_round_ref(
i64::from(base.trailing_zeros()),
prec,
rm,
);
}
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),
_ => log_base_prec_round_normal(self, base, prec, rm),
}
}
#[inline]
pub fn log_base_prec(self, base: u64, prec: u64) -> (Self, Ordering) {
self.log_base_prec_round(base, prec, Nearest)
}
#[inline]
pub fn log_base_prec_ref(&self, base: u64, prec: u64) -> (Self, Ordering) {
self.log_base_prec_round_ref(base, prec, Nearest)
}
#[inline]
pub fn log_base_round(self, base: u64, rm: RoundingMode) -> (Self, Ordering) {
let prec = self.significant_bits();
self.log_base_prec_round(base, prec, rm)
}
#[inline]
pub fn log_base_round_ref(&self, base: u64, rm: RoundingMode) -> (Self, Ordering) {
self.log_base_prec_round_ref(base, self.significant_bits(), rm)
}
#[inline]
pub fn log_base_prec_round_assign(
&mut self,
base: u64,
prec: u64,
rm: RoundingMode,
) -> Ordering {
let (result, o) = core::mem::take(self).log_base_prec_round(base, prec, rm);
*self = result;
o
}
#[inline]
pub fn log_base_prec_assign(&mut self, base: u64, prec: u64) -> Ordering {
self.log_base_prec_round_assign(base, prec, Nearest)
}
#[inline]
pub fn log_base_round_assign(&mut self, base: u64, rm: RoundingMode) -> Ordering {
let prec = self.significant_bits();
self.log_base_prec_round_assign(base, prec, rm)
}
#[allow(clippy::needless_pass_by_value)]
#[inline]
pub fn log_base_rational_prec_round(
x: Rational,
base: u64,
prec: u64,
rm: RoundingMode,
) -> (Self, Ordering) {
Self::log_base_rational_prec_round_ref(&x, base, prec, rm)
}
pub fn log_base_rational_prec_round_ref(
x: &Rational,
base: u64,
prec: u64,
rm: RoundingMode,
) -> (Self, Ordering) {
assert_ne!(prec, 0);
assert!(base > 1, "Logarithm base must be greater than 1");
if base.is_power_of_2() {
return Self::log_base_power_of_2_rational_prec_round_ref(
x,
i64::from(base.trailing_zeros()),
prec,
rm,
);
}
match x.sign() {
Equal => return (float_negative_infinity!(), Equal),
Less => return (float_nan!(), Equal),
Greater => {}
}
if let Some(q) = rational_log_base_of_rational(x, base) {
return Self::from_rational_prec_round(q, prec, rm);
}
assert_ne!(rm, Exact, "Inexact log_base");
log_base_rational_prec_round_helper(x, base, prec, rm)
}
#[inline]
pub fn log_base_rational_prec(x: Rational, base: u64, prec: u64) -> (Self, Ordering) {
Self::log_base_rational_prec_round(x, base, prec, Nearest)
}
#[inline]
pub fn log_base_rational_prec_ref(x: &Rational, base: u64, prec: u64) -> (Self, Ordering) {
Self::log_base_rational_prec_round_ref(x, base, prec, Nearest)
}
}
impl LogBase<u64> for Float {
type Output = Self;
#[inline]
fn log_base(self, base: u64) -> Self {
let prec = self.significant_bits();
self.log_base_prec_round(base, prec, Nearest).0
}
}
impl LogBase<u64> for &Float {
type Output = Float;
#[inline]
fn log_base(self, base: u64) -> Float {
self.log_base_prec_round_ref(base, self.significant_bits(), Nearest)
.0
}
}
impl LogBaseAssign<u64> for Float {
#[inline]
fn log_base_assign(&mut self, base: u64) {
let prec = self.significant_bits();
self.log_base_prec_round_assign(base, prec, Nearest);
}
}
#[inline]
#[allow(clippy::type_repetition_in_bounds)]
pub fn primitive_float_log_base<T: PrimitiveFloat>(x: T, base: u64) -> T
where
Float: From<T> + PartialOrd<T>,
for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
{
emulate_float_to_float_fn(|x, prec| Float::log_base_prec(x, base, prec), x)
}
#[inline]
#[allow(clippy::type_repetition_in_bounds)]
pub fn primitive_float_log_base_rational<T: PrimitiveFloat>(x: &Rational, base: u64) -> T
where
Float: PartialOrd<T>,
for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
{
emulate_rational_to_float_fn(
|x, prec| Float::log_base_rational_prec_ref(x, base, prec),
x,
)
}