malachite_float/arithmetic/log_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, Zero};
10use crate::{Float, emulate_float_to_float_fn, float_infinity, float_nan, float_negative_infinity};
11use core::cmp::Ordering::{self, *};
12use malachite_base::num::arithmetic::traits::{
13 CeilingLogBase2, CheckedLogBase, IsPowerOf2, LogBaseOf1PlusX, LogBaseOf1PlusXAssign,
14};
15use malachite_base::num::basic::floats::PrimitiveFloat;
16use malachite_base::num::basic::integers::PrimitiveInt;
17use malachite_base::num::basic::traits::{One as OneTrait, Zero as ZeroTrait};
18use malachite_base::num::comparison::traits::PartialOrdAbs;
19use malachite_base::num::conversion::traits::{ExactFrom, RoundingFrom};
20use malachite_base::num::factorization::traits::ExpressAsPower;
21use malachite_base::num::logic::traits::SignificantBits;
22use malachite_base::rounding_modes::RoundingMode::{self, *};
23use malachite_nz::natural::Natural;
24use malachite_nz::natural::arithmetic::float_extras::float_can_round;
25use malachite_nz::platform::Limb;
26use malachite_q::Rational;
27
28// Returns `Some(m / e_base)` -- the value of `log_base(1 + x)` -- when `1 + x = g^m` for the root
29// `g` of `base` (so `base = g^e_base` and `log_base(1 + x)` is rational), and `None` when it is
30// irrational. The input `x` must be finite and greater than -1, and `base > 1` must not be a power
31// of 2.
32//
33// For a non-power-of-2 base, `g` is not a perfect power, so `1 + x = g^m` is an exact (dyadic)
34// `Float` value only when `m >= 0`: `m = 0` gives `x = 0`, and `m >= 1` gives `1 + x` a positive
35// integer power of `g` (so `x` is a positive integer). A negative `m` would make `g^m` a non-dyadic
36// fraction, impossible for the dyadic `1 + x`.
37//
38// Detecting these rational results up front is essential: the Ziv loop in
39// `log_base_1_plus_x_prec_round_normal` could never certify an exactly-representable one. The check
40// is balloon-safe, materializing `x` as an integer only when its exponent is within `64 * prec` of
41// being a representable `g^m - 1`.
42pub(crate) fn log_base_1_plus_x_rational(x: &Float, base: u64) -> Option<Rational> {
43 if *x == 0u32 {
44 return Some(Rational::ZERO);
45 }
46 let e = i64::from(x.get_exponent()?);
47 if e < 1 || u64::exact_from(e) > x.get_prec()?.saturating_mul(64) {
48 return None;
49 }
50 // `Natural::try_from` fails unless `x` is a nonnegative integer.
51 let n = Natural::try_from(x).ok()?;
52 let (g, e_base) = base.express_as_power().unwrap_or((base, 1));
53 let m = (n + Natural::ONE).checked_log_base(&Natural::from(g))?;
54 Some(Rational::from_unsigneds(m, e_base))
55}
56
57// The computation of log_base_1_plus_x(x, base) is done by log_base(1 + x) = log_2(1 + x) /
58// log_2(base). The input is finite and greater than -1, and `base > 1` is not a power of 2.
59//
60// Routing through `log_base_2_1_plus_x` (rather than computing `log_base(1 + x)` from `1 + x`
61// directly) preserves accuracy when x is near 0, where `1 + x` would lose precision. Unlike the
62// power-of-2 case, no near-power-of-2 special handling is needed: `log_2(base)` is irrational, so
63// every non-rational result is strictly between `Float`s and the Ziv loop converges. (The rational
64// results, where `1 + x = g^m`, are detected up front.)
65fn log_base_1_plus_x_prec_round_normal(
66 x: &Float,
67 base: u64,
68 prec: u64,
69 rm: RoundingMode,
70) -> (Float, Ordering) {
71 // log_base(1 + x) is undefined for x < -1.
72 match x.partial_cmp(&-1i32).unwrap() {
73 // 1 + x = 0, so log_base(1 + x) = -infinity (base > 1).
74 Equal => return (float_negative_infinity!(), Equal),
75 Less => return (float_nan!(), Equal),
76 _ => {}
77 }
78 // If 1 + x = g^m, then log_base(1 + x) = m / e_base is rational and exact.
79 if let Some(q) = log_base_1_plus_x_rational(x, base) {
80 return Float::from_rational_prec_round(q, prec, rm);
81 }
82 // The result is irrational, so it is never exactly representable.
83 assert_ne!(rm, Exact, "Inexact log_base_1_plus_x");
84 let base_float = Float::from(base);
85 let min_exp = i64::from(Float::MIN_EXPONENT);
86 let mut working_prec = prec + 4 + prec.ceiling_log_base_2();
87 let mut increment = Limb::WIDTH;
88 loop {
89 // log_2(1 + x), correctly rounded to working_prec; always within the Float exponent range.
90 let num = x.log_base_2_1_plus_x_prec_ref(working_prec).0;
91 // log_2(base) > 1, correctly rounded to working_prec.
92 let den = base_float.log_base_2_prec_ref(working_prec).0;
93 // Dividing by log_2(base) > 1 only shrinks the magnitude (overflow is impossible), but can
94 // push the result below MIN_EXPONENT. When it underflows, the Ziv test below could never
95 // resolve it (the quotient clamps), so hand the rounding to div_prec_round, which clamps to
96 // zero or the minimum positive value per the rounding mode. The exact quotient exponent is
97 // only resolved in the narrow band where the cheap exponent bound is inconclusive (then
98 // e_num - e_den == min_exp - 1, so the result underflows iff |log_2(1 + x)| * 2^(1 -
99 // min_exp) < log_2(base)). The left shift only adjusts the exponent, avoiding a huge
100 // Rational conversion.
101 let e_num = i64::from(num.get_exponent().unwrap());
102 let e_den = i64::from(den.get_exponent().unwrap());
103 if e_num - e_den + 1 < min_exp
104 || (e_num - e_den < min_exp && (&num << u64::exact_from(1 - min_exp)).lt_abs(&den))
105 {
106 return num.div_prec_round(den, prec, rm);
107 }
108 // log_2(1 + x) / log_2(base), with three correctly-rounded operations (log_base_2_1_plus_x,
109 // log_base_2, and the division, each at most 1/2 ulp), so the relative error is below 2^(2
110 // - working_prec) and working_prec - 4 correct bits suffice for rounding.
111 let t = num.div_prec(den, working_prec).0;
112 if float_can_round(t.significand_ref().unwrap(), working_prec - 4, prec, rm) {
113 return Float::from_float_prec_round(t, prec, rm);
114 }
115 // Increase the precision.
116 working_prec += increment;
117 increment = working_prec >> 1;
118 }
119}
120
121impl Float {
122 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
123 /// the result to the specified precision and with the specified rounding mode. The [`Float`] is
124 /// taken by value. An [`Ordering`] is also returned, indicating whether the rounded value is
125 /// less than, equal to, or greater than the exact value. Although `NaN`s are not comparable to
126 /// any [`Float`], whenever this function returns a `NaN` it also returns `Equal`.
127 ///
128 /// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned.
129 ///
130 /// When `base` is a power of 2, this function delegates to
131 /// [`Float::log_base_power_of_2_1_plus_x_prec_round`]; otherwise it computes $\log_2(1+x) /
132 /// \log_2 b$, preserving accuracy for $x$ near 0.
133 ///
134 /// See [`RoundingMode`] for a description of the possible rounding modes.
135 ///
136 /// $$
137 /// f(x,b,p,m) = \log_b(1+x)+\varepsilon.
138 /// $$
139 /// - If $\log_b(1+x)$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to
140 /// be 0.
141 /// - If $\log_b(1+x)$ is finite and nonzero, and $m$ is not `Nearest`, then $|\varepsilon| <
142 /// 2^{\lfloor\log_2 |\log_b(1+x)|\rfloor-p+1}$.
143 /// - If $\log_b(1+x)$ is finite and nonzero, and $m$ is `Nearest`, then $|\varepsilon| \leq
144 /// 2^{\lfloor\log_2 |\log_b(1+x)|\rfloor-p}$.
145 ///
146 /// If the output has a precision, it is `prec`.
147 ///
148 /// Special cases:
149 /// - $f(\text{NaN},b,p,m)=\text{NaN}$
150 /// - $f(\infty,b,p,m)=\infty$
151 /// - $f(-\infty,b,p,m)=\text{NaN}$
152 /// - $f(\pm0.0,b,p,m)=\pm0.0$
153 /// - $f(-1.0,b,p,m)=-\infty$
154 /// - $f(x,b,p,m)=\text{NaN}$ for $x<-1$
155 /// - $f(x,b,p,m)=m/e$ when $1+x=g^m$, where $g$ is the smallest integer of which $b$ is a power
156 /// and $b=g^e$, rounded to precision $p$; the result is exact if and only if $m/e$ is
157 /// representable with precision $p$ (for example $\log_9(1+8)=1$ when $x=8$ is exact)
158 ///
159 /// This function cannot overflow, but it can underflow.
160 ///
161 /// If you know you'll be using `Nearest`, consider using [`Float::log_base_1_plus_x_prec`]
162 /// instead. If you know that your target precision is the precision of the input, consider
163 /// using [`Float::log_base_1_plus_x_round`] instead. If both of these things are true, consider
164 /// using `(&Float).log_base_1_plus_x()` instead.
165 ///
166 /// # Worst-case complexity
167 /// $T(n) = O(n (\log n)^2 \log\log n)$
168 ///
169 /// $M(n) = O(n (\log n)^2)$
170 ///
171 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
172 ///
173 /// # Panics
174 /// Panics if `prec` is zero, if `base` is less than 2, or if `rm` is `Exact` but the result
175 /// cannot be represented exactly with the given precision.
176 ///
177 /// # Examples
178 /// ```
179 /// use malachite_base::rounding_modes::RoundingMode::*;
180 /// use malachite_float::Float;
181 /// use std::cmp::Ordering::*;
182 ///
183 /// let (log, o) = Float::from(8).log_base_1_plus_x_prec_round(9, 10, Exact);
184 /// assert_eq!(log.to_string(), "1.0"); // log_9(9) = 1
185 /// assert_eq!(o, Equal);
186 ///
187 /// let (log, o) = Float::from(1).log_base_1_plus_x_prec_round(3, 20, Nearest);
188 /// assert_eq!(log.to_string(), "0.63093"); // log_3(2)
189 /// assert_eq!(o, Greater);
190 /// ```
191 #[inline]
192 pub fn log_base_1_plus_x_prec_round(
193 self,
194 base: u64,
195 prec: u64,
196 rm: RoundingMode,
197 ) -> (Self, Ordering) {
198 assert_ne!(prec, 0);
199 assert!(base > 1, "Logarithm base must be greater than 1");
200 if base.is_power_of_2() {
201 return self.log_base_power_of_2_1_plus_x_prec_round(
202 i64::from(base.trailing_zeros()),
203 prec,
204 rm,
205 );
206 }
207 match self {
208 Self(NaN | Infinity { sign: false }) => (float_nan!(), Equal),
209 float_infinity!() => (float_infinity!(), Equal),
210 Self(Zero { .. }) => (self, Equal),
211 _ => log_base_1_plus_x_prec_round_normal(&self, base, prec, rm),
212 }
213 }
214
215 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
216 /// the result to the specified precision and with the specified rounding mode. The [`Float`] is
217 /// taken by reference. An [`Ordering`] is also returned, indicating whether the rounded value
218 /// is less than, equal to, or greater than the exact value. Although `NaN`s are not comparable
219 /// to any [`Float`], whenever this function returns a `NaN` it also returns `Equal`.
220 ///
221 /// See [`Float::log_base_1_plus_x_prec_round`] for details, special cases, and a description of
222 /// the rounding behavior.
223 ///
224 /// # Worst-case complexity
225 /// $T(n) = O(n (\log n)^2 \log\log n)$
226 ///
227 /// $M(n) = O(n (\log n)^2)$
228 ///
229 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
230 ///
231 /// # Panics
232 /// Panics if `prec` is zero, if `base` is less than 2, or if `rm` is `Exact` but the result
233 /// cannot be represented exactly with the given precision.
234 ///
235 /// # Examples
236 /// ```
237 /// use malachite_base::rounding_modes::RoundingMode::*;
238 /// use malachite_float::Float;
239 /// use std::cmp::Ordering::*;
240 ///
241 /// let (log, o) = Float::from(8).log_base_1_plus_x_prec_round_ref(3, 10, Exact);
242 /// assert_eq!(log.to_string(), "2.0"); // log_3(9) = 2
243 /// assert_eq!(o, Equal);
244 ///
245 /// let (log, o) = Float::from(1).log_base_1_plus_x_prec_round_ref(3, 20, Floor);
246 /// assert_eq!(log.to_string(), "0.630929"); // log_3(2), rounded down
247 /// assert_eq!(o, Less);
248 /// ```
249 #[inline]
250 pub fn log_base_1_plus_x_prec_round_ref(
251 &self,
252 base: u64,
253 prec: u64,
254 rm: RoundingMode,
255 ) -> (Self, Ordering) {
256 assert_ne!(prec, 0);
257 assert!(base > 1, "Logarithm base must be greater than 1");
258 if base.is_power_of_2() {
259 return self.log_base_power_of_2_1_plus_x_prec_round_ref(
260 i64::from(base.trailing_zeros()),
261 prec,
262 rm,
263 );
264 }
265 match self {
266 Self(NaN | Infinity { sign: false }) => (float_nan!(), Equal),
267 float_infinity!() => (float_infinity!(), Equal),
268 Self(Zero { .. }) => (self.clone(), Equal),
269 _ => log_base_1_plus_x_prec_round_normal(self, base, prec, rm),
270 }
271 }
272
273 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
274 /// the result to the nearest value of the specified precision. The [`Float`] is taken by value.
275 /// An [`Ordering`] is also returned, indicating whether the rounded value is less than, equal
276 /// to, or greater than the exact value.
277 ///
278 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
279 ///
280 /// # Worst-case complexity
281 /// $T(n) = O(n (\log n)^2 \log\log n)$
282 ///
283 /// $M(n) = O(n (\log n)^2)$
284 ///
285 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
286 ///
287 /// # Panics
288 /// Panics if `prec` is zero or if `base` is less than 2.
289 ///
290 /// # Examples
291 /// ```
292 /// use malachite_float::Float;
293 /// use std::cmp::Ordering::*;
294 ///
295 /// let (log, o) = Float::from(8).log_base_1_plus_x_prec(9, 10);
296 /// assert_eq!(log.to_string(), "1.0"); // log_9(9) = 1
297 /// assert_eq!(o, Equal);
298 ///
299 /// let (log, o) = Float::from(1).log_base_1_plus_x_prec(3, 20);
300 /// assert_eq!(log.to_string(), "0.63093"); // log_3(2)
301 /// assert_eq!(o, Greater);
302 /// ```
303 #[inline]
304 pub fn log_base_1_plus_x_prec(self, base: u64, prec: u64) -> (Self, Ordering) {
305 self.log_base_1_plus_x_prec_round(base, prec, Nearest)
306 }
307
308 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
309 /// the result to the nearest value of the specified precision. The [`Float`] is taken by
310 /// reference. An [`Ordering`] is also returned, indicating whether the rounded value is less
311 /// than, equal to, or greater than the exact value.
312 ///
313 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
314 ///
315 /// # Worst-case complexity
316 /// $T(n) = O(n (\log n)^2 \log\log n)$
317 ///
318 /// $M(n) = O(n (\log n)^2)$
319 ///
320 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
321 ///
322 /// # Panics
323 /// Panics if `prec` is zero or if `base` is less than 2.
324 ///
325 /// # Examples
326 /// ```
327 /// use malachite_float::Float;
328 /// use std::cmp::Ordering::*;
329 ///
330 /// let (log, o) = (&Float::from(2)).log_base_1_plus_x_prec_ref(9, 10);
331 /// assert_eq!(log.to_string(), "0.5"); // log_9(3) = 1/2
332 /// assert_eq!(o, Equal);
333 ///
334 /// let (log, o) = (&Float::from(7)).log_base_1_plus_x_prec_ref(5, 30);
335 /// assert_eq!(log.to_string(), "1.292029675"); // log_5(8)
336 /// assert_eq!(o, Greater);
337 /// ```
338 #[inline]
339 pub fn log_base_1_plus_x_prec_ref(&self, base: u64, prec: u64) -> (Self, Ordering) {
340 self.log_base_1_plus_x_prec_round_ref(base, prec, Nearest)
341 }
342
343 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
344 /// the result to the precision of the input and with the specified rounding mode. The [`Float`]
345 /// is taken by value. An [`Ordering`] is also returned, indicating whether the rounded value is
346 /// less than, equal to, or greater than the exact value.
347 ///
348 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
349 ///
350 /// # Worst-case complexity
351 /// $T(n) = O(n (\log n)^2 \log\log n)$
352 ///
353 /// $M(n) = O(n (\log n)^2)$
354 ///
355 /// where $T$ is time, $M$ is additional memory, and $n$ is the precision of the input.
356 ///
357 /// # Panics
358 /// Panics if `base` is less than 2, or if `rm` is `Exact` but the result cannot be represented
359 /// exactly with the input's precision.
360 ///
361 /// # Examples
362 /// ```
363 /// use malachite_base::rounding_modes::RoundingMode::*;
364 /// use malachite_float::Float;
365 /// use std::cmp::Ordering::*;
366 ///
367 /// let (log, o) = Float::from(8).log_base_1_plus_x_round(9, Exact);
368 /// assert_eq!(log.to_string(), "1.0"); // log_9(9) = 1
369 /// assert_eq!(o, Equal);
370 ///
371 /// let (log, o) = Float::from(8).log_base_1_plus_x_round(3, Exact);
372 /// assert_eq!(log.to_string(), "2.0"); // log_3(9) = 2
373 /// assert_eq!(o, Equal);
374 /// ```
375 #[inline]
376 pub fn log_base_1_plus_x_round(self, base: u64, rm: RoundingMode) -> (Self, Ordering) {
377 let prec = self.significant_bits();
378 self.log_base_1_plus_x_prec_round(base, prec, rm)
379 }
380
381 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
382 /// the result to the precision of the input and with the specified rounding mode. The [`Float`]
383 /// is taken by reference. An [`Ordering`] is also returned, indicating whether the rounded
384 /// value is less than, equal to, or greater than the exact value.
385 ///
386 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
387 ///
388 /// # Worst-case complexity
389 /// $T(n) = O(n (\log n)^2 \log\log n)$
390 ///
391 /// $M(n) = O(n (\log n)^2)$
392 ///
393 /// where $T$ is time, $M$ is additional memory, and $n$ is the precision of the input.
394 ///
395 /// # Panics
396 /// Panics if `base` is less than 2, or if `rm` is `Exact` but the result cannot be represented
397 /// exactly with the input's precision.
398 ///
399 /// # Examples
400 /// ```
401 /// use malachite_base::rounding_modes::RoundingMode::*;
402 /// use malachite_float::Float;
403 /// use std::cmp::Ordering::*;
404 ///
405 /// let (log, o) = (&Float::from(8)).log_base_1_plus_x_round_ref(3, Exact);
406 /// assert_eq!(log.to_string(), "2.0"); // log_3(9) = 2
407 /// assert_eq!(o, Equal);
408 ///
409 /// let (log, o) = (&Float::from(2)).log_base_1_plus_x_round_ref(9, Exact);
410 /// assert_eq!(log.to_string(), "0.5"); // log_9(3) = 1/2
411 /// assert_eq!(o, Equal);
412 /// ```
413 #[inline]
414 pub fn log_base_1_plus_x_round_ref(&self, base: u64, rm: RoundingMode) -> (Self, Ordering) {
415 self.log_base_1_plus_x_prec_round_ref(base, self.significant_bits(), rm)
416 }
417
418 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, in
419 /// place, rounding the result to the specified precision and with the specified rounding mode.
420 /// An [`Ordering`] is returned, indicating whether the rounded value is less than, equal to, or
421 /// greater than the exact value.
422 ///
423 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
424 ///
425 /// # Worst-case complexity
426 /// $T(n) = O(n (\log n)^2 \log\log n)$
427 ///
428 /// $M(n) = O(n (\log n)^2)$
429 ///
430 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
431 ///
432 /// # Panics
433 /// Panics if `prec` is zero, if `base` is less than 2, or if `rm` is `Exact` but the result
434 /// cannot be represented exactly with the given precision.
435 ///
436 /// # Examples
437 /// ```
438 /// use malachite_base::rounding_modes::RoundingMode::*;
439 /// use malachite_float::Float;
440 /// use std::cmp::Ordering::*;
441 ///
442 /// let mut x = Float::from(8);
443 /// assert_eq!(x.log_base_1_plus_x_prec_round_assign(9, 10, Exact), Equal);
444 /// assert_eq!(x.to_string(), "1.0"); // log_9(9) = 1
445 ///
446 /// let mut x = Float::from(1);
447 /// assert_eq!(x.log_base_1_plus_x_prec_round_assign(3, 20, Floor), Less);
448 /// assert_eq!(x.to_string(), "0.630929"); // log_3(2), rounded down
449 /// ```
450 #[inline]
451 pub fn log_base_1_plus_x_prec_round_assign(
452 &mut self,
453 base: u64,
454 prec: u64,
455 rm: RoundingMode,
456 ) -> Ordering {
457 let (result, o) = core::mem::take(self).log_base_1_plus_x_prec_round(base, prec, rm);
458 *self = result;
459 o
460 }
461
462 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, in
463 /// place, rounding the result to the nearest value of the specified precision. An [`Ordering`]
464 /// is returned, indicating whether the rounded value is less than, equal to, or greater than
465 /// the exact value.
466 ///
467 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
468 ///
469 /// # Worst-case complexity
470 /// $T(n) = O(n (\log n)^2 \log\log n)$
471 ///
472 /// $M(n) = O(n (\log n)^2)$
473 ///
474 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
475 ///
476 /// # Panics
477 /// Panics if `prec` is zero or if `base` is less than 2.
478 ///
479 /// # Examples
480 /// ```
481 /// use malachite_float::Float;
482 ///
483 /// let mut x = Float::from(8);
484 /// x.log_base_1_plus_x_prec_assign(3, 10);
485 /// assert_eq!(x.to_string(), "2.0"); // log_3(9) = 2
486 ///
487 /// let mut x = Float::from(2);
488 /// x.log_base_1_plus_x_prec_assign(9, 10);
489 /// assert_eq!(x.to_string(), "0.5"); // log_9(3) = 1/2
490 /// ```
491 #[inline]
492 pub fn log_base_1_plus_x_prec_assign(&mut self, base: u64, prec: u64) -> Ordering {
493 self.log_base_1_plus_x_prec_round_assign(base, prec, Nearest)
494 }
495
496 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, in
497 /// place, rounding the result to the precision of the input and with the specified rounding
498 /// mode. An [`Ordering`] is returned, indicating whether the rounded value is less than, equal
499 /// to, or greater than the exact value.
500 ///
501 /// See [`Float::log_base_1_plus_x_prec_round`] for details and special cases.
502 ///
503 /// # Worst-case complexity
504 /// $T(n) = O(n (\log n)^2 \log\log n)$
505 ///
506 /// $M(n) = O(n (\log n)^2)$
507 ///
508 /// where $T$ is time, $M$ is additional memory, and $n$ is the precision of the input.
509 ///
510 /// # Panics
511 /// Panics if `base` is less than 2, or if `rm` is `Exact` but the result cannot be represented
512 /// exactly with the input's precision.
513 ///
514 /// # Examples
515 /// ```
516 /// use malachite_base::rounding_modes::RoundingMode::*;
517 /// use malachite_float::Float;
518 ///
519 /// let mut x = Float::from(8);
520 /// x.log_base_1_plus_x_round_assign(9, Exact);
521 /// assert_eq!(x.to_string(), "1.0"); // log_9(9) = 1
522 ///
523 /// let mut x = Float::from(8);
524 /// x.log_base_1_plus_x_round_assign(3, Exact);
525 /// assert_eq!(x.to_string(), "2.0"); // log_3(9) = 2
526 /// ```
527 #[inline]
528 pub fn log_base_1_plus_x_round_assign(&mut self, base: u64, rm: RoundingMode) -> Ordering {
529 let prec = self.significant_bits();
530 self.log_base_1_plus_x_prec_round_assign(base, prec, rm)
531 }
532}
533
534impl LogBaseOf1PlusX<u64> for Float {
535 type Output = Self;
536
537 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
538 /// the result to the nearest value of the input's precision. The [`Float`] is taken by value.
539 ///
540 /// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned. See
541 /// [`Float::log_base_1_plus_x_prec_round`] for the other special cases.
542 ///
543 /// # Worst-case complexity
544 /// $T(n) = O(n (\log n)^2 \log\log n)$
545 ///
546 /// $M(n) = O(n (\log n)^2)$
547 ///
548 /// where $T$ is time, $M$ is additional memory, and $n$ is the precision of the input.
549 ///
550 /// # Panics
551 /// Panics if `base` is less than 2.
552 ///
553 /// # Examples
554 /// ```
555 /// use malachite_base::num::arithmetic::traits::LogBaseOf1PlusX;
556 /// use malachite_float::Float;
557 ///
558 /// assert_eq!(Float::from(8).log_base_1_plus_x(9).to_string(), "1.0"); // log_9(9) = 1
559 /// assert_eq!(Float::from(8).log_base_1_plus_x(3).to_string(), "2.0"); // log_3(9) = 2
560 /// ```
561 #[inline]
562 fn log_base_1_plus_x(self, base: u64) -> Self {
563 let prec = self.significant_bits();
564 self.log_base_1_plus_x_prec_round(base, prec, Nearest).0
565 }
566}
567
568impl LogBaseOf1PlusX<u64> for &Float {
569 type Output = Float;
570
571 /// Computes $\log_b(1+x)$, where $x$ is a [`Float`] and $b$ is a `u64` greater than 1, rounding
572 /// the result to the nearest value of the input's precision. The [`Float`] is taken by
573 /// reference.
574 ///
575 /// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned. See
576 /// [`Float::log_base_1_plus_x_prec_round`] for the other special cases.
577 ///
578 /// # Worst-case complexity
579 /// $T(n) = O(n (\log n)^2 \log\log n)$
580 ///
581 /// $M(n) = O(n (\log n)^2)$
582 ///
583 /// where $T$ is time, $M$ is additional memory, and $n$ is the precision of the input.
584 ///
585 /// # Panics
586 /// Panics if `base` is less than 2.
587 ///
588 /// # Examples
589 /// ```
590 /// use malachite_base::num::arithmetic::traits::LogBaseOf1PlusX;
591 /// use malachite_float::Float;
592 ///
593 /// assert_eq!((&Float::from(8)).log_base_1_plus_x(9).to_string(), "1.0"); // log_9(9) = 1
594 /// assert_eq!((&Float::from(2)).log_base_1_plus_x(9).to_string(), "0.5"); // log_9(3) = 1/2
595 /// ```
596 #[inline]
597 fn log_base_1_plus_x(self, base: u64) -> Float {
598 self.log_base_1_plus_x_prec_round_ref(base, self.significant_bits(), Nearest)
599 .0
600 }
601}
602
603impl LogBaseOf1PlusXAssign<u64> for Float {
604 /// Replaces a [`Float`] $x$ with $\log_b(1+x)$, where $b$ is a `u64` greater than 1, rounding
605 /// the result to the nearest value of the input's precision.
606 ///
607 /// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned. See
608 /// [`Float::log_base_1_plus_x_prec_round`] for the other special cases.
609 ///
610 /// # Worst-case complexity
611 /// $T(n) = O(n (\log n)^2 \log\log n)$
612 ///
613 /// $M(n) = O(n (\log n)^2)$
614 ///
615 /// where $T$ is time, $M$ is additional memory, and $n$ is the precision of the input.
616 ///
617 /// # Panics
618 /// Panics if `base` is less than 2.
619 ///
620 /// # Examples
621 /// ```
622 /// use malachite_base::num::arithmetic::traits::LogBaseOf1PlusXAssign;
623 /// use malachite_float::Float;
624 ///
625 /// let mut x = Float::from(8);
626 /// x.log_base_1_plus_x_assign(9);
627 /// assert_eq!(x.to_string(), "1.0"); // log_9(9) = 1
628 ///
629 /// let mut x = Float::from(8);
630 /// x.log_base_1_plus_x_assign(3);
631 /// assert_eq!(x.to_string(), "2.0"); // log_3(9) = 2
632 /// ```
633 #[inline]
634 fn log_base_1_plus_x_assign(&mut self, base: u64) {
635 let prec = self.significant_bits();
636 self.log_base_1_plus_x_prec_round_assign(base, prec, Nearest);
637 }
638}
639
640/// Computes $\log_b(1+x)$, the base-$b$ logarithm of one plus a primitive float, where $b$ is a
641/// `u64` greater than 1. Using this function is more accurate than computing the logarithm using
642/// the standard library, both because $1+x$ may not be representable as a primitive float and
643/// because the standard library's `log` is not always correctly rounded.
644///
645/// $\log_b(1+x)$ is undefined for $x<-1$, so whenever $x<-1$, `NaN` is returned.
646///
647/// $$
648/// f(x,b) = \log_b(1+x)+\varepsilon.
649/// $$
650/// - If $\log_b(1+x)$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
651/// - If $\log_b(1+x)$ is finite and nonzero, then $|\varepsilon| < 2^{\lfloor\log_2
652/// |\log_b(1+x)|\rfloor-p}$, where $p$ is precision of the output (typically 24 if `T` is a
653/// [`f32`] and 53 if `T` is a [`f64`], but less if the output is subnormal).
654///
655/// Special cases:
656/// - $f(\text{NaN},b)=\text{NaN}$
657/// - $f(\infty,b)=\infty$
658/// - $f(-\infty,b)=\text{NaN}$
659/// - $f(\pm0.0,b)=\pm0.0$
660/// - $f(-1.0,b)=-\infty$
661/// - $f(x,b)=\text{NaN}$ for $x<-1$
662///
663/// This function can underflow (to a subnormal or zero) when $x$ is close to zero and $b$ is large,
664/// but it cannot overflow.
665///
666/// # Worst-case complexity
667/// Constant time and additional memory.
668///
669/// # Panics
670/// Panics if `base` is less than 2.
671///
672/// # Examples
673/// ```
674/// use malachite_base::num::basic::traits::NegativeInfinity;
675/// use malachite_base::num::float::NiceFloat;
676/// use malachite_float::arithmetic::log_base_1_plus_x::primitive_float_log_base_1_plus_x;
677///
678/// assert!(primitive_float_log_base_1_plus_x(f32::NAN, 10).is_nan());
679/// assert_eq!(
680/// NiceFloat(primitive_float_log_base_1_plus_x(f32::INFINITY, 10)),
681/// NiceFloat(f32::INFINITY)
682/// );
683/// assert_eq!(
684/// NiceFloat(primitive_float_log_base_1_plus_x(-1.0f32, 10)),
685/// NiceFloat(f32::NEGATIVE_INFINITY)
686/// );
687/// assert!(primitive_float_log_base_1_plus_x(-2.0f32, 10).is_nan());
688/// // log_10(1 + 999) = log_10(1000) = 3
689/// assert_eq!(
690/// NiceFloat(primitive_float_log_base_1_plus_x(999.0f32, 10)),
691/// NiceFloat(3.0)
692/// );
693/// // log_9(1 + 8) = log_9(9) = 1
694/// assert_eq!(
695/// NiceFloat(primitive_float_log_base_1_plus_x(8.0f32, 9)),
696/// NiceFloat(1.0)
697/// );
698/// // log_3(1 + 1) = log_3(2)
699/// assert_eq!(
700/// NiceFloat(primitive_float_log_base_1_plus_x(1.0f32, 3)),
701/// NiceFloat(0.63092977)
702/// );
703/// ```
704#[inline]
705#[allow(clippy::type_repetition_in_bounds)]
706pub fn primitive_float_log_base_1_plus_x<T: PrimitiveFloat>(x: T, base: u64) -> T
707where
708 Float: From<T> + PartialOrd<T>,
709 for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
710{
711 emulate_float_to_float_fn(|x, prec| Float::log_base_1_plus_x_prec(x, base, prec), x)
712}