malachite_float/float/arithmetic/log_base_float_base_1_plus_x.rs
1// Copyright © 2026 Mikhail Hogrefe
2//
3// This file is part of Malachite.
4//
5// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
6// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
7// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
8
9use crate::InnerFloat::{Infinity, NaN};
10use crate::float::arithmetic::log_base::{
11 dyadic_1p_log_of_root, dyadic_primitive_root, odd_significand_and_exponent,
12};
13use crate::float::basic::extended::ExtendedFloat;
14use crate::{
15 Float, emulate_float_float_to_float_fn, float_infinity, float_nan, float_negative_infinity,
16};
17use core::cmp::Ordering::{self, *};
18use malachite_base::num::arithmetic::traits::{
19 CeilingLogBase2, LogBaseOf1PlusX, LogBaseOf1PlusXAssign,
20};
21use malachite_base::num::basic::floats::PrimitiveFloat;
22use malachite_base::num::basic::integers::PrimitiveInt;
23use malachite_base::num::basic::traits::{NegativeZero, Zero as ZeroTrait};
24use malachite_base::num::conversion::traits::{ExactFrom, RoundingFrom};
25use malachite_base::num::logic::traits::SignificantBits;
26use malachite_base::rounding_modes::RoundingMode::{self, *};
27use malachite_nz::natural::arithmetic::float::round::float_can_round;
28use malachite_nz::platform::Limb;
29use malachite_q::Rational;
30
31// Returns `Some(log_base(1 + x))` when it is rational, and `None` when it is irrational. `x` must
32// be a finite [`Float`] in (-1, 0) or positive (not 0), and `base` a finite positive [`Float`] not
33// equal to 1.
34//
35// `log_base(1 + x)` is rational exactly when `1 + x` and `base` are commensurable. `1 + x` is
36// formed exactly as a `Rational`, and `base` is dyadic, so this reuses
37// `rational_log_base_rational_rational_base`; for a base in (0, 1) -- where
38// `Rational::checked_log_base` requires a base above 1 -- the identity `log_b(y) = -log_{1/b}(y)`
39// reduces to a base above 1. Balloon-safe via the `64 * prec` size bound on `x`'s exponent and the
40// operands' precisions (an `x` near -1 has a near-zero exponent and is materialized, but `1 + x` is
41// then bounded by `x`'s precision).
42pub(crate) fn log_base_float_base_1_plus_x_rational(x: &Float, base: &Float) -> Option<Rational> {
43 if *x == 0u32 {
44 return Some(Rational::ZERO);
45 }
46 // The base's primitive root comes from its odd significand and exponent without materializing
47 // it, and `1 + x` is matched against that root implicitly (see `dyadic_1p_log_of_root`): its
48 // integer form can be enormous even when `x` has few bits, so it is only materialized to verify
49 // a match that a cheap congruence filter has already endorsed. No size cutoff: skipping the
50 // check when the result is exactly representable would leave the Ziv loop unable to terminate.
51 let (s_b, t_b) = odd_significand_and_exponent(base);
52 let (z, h, e_base) = dyadic_primitive_root(&s_b, t_b);
53 let m = dyadic_1p_log_of_root(x, z, &h)?;
54 Some(Rational::from_signeds(m, i64::exact_from(e_base)))
55}
56
57// The computation of log_base(1 + x) for a `Float` base is done by log_base(1 + x) = log_2(1 + x) /
58// log_2(base). The inputs are a finite `Float` `x` in (-1, 0) or positive (not 0), and a finite
59// positive `Float` base not equal to 1.
60//
61// Both logs are ordinary native `Float`s. `log_2(1 + x)` is routed through `log_base_2_1_plus_x` to
62// preserve accuracy for x near 0; it cannot underflow (|x| is at least the smallest positive
63// `Float`). Their quotient can overflow (base near 1) or underflow (x near 0), so the operands are
64// wrapped as `ExtendedFloat`s, divided in the extended range, and converted back with a single
65// `into_float_helper` clamp. A base in (0, 1) gives a negative `log_2(base)`, so the division
66// yields the (sign-flipped) result for free.
67fn log_base_float_base_1_plus_x_normal(
68 x: &Float,
69 base: &Float,
70 prec: u64,
71 rm: RoundingMode,
72) -> (Float, Ordering) {
73 // If log_base(1 + x) is rational -- 1 + x and base commensurable -- compute it directly.
74 if let Some(q) = log_base_float_base_1_plus_x_rational(x, base) {
75 return Float::from_rational_prec_round(q, prec, rm);
76 }
77 // The result is irrational, so it is never exactly representable.
78 assert_ne!(rm, Exact, "Inexact log_base_float_base_1_plus_x");
79 // The initial slack keeps working_prec at least 7, so the working_prec - 6 below stays
80 // positive.
81 let mut working_prec = prec + 6 + prec.ceiling_log_base_2();
82 let mut increment = Limb::WIDTH;
83 loop {
84 // log_2(1 + x) and log_2(base), correctly rounded and wrapped; both finite and nonzero (x
85 // is not 0 and base is not 1), neither underflowing.
86 let num = ExtendedFloat::from(x.log_base_2_1_plus_x_prec_ref(working_prec).0);
87 let den = ExtendedFloat::from(base.log_base_2_prec_ref(working_prec).0);
88 // log_2(1 + x) / log_2(base) in the extended range; cannot overflow or underflow here.
89 let quotient = num.div_prec_val_ref(&den, working_prec).0;
90 // Two correctly-rounded logs (<= 1/2 ulp each) and the division (<= 1/2 ulp) give under 2
91 // ulps total; working_prec - 6 correct bits comfortably suffice for the rounding test.
92 if float_can_round(
93 quotient.x.significand_ref().unwrap(),
94 working_prec - 6,
95 prec,
96 rm,
97 ) {
98 // Round the mantissa to prec, then place the extended exponent, clamping once to the
99 // Float range as the rounding mode dictates.
100 let (rounded, o) = Float::from_float_prec_round(quotient.x, prec, rm);
101 let mut result = ExtendedFloat::from(rounded);
102 result.exp = result.exp.checked_add(quotient.exp).unwrap();
103 return result.into_float_helper(prec, rm, o);
104 }
105 // Increase the precision.
106 working_prec += increment;
107 increment = working_prec >> 1;
108 }
109}
110
111// Computes log_base(1 + x) = ln(1 + x) / ln(base) for `Float` `x` and `base`, following IEEE
112// division of the natural logs for every special case (so the function is total: no input value
113// panics). `ln(1 + x)` uses the sign-preserving `ln_1p` convention, so `ln_1p(x)` has the sign of
114// `x` for `x` in (-1, infinity].
115fn log_base_float_base_1_plus_x_helper(
116 x: &Float,
117 base: &Float,
118 prec: u64,
119 rm: RoundingMode,
120) -> (Float, Ordering) {
121 // ln(1 + x) or ln(base) is NaN: x or base is NaN, base is negative, or 1 + x < 0.
122 if x.is_nan() || base.is_nan() {
123 return (float_nan!(), Equal);
124 }
125 if *base < 0u32 {
126 return (float_nan!(), Equal); // base negative finite or -infinity
127 }
128 if *x < -1i32 {
129 return (float_nan!(), Equal); // 1 + x < 0 (including x = -infinity)
130 }
131 // x is in [-1, infinity] and not NaN; base is +infinity, zero, or positive finite. ln_1p(x) has
132 // the sign of x (negative for x in [-1, 0) including -0.0).
133 let x_neg = x.is_sign_negative();
134 if base.is_infinite() {
135 // ln(base) = +infinity. ln_1p(x) / +infinity = 0 for finite ln_1p(x) (NaN when it is
136 // +-infinity, i.e. x = +infinity or x = -1), with the sign of ln_1p(x).
137 if x.is_infinite() || *x == -1i32 {
138 return (float_nan!(), Equal);
139 }
140 return if x_neg {
141 (Float::NEGATIVE_ZERO, Equal)
142 } else {
143 (Float::ZERO, Equal)
144 };
145 }
146 if *base == 0u32 {
147 // ln(base) = -infinity. Sign-flipped from the +infinity case.
148 if x.is_infinite() || *x == -1i32 {
149 return (float_nan!(), Equal);
150 }
151 return if x_neg {
152 (Float::ZERO, Equal)
153 } else {
154 (Float::NEGATIVE_ZERO, Equal)
155 };
156 }
157 if *base == 1u32 {
158 // ln(base) = +0. ln_1p(x) / +0 = +-infinity by the sign of ln_1p(x), or NaN for ln_1p(x) =
159 // +-0 (x = +-0).
160 if *x == 0u32 {
161 return (float_nan!(), Equal);
162 }
163 return if x_neg {
164 (float_negative_infinity!(), Equal)
165 } else {
166 (float_infinity!(), Equal)
167 };
168 }
169 // base is positive finite and not 1.
170 if x.is_infinite() {
171 // ln_1p(+infinity) = +infinity. +infinity / ln(base): +infinity for base > 1, -infinity for
172 // base < 1.
173 return if *base < 1u32 {
174 (float_negative_infinity!(), Equal)
175 } else {
176 (float_infinity!(), Equal)
177 };
178 }
179 if *x == -1i32 {
180 // ln_1p(-1) = ln(0) = -infinity. -infinity / ln(base): -infinity for base > 1, +infinity
181 // for base < 1.
182 return if *base < 1u32 {
183 (float_infinity!(), Equal)
184 } else {
185 (float_negative_infinity!(), Equal)
186 };
187 }
188 if *x == 0u32 {
189 // ln_1p(+-0) = +-0. +-0 / ln(base) = 0, with the sign of ln_1p(x) times the sign of
190 // ln(base) (positive for base > 1, negative for base < 1).
191 return if x_neg == (*base < 1u32) {
192 (Float::ZERO, Equal)
193 } else {
194 (Float::NEGATIVE_ZERO, Equal)
195 };
196 }
197 // x is finite in (-1, 0) or positive (not 0), and base is positive finite and not 1.
198 log_base_float_base_1_plus_x_normal(x, base, prec, rm)
199}
200
201impl Float {
202 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
203 /// to the specified precision and with the specified rounding mode. The [`Float`] is taken by
204 /// value and the base by reference. An [`Ordering`] is also returned, indicating whether the
205 /// rounded value is less than, equal to, or greater than the exact value. Although `NaN`s are
206 /// not comparable to any [`Float`], whenever this function returns a `NaN` it also returns
207 /// `Equal`.
208 ///
209 /// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned. Otherwise the
210 /// base may be any [`Float`]: the function is defined as $\ln(1+x) / \ln b$ for every pair,
211 /// applying IEEE division to the natural logs, and never panics on an input value. In
212 /// particular a base in $(0,1)$ gives a (sign-flipped) logarithm, and the non-normal and
213 /// degenerate bases follow the limits below.
214 ///
215 /// This computes $\log_2(1+x) / \log_2 b$, routing through
216 /// [`Float::log_base_2_1_plus_x_prec_ref`] to preserve accuracy for $x$ near 0, and wrapping
217 /// the quotient so it may overflow (base near 1) or underflow (x near 0) and be clamped exactly
218 /// once.
219 ///
220 /// See [`RoundingMode`] for a description of the possible rounding modes.
221 ///
222 /// $$
223 /// f(x,b,p,m) = \log_b(1+x)+\varepsilon.
224 /// $$
225 /// - If $\log_b(1+x)$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to
226 /// be 0.
227 /// - If $\log_b(1+x)$ is finite and nonzero, and $m$ is not `Nearest`, then $|\varepsilon| <
228 /// 2^{\lfloor\log_2 |\log_b(1+x)|\rfloor-p+1}$.
229 /// - If $\log_b(1+x)$ is finite and nonzero, and $m$ is `Nearest`, then $|\varepsilon| \leq
230 /// 2^{\lfloor\log_2 |\log_b(1+x)|\rfloor-p}$.
231 ///
232 /// If the output has a precision, it is `prec`.
233 ///
234 /// Special cases (with $b$ the base):
235 /// - $f(\text{NaN},b,p,m)=\text{NaN}$, and $f(x,\text{NaN},p,m)=\text{NaN}$
236 /// - $f(x,b,p,m)=\text{NaN}$ for $x<-1$ or $b<0$
237 /// - $f(\infty,b,p,m)=\infty$ for $b>1$, and $-\infty$ for $0\leq b<1$
238 /// - $f(-1.0,b,p,m)=-\infty$ for $b>1$, and $\infty$ for $0<b<1$
239 /// - $f(\pm0.0,b,p,m)=0$ (the sign of $\pm0.0$ times the sign of $1/\ln b$)
240 /// - $f(x,\infty,p,m)=0$ for finite $x>-1$ with $x\neq0$ (and $\text{NaN}$ for
241 /// $x\in\{\infty,-1\}$)
242 /// - $f(x,\pm0.0,p,m)=0$ for finite $x>-1$ with $x\neq0$ (and $\text{NaN}$ for
243 /// $x\in\{\infty,-1\}$)
244 /// - $f(x,1.0,p,m)=\infty$ for $x>0$ or $x=\infty$, $-\infty$ for $-1\leq x<0$, and
245 /// $\text{NaN}$ for $x=\pm0.0$
246 /// - $f(g^a-1,g^e,p,m)=a/e$ for a common rational $g$, rounded to precision $p$; the result is
247 /// exact if and only if $a/e$ is representable with precision $p$ (for example
248 /// $\log_4(1+1)=1/2$)
249 ///
250 /// This function can both overflow (for a base near 1) and underflow (for an $x$ near 0).
251 ///
252 /// # Worst-case complexity
253 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
254 ///
255 /// $M(n, m) = O(n \log n + m \log m)$
256 ///
257 /// where $T$ is time, $M$ is additional memory, $n$ is `prec`, and $m$ is
258 /// `max(self.significant_bits(), base.significant_bits())`.
259 ///
260 /// # Panics
261 /// Panics if `prec` is zero, or if `rm` is `Exact` but the result cannot be represented exactly
262 /// with the given precision.
263 ///
264 /// # Examples
265 /// ```
266 /// use malachite_base::rounding_modes::RoundingMode::*;
267 /// use malachite_float::Float;
268 /// use std::cmp::Ordering::*;
269 ///
270 /// let (log, o) =
271 /// Float::from(8).log_base_float_base_1_plus_x_prec_round(&Float::from(3), 10, Exact);
272 /// assert_eq!(log.to_string(), "2.0000"); // log_3(1 + 8) = log_3(9) = 2
273 /// assert_eq!(o, Equal);
274 ///
275 /// let (log, o) =
276 /// Float::from(3).log_base_float_base_1_plus_x_prec_round(&Float::from(0.5), 10, Exact);
277 /// assert_eq!(log.to_string(), "-2.0000"); // log_{1/2}(1 + 3) = log_{1/2}(4) = -2
278 /// assert_eq!(o, Equal);
279 /// ```
280 #[inline]
281 pub fn log_base_float_base_1_plus_x_prec_round(
282 self,
283 base: &Self,
284 prec: u64,
285 rm: RoundingMode,
286 ) -> (Self, Ordering) {
287 assert_ne!(prec, 0);
288 log_base_float_base_1_plus_x_helper(&self, base, prec, rm)
289 }
290
291 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
292 /// to the specified precision and with the specified rounding mode. Both are taken by
293 /// reference. An [`Ordering`] is also returned, indicating whether the rounded value is less
294 /// than, equal to, or greater than the exact value.
295 ///
296 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details, special cases, and a
297 /// description of the rounding behavior.
298 ///
299 /// # Worst-case complexity
300 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
301 ///
302 /// $M(n, m) = O(n \log n + m \log m)$
303 ///
304 /// where $T$ is time, $M$ is additional memory, $n$ is `prec`, and $m$ is
305 /// `max(self.significant_bits(), base.significant_bits())`.
306 ///
307 /// # Panics
308 /// Panics if `prec` is zero, or if `rm` is `Exact` but the result cannot be represented exactly
309 /// with the given precision.
310 ///
311 /// # Examples
312 /// ```
313 /// use malachite_base::num::basic::traits::{One, Two};
314 /// use malachite_base::rounding_modes::RoundingMode::*;
315 /// use malachite_float::Float;
316 /// use std::cmp::Ordering::*;
317 ///
318 /// let x = Float::from(7);
319 /// let (log, o) = x.log_base_float_base_1_plus_x_prec_round_ref(&Float::TWO, 10, Exact);
320 /// assert_eq!(log.to_string(), "3.0000"); // log_2(1 + 7) = log_2(8) = 3
321 /// assert_eq!(o, Equal);
322 ///
323 /// let x = Float::ONE;
324 /// let (log, o) = x.log_base_float_base_1_plus_x_prec_round_ref(&Float::from(3), 20, Floor);
325 /// assert_eq!(log.to_string(), "0.63092899"); // log_3(2), rounded down
326 /// assert_eq!(o, Less);
327 /// ```
328 pub fn log_base_float_base_1_plus_x_prec_round_ref(
329 &self,
330 base: &Self,
331 prec: u64,
332 rm: RoundingMode,
333 ) -> (Self, Ordering) {
334 assert_ne!(prec, 0);
335 log_base_float_base_1_plus_x_helper(self, base, prec, rm)
336 }
337
338 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
339 /// to the nearest value of the specified precision. The [`Float`] is taken by value and the
340 /// base by reference. An [`Ordering`] is also returned.
341 ///
342 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
343 ///
344 /// # Worst-case complexity
345 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
346 ///
347 /// $M(n, m) = O(n \log n + m \log m)$
348 ///
349 /// where $T$ is time, $M$ is additional memory, $n$ is `prec`, and $m$ is
350 /// `max(self.significant_bits(), base.significant_bits())`.
351 ///
352 /// # Panics
353 /// Panics if `prec` is zero.
354 ///
355 /// # Examples
356 /// ```
357 /// use malachite_float::Float;
358 /// use std::cmp::Ordering::*;
359 ///
360 /// let (log, o) = Float::from(8).log_base_float_base_1_plus_x_prec(&Float::from(3), 10);
361 /// assert_eq!(log.to_string(), "2.0000"); // log_3(1 + 8) = log_3(9) = 2
362 /// assert_eq!(o, Equal);
363 /// ```
364 #[inline]
365 pub fn log_base_float_base_1_plus_x_prec(self, base: &Self, prec: u64) -> (Self, Ordering) {
366 self.log_base_float_base_1_plus_x_prec_round(base, prec, Nearest)
367 }
368
369 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
370 /// to the nearest value of the specified precision. Both are taken by reference. An
371 /// [`Ordering`] is also returned.
372 ///
373 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
374 ///
375 /// # Worst-case complexity
376 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
377 ///
378 /// $M(n, m) = O(n \log n + m \log m)$
379 ///
380 /// where $T$ is time, $M$ is additional memory, $n$ is `prec`, and $m$ is
381 /// `max(self.significant_bits(), base.significant_bits())`.
382 ///
383 /// # Panics
384 /// Panics if `prec` is zero.
385 ///
386 /// # Examples
387 /// ```
388 /// use malachite_float::Float;
389 /// use std::cmp::Ordering::*;
390 ///
391 /// let (log, o) = (&Float::from(8)).log_base_float_base_1_plus_x_prec_ref(&Float::from(3), 10);
392 /// assert_eq!(log.to_string(), "2.0000"); // log_3(1 + 8) = log_3(9) = 2
393 /// assert_eq!(o, Equal);
394 /// ```
395 #[inline]
396 pub fn log_base_float_base_1_plus_x_prec_ref(
397 &self,
398 base: &Self,
399 prec: u64,
400 ) -> (Self, Ordering) {
401 self.log_base_float_base_1_plus_x_prec_round_ref(base, prec, Nearest)
402 }
403
404 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
405 /// to the precision of the input and with the specified rounding mode. The [`Float`] is taken
406 /// by value and the base by reference. An [`Ordering`] is also returned.
407 ///
408 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
409 ///
410 /// # Worst-case complexity
411 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
412 ///
413 /// $M(n, m) = O(n \log n + m \log m)$
414 ///
415 /// where $T$ is time, $M$ is additional memory, $n$ is the precision of the input, and $m$ is
416 /// `base.significant_bits()`.
417 ///
418 /// # Panics
419 /// Panics if `rm` is `Exact` but the result cannot be represented exactly with the input's
420 /// precision.
421 ///
422 /// # Examples
423 /// ```
424 /// use malachite_base::rounding_modes::RoundingMode::*;
425 /// use malachite_float::Float;
426 /// use std::cmp::Ordering::*;
427 ///
428 /// let (log, o) = Float::from(8).log_base_float_base_1_plus_x_round(&Float::from(3), Exact);
429 /// assert_eq!(log.to_string(), "2.0"); // log_3(1 + 8) = log_3(9) = 2
430 /// assert_eq!(o, Equal);
431 /// ```
432 #[inline]
433 pub fn log_base_float_base_1_plus_x_round(
434 self,
435 base: &Self,
436 rm: RoundingMode,
437 ) -> (Self, Ordering) {
438 let prec = self.significant_bits();
439 self.log_base_float_base_1_plus_x_prec_round(base, prec, rm)
440 }
441
442 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
443 /// to the precision of the input and with the specified rounding mode. Both are taken by
444 /// reference. An [`Ordering`] is also returned.
445 ///
446 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
447 ///
448 /// # Worst-case complexity
449 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
450 ///
451 /// $M(n, m) = O(n \log n + m \log m)$
452 ///
453 /// where $T$ is time, $M$ is additional memory, $n$ is the precision of the input, and $m$ is
454 /// `base.significant_bits()`.
455 ///
456 /// # Panics
457 /// Panics if `rm` is `Exact` but the result cannot be represented exactly with the input's
458 /// precision.
459 ///
460 /// # Examples
461 /// ```
462 /// use malachite_base::rounding_modes::RoundingMode::*;
463 /// use malachite_float::Float;
464 /// use std::cmp::Ordering::*;
465 ///
466 /// let (log, o) =
467 /// (&Float::from(8)).log_base_float_base_1_plus_x_round_ref(&Float::from(3), Exact);
468 /// assert_eq!(log.to_string(), "2.0"); // log_3(1 + 8) = log_3(9) = 2
469 /// assert_eq!(o, Equal);
470 /// ```
471 #[inline]
472 pub fn log_base_float_base_1_plus_x_round_ref(
473 &self,
474 base: &Self,
475 rm: RoundingMode,
476 ) -> (Self, Ordering) {
477 self.log_base_float_base_1_plus_x_prec_round_ref(base, self.significant_bits(), rm)
478 }
479
480 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, in place, rounding
481 /// the result to the specified precision and with the specified rounding mode. The base is
482 /// taken by reference. An [`Ordering`] is returned.
483 ///
484 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
485 ///
486 /// # Worst-case complexity
487 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
488 ///
489 /// $M(n, m) = O(n \log n + m \log m)$
490 ///
491 /// where $T$ is time, $M$ is additional memory, $n$ is `prec`, and $m$ is
492 /// `max(self.significant_bits(), base.significant_bits())`.
493 ///
494 /// # Panics
495 /// Panics if `prec` is zero, or if `rm` is `Exact` but the result cannot be represented exactly
496 /// with the given precision.
497 ///
498 /// # Examples
499 /// ```
500 /// use malachite_base::rounding_modes::RoundingMode::*;
501 /// use malachite_float::Float;
502 /// use std::cmp::Ordering::*;
503 ///
504 /// let mut x = Float::from(8);
505 /// assert_eq!(
506 /// x.log_base_float_base_1_plus_x_prec_round_assign(&Float::from(3), 10, Exact),
507 /// Equal
508 /// );
509 /// assert_eq!(x.to_string(), "2.0000"); // log_3(1 + 8) = log_3(9) = 2
510 /// ```
511 #[inline]
512 pub fn log_base_float_base_1_plus_x_prec_round_assign(
513 &mut self,
514 base: &Self,
515 prec: u64,
516 rm: RoundingMode,
517 ) -> Ordering {
518 let (result, o) =
519 core::mem::take(self).log_base_float_base_1_plus_x_prec_round(base, prec, rm);
520 *self = result;
521 o
522 }
523
524 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, in place, rounding
525 /// the result to the nearest value of the specified precision. The base is taken by reference.
526 /// An [`Ordering`] is returned.
527 ///
528 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
529 ///
530 /// # Worst-case complexity
531 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
532 ///
533 /// $M(n, m) = O(n \log n + m \log m)$
534 ///
535 /// where $T$ is time, $M$ is additional memory, $n$ is `prec`, and $m$ is
536 /// `max(self.significant_bits(), base.significant_bits())`.
537 ///
538 /// # Panics
539 /// Panics if `prec` is zero.
540 ///
541 /// # Examples
542 /// ```
543 /// use malachite_float::Float;
544 ///
545 /// let mut x = Float::from(8);
546 /// x.log_base_float_base_1_plus_x_prec_assign(&Float::from(3), 10);
547 /// assert_eq!(x.to_string(), "2.0000"); // log_3(1 + 8) = log_3(9) = 2
548 /// ```
549 #[inline]
550 pub fn log_base_float_base_1_plus_x_prec_assign(&mut self, base: &Self, prec: u64) -> Ordering {
551 self.log_base_float_base_1_plus_x_prec_round_assign(base, prec, Nearest)
552 }
553
554 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, in place, rounding
555 /// the result to the precision of the input and with the specified rounding mode. The base is
556 /// taken by reference. An [`Ordering`] is returned.
557 ///
558 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for details and special cases.
559 ///
560 /// # Worst-case complexity
561 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
562 ///
563 /// $M(n, m) = O(n \log n + m \log m)$
564 ///
565 /// where $T$ is time, $M$ is additional memory, $n$ is the precision of the input, and $m$ is
566 /// `base.significant_bits()`.
567 ///
568 /// # Panics
569 /// Panics if `rm` is `Exact` but the result cannot be represented exactly with the input's
570 /// precision.
571 ///
572 /// # Examples
573 /// ```
574 /// use malachite_base::rounding_modes::RoundingMode::*;
575 /// use malachite_float::Float;
576 ///
577 /// let mut x = Float::from(8);
578 /// x.log_base_float_base_1_plus_x_round_assign(&Float::from(3), Exact);
579 /// assert_eq!(x.to_string(), "2.0"); // log_3(1 + 8) = log_3(9) = 2
580 /// ```
581 #[inline]
582 pub fn log_base_float_base_1_plus_x_round_assign(
583 &mut self,
584 base: &Self,
585 rm: RoundingMode,
586 ) -> Ordering {
587 let prec = self.significant_bits();
588 self.log_base_float_base_1_plus_x_prec_round_assign(base, prec, rm)
589 }
590}
591
592impl LogBaseOf1PlusX<Self> for Float {
593 type Output = Self;
594
595 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
596 /// to the nearest value of the input's precision. Both are taken by value.
597 ///
598 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for special cases.
599 ///
600 /// # Worst-case complexity
601 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
602 ///
603 /// $M(n, m) = O(n \log n + m \log m)$
604 ///
605 /// where $T$ is time, $M$ is additional memory, $n$ is the precision of the input, and $m$ is
606 /// `base.significant_bits()`.
607 ///
608 /// # Examples
609 /// ```
610 /// use malachite_base::num::arithmetic::traits::LogBaseOf1PlusX;
611 /// use malachite_float::Float;
612 ///
613 /// // log_3(1 + 8) = log_3(9) = 2
614 /// assert_eq!(
615 /// Float::from(8).log_base_1_plus_x(Float::from(3)).to_string(),
616 /// "2.0"
617 /// );
618 /// ```
619 #[inline]
620 fn log_base_1_plus_x(self, base: Self) -> Self {
621 let prec = self.significant_bits();
622 self.log_base_float_base_1_plus_x_prec_round(&base, prec, Nearest)
623 .0
624 }
625}
626
627impl LogBaseOf1PlusX<&Float> for &Float {
628 type Output = Float;
629
630 /// Computes $\log_b(1+x)$, where $x$ and the base $b$ are both [`Float`]s, rounding the result
631 /// to the nearest value of the input's precision. Both are taken by reference.
632 ///
633 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for special cases.
634 ///
635 /// # Worst-case complexity
636 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
637 ///
638 /// $M(n, m) = O(n \log n + m \log m)$
639 ///
640 /// where $T$ is time, $M$ is additional memory, $n$ is the precision of the input, and $m$ is
641 /// `base.significant_bits()`.
642 ///
643 /// # Examples
644 /// ```
645 /// use malachite_base::num::arithmetic::traits::LogBaseOf1PlusX;
646 /// use malachite_float::Float;
647 ///
648 /// // log_3(1 + 8) = log_3(9) = 2
649 /// assert_eq!(
650 /// (&Float::from(8))
651 /// .log_base_1_plus_x(&Float::from(3))
652 /// .to_string(),
653 /// "2.0"
654 /// );
655 /// ```
656 #[inline]
657 fn log_base_1_plus_x(self, base: &Float) -> Float {
658 self.log_base_float_base_1_plus_x_prec_round_ref(base, self.significant_bits(), Nearest)
659 .0
660 }
661}
662
663impl LogBaseOf1PlusXAssign<&Self> for Float {
664 /// Replaces a [`Float`] $x$ with $\log_b(1+x)$, where the base $b$ is a [`Float`], rounding the
665 /// result to the nearest value of the input's precision. The base is taken by reference.
666 ///
667 /// See [`Float::log_base_float_base_1_plus_x_prec_round`] for special cases.
668 ///
669 /// # Worst-case complexity
670 /// $T(n, m) = O(n (\log n)^2 \log\log n + m \log m \log\log m)$
671 ///
672 /// $M(n, m) = O(n \log n + m \log m)$
673 ///
674 /// where $T$ is time, $M$ is additional memory, $n$ is the precision of the input, and $m$ is
675 /// `base.significant_bits()`.
676 ///
677 /// # Examples
678 /// ```
679 /// use malachite_base::num::arithmetic::traits::LogBaseOf1PlusXAssign;
680 /// use malachite_float::Float;
681 ///
682 /// let mut x = Float::from(8);
683 /// x.log_base_1_plus_x_assign(&Float::from(3));
684 /// assert_eq!(x.to_string(), "2.0"); // log_3(1 + 8) = log_3(9) = 2
685 /// ```
686 #[inline]
687 fn log_base_1_plus_x_assign(&mut self, base: &Self) {
688 let prec = self.significant_bits();
689 self.log_base_float_base_1_plus_x_prec_round_assign(base, prec, Nearest);
690 }
691}
692
693/// Computes $\log_b(1+x)$, the base-$b$ logarithm of one plus a primitive float, where the base $b$
694/// is also a primitive float, returning a primitive float result. Using this function is more
695/// accurate than computing the logarithm using the standard library, both because $1+x$ may not be
696/// representable as a primitive float and because the standard library's `log` is not always
697/// correctly rounded.
698///
699/// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned. Otherwise the base
700/// may be any primitive float: the function is defined as $\ln(1+x) / \ln b$ and never panics on an
701/// input value. A base in $(0,1)$ gives a (sign-flipped) logarithm, and the non-normal and
702/// degenerate bases follow the limits below.
703///
704/// $$
705/// f(x,b) = \log_b(1+x)+\varepsilon.
706/// $$
707/// - If $\log_b(1+x)$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
708/// - If $\log_b(1+x)$ is finite and nonzero, then $|\varepsilon| < 2^{\lfloor\log_2
709/// |\log_b(1+x)|\rfloor-p}$, where $p$ is precision of the output (typically 24 if `T` is a
710/// [`f32`] and 53 if `T` is a [`f64`], but less if the output is subnormal).
711///
712/// Special cases (with $b$ the base):
713/// - $f(\text{NaN},b)=\text{NaN}$, and $f(x,\text{NaN})=\text{NaN}$
714/// - $f(x,b)=\text{NaN}$ for $x<-1$ or $b<0$
715/// - $f(\infty,b)=\infty$ for $b>1$, and $-\infty$ for $0\leq b<1$
716/// - $f(-1.0,b)=-\infty$ for $b>1$, and $\infty$ for $0<b<1$
717/// - $f(\pm0.0,b)=0$ (the sign of $\pm0.0$ times the sign of $1/\ln b$)
718/// - $f(x,\infty)=0$ for finite $x>-1$ with $x\neq0$ (and $\text{NaN}$ for $x\in\{\infty,-1\}$)
719/// - $f(x,\pm0.0)=0$ for finite $x>-1$ with $x\neq0$ (and $\text{NaN}$ for $x\in\{\infty,-1\}$)
720/// - $f(x,1.0)=\infty$ for $x>0$ or $x=\infty$, $-\infty$ for $-1\leq x<0$, and $\text{NaN}$ for
721/// $x=\pm0.0$
722///
723/// This function can both overflow (for a base near 1) and underflow (for an $x$ near 0).
724///
725/// # Worst-case complexity
726/// Constant time and additional memory.
727///
728/// # Examples
729/// ```
730/// use malachite_base::num::basic::traits::NegativeInfinity;
731/// use malachite_base::num::float::NiceFloat;
732/// use malachite_float::float::arithmetic::log_base_float_base_1_plus_x::*;
733///
734/// // log_4(1 + 3) = log_4(4) = 1
735/// assert_eq!(
736/// NiceFloat(primitive_float_log_base_float_base_1_plus_x(3.0f32, 4.0)),
737/// NiceFloat(1.0)
738/// );
739/// // log_4(1 + 1) = log_4(2) = 1/2
740/// assert_eq!(
741/// NiceFloat(primitive_float_log_base_float_base_1_plus_x(1.0f32, 4.0)),
742/// NiceFloat(0.5)
743/// );
744/// // log_(1/2)(1 + 3) = log_(1/2)(4) = -2
745/// assert_eq!(
746/// NiceFloat(primitive_float_log_base_float_base_1_plus_x(3.0f32, 0.5)),
747/// NiceFloat(-2.0)
748/// );
749/// assert_eq!(
750/// NiceFloat(primitive_float_log_base_float_base_1_plus_x(-1.0f32, 10.0)),
751/// NiceFloat(f32::NEGATIVE_INFINITY)
752/// );
753/// assert!(primitive_float_log_base_float_base_1_plus_x(-2.0f32, 10.0).is_nan());
754/// assert!(primitive_float_log_base_float_base_1_plus_x(3.0f32, f32::NAN).is_nan());
755/// ```
756#[inline]
757#[allow(clippy::type_repetition_in_bounds)]
758pub fn primitive_float_log_base_float_base_1_plus_x<T: PrimitiveFloat>(x: T, base: T) -> T
759where
760 Float: From<T> + PartialOrd<T>,
761 for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
762{
763 emulate_float_float_to_float_fn(
764 |x, base, prec| x.log_base_float_base_1_plus_x_prec(&base, prec),
765 x,
766 base,
767 )
768}