malachite_float/float/arithmetic/compound.rs
1// Copyright © 2026 Mikhail Hogrefe
2//
3// Uses code adopted from the GNU MPFR Library.
4//
5// Copyright 2021-2025 Free Software Foundation, Inc.
6//
7// Contributed by the Pascaline and Caramba projects, INRIA.
8//
9// This file is part of Malachite.
10//
11// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
12// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
13// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
14
15use crate::Float;
16use crate::InnerFloat::{Infinity, NaN, Zero};
17use crate::emulate_float_to_float_fn;
18use crate::float::arithmetic::exp::{exp_overflow, exp_underflow};
19use core::cmp::Ordering::{self, *};
20use malachite_base::num::arithmetic::traits::{CeilingLogBase2, Compound, CompoundAssign, Sign};
21use malachite_base::num::basic::floats::PrimitiveFloat;
22use malachite_base::num::basic::integers::PrimitiveInt;
23use malachite_base::num::basic::traits::{
24 Infinity as InfinityTrait, NaN as NaNTrait, One, Zero as ZeroTrait,
25};
26use malachite_base::num::conversion::traits::{ExactFrom, RoundingFrom};
27use malachite_base::num::logic::traits::SignificantBits;
28use malachite_base::rounding_modes::RoundingMode::{self, *};
29use malachite_nz::natural::arithmetic::float::round::float_can_round;
30use malachite_nz::platform::Limb;
31
32// Rounds (1+x)^n to `prec` bits, assuming |(1+x)^n - 1| < (1/4)ulp(1) = 2^(-prec-2), where `s_pos`
33// is the sign of n*log2(1+x) (true if positive; that quantity is nonzero here).
34//
35// This is mpfr_compound_near_one from compound.c, MPFR 4.2.2.
36fn compound_near_one(prec: u64, s_pos: bool, rm: RoundingMode) -> (Float, Ordering) {
37 let mut y = Float::one_prec(prec);
38 match rm {
39 Exact => panic!("compound: Exact rounding was requested, but the result is inexact"),
40 // round toward 1
41 Nearest => (y, if s_pos { Less } else { Greater }),
42 Down | Floor if s_pos => (y, Less),
43 Up | Ceiling if !s_pos => (y, Greater),
44 // round toward +Inf
45 Up | Ceiling => {
46 y.increment();
47 (y, Greater)
48 }
49 // necessarily Down or Floor with a negative sign; round toward 0
50 _ => {
51 y.decrement();
52 (y, Less)
53 }
54 }
55}
56
57// A shortcut for cases where Ziv's strategy may take too much memory and be too long, i.e. when x^n
58// fits in the target precision (+ 1 additional bit for rounding to nearest) and the exact result
59// (1+x)^n is very close to x^n. Necessarily, x is a large even integer and n > 1. The kx < ex test
60// checks that x is an even integer (iff its least bit 1 has exponent >= 1), and the test after it
61// is a simple condition that implies that x^n fits in the target precision. Here are the details:
62// let k be the minimum length of the significand of x, and x' the odd (integer) significand of x.
63// This means that 2^(k-1) <= x' < 2^k. Thus 2^(n*(k-1)) <= (x')^n < 2^(k*n), and x^n has between
64// n*(k-1)+1 and k*n bits. So x^n can fit into p bits only if p >= n*(k-1)+1, i.e. n*(k-1) <= p-1.
65//
66// This is the "check if x^n fits" portion of mpfr_compound_si from compound.c, MPFR 4.2.2.
67fn compound_x_n_fits(
68 x: &Float,
69 n: i64,
70 prec: u64,
71 rm: RoundingMode,
72 wprec: u64,
73) -> Option<(Float, Ordering)> {
74 let ex = i64::from(x.get_exponent().unwrap());
75 if ex < 17 {
76 return None;
77 }
78 let kx = x.get_min_prec().unwrap();
79 let p = prec + u64::from(rm == Nearest);
80 if kx >= u64::exact_from(ex)
81 || u128::from(n.unsigned_abs()) * u128::from(kx - 1) > u128::from(p - 1)
82 {
83 return None;
84 }
85 // Check whether x^n really fits into p bits.
86 let (v, o_v) = x.pow_u_prec_round_ref(u64::exact_from(n), p, Down);
87 if o_v != Equal {
88 return None;
89 }
90 // (x+1)^n = x^n * (1 + 1/x)^n For directed rounding, we can round when (1 + 1/x)^n < 1 + 2^-p,
91 // and then the result is x^n, except for rounding up. Indeed, if (1 + 1/x)^n < 1 + 2^-p, 1 <=
92 // (x+1)^n < x^n * (1 + 2^-p) = x^n + x^n/2^p < x^n + ulp(x^n). For rounding to nearest, we can
93 // round when (1 + 1/x)^n < 1 + 2^-p, and then the result is x^n when x^n fits into p-1 bits,
94 // and nextabove(x^n) otherwise.
95 let mut r = x.reciprocal_prec_round_ref(wprec, Up).0;
96 r.add_prec_round_assign(Float::ONE, wprec, Up);
97 r.pow_u_round_assign(u64::exact_from(n), Up);
98 r.sub_prec_round_assign(Float::ONE, wprec, Up);
99 // r cannot be zero
100 if i64::from(r.get_exponent().unwrap()) >= -i64::exact_from(prec) {
101 return None;
102 }
103 let v_min_prec = v.get_min_prec().unwrap();
104 let mut y = Float::from_float_prec_round(v, prec, Down).0;
105 Some(
106 if (rm == Nearest && v_min_prec == p) || rm == Up || rm == Ceiling {
107 // round up
108 y.increment();
109 (y, Greater)
110 } else {
111 (y, Less)
112 },
113 )
114}
115
116// This is mpfr_compound_si from compound.c, MPFR 4.2.2, with two corrections taken from the MPFR
117// development sources: log2p1 is rounded toward zero unconditionally (4.2.2 chooses the direction
118// from the signs of x and n, which is backwards for negative n and can yield a result off by one
119// ulp in the min_prec escape below -- confirmed against 4.2.2 via rug and against exact rational
120// arithmetic), and the rounding tests are skipped when e >= precu (when the error bound on u is too
121// large to say anything). MPFR also runs the computation in its extended exponent range and maps
122// back at the end via mpfr_check_range; we instead cut overflow and underflow against the real
123// exponent range up front. This is safe because u is rounded toward zero (making the cuts sound),
124// and because 2^u is rounded toward 1, which keeps the intermediate t representable whenever u
125// survives the cuts.
126fn compound_prec_round_helper(x: &Float, n: i64, prec: u64, rm: RoundingMode) -> (Float, Ordering) {
127 assert_ne!(prec, 0);
128 // Special cases
129 match x {
130 // compound(-Inf, n) is NaN, even for n == 0
131 Float(Infinity { sign: false }) => return (Float::NAN, Equal),
132 // compound(NaN, 0) = 1, like compound(x, 0) for any x >= -1; otherwise NaN propagates
133 Float(NaN) => {
134 return if n == 0 {
135 (Float::one_prec(prec), Equal)
136 } else {
137 (Float::NAN, Equal)
138 };
139 }
140 // compound(0, n) = 1
141 Float(Zero { .. }) => return (Float::one_prec(prec), Equal),
142 // compound(+Inf, 0) = 1, and otherwise (1 + Inf)^n is +0 for n < 0 and +Inf for n > 0
143 Float(Infinity { .. }) => {
144 return match n.sign() {
145 Equal => (Float::one_prec(prec), Equal),
146 Less => (Float::ZERO, Equal),
147 Greater => (Float::INFINITY, Equal),
148 };
149 }
150 _ => {}
151 }
152 // (1+x)^n = NaN for x < -1
153 let compared = x.partial_cmp(&-1i32).unwrap();
154 if compared == Less {
155 return (Float::NAN, Equal);
156 }
157 // compound(x, 0) gives 1 for x >= -1
158 if n == 0 {
159 return (Float::one_prec(prec), Equal);
160 }
161 if compared == Equal {
162 return if n < 0 {
163 // compound(-1, n) = +Inf (MPFR also raises the divide-by-zero exception)
164 (Float::INFINITY, Equal)
165 } else {
166 // compound(-1, n) = +0
167 (Float::ZERO, Equal)
168 };
169 }
170 if n == 1 {
171 return x.add_prec_round_ref_val(Float::ONE, prec, rm);
172 }
173 let mut wprec = prec + prec.ceiling_log_base_2() + 6;
174 // |n| <= 2^k
175 let k = i64::exact_from(n.unsigned_abs().ceiling_log_base_2());
176 let nf = Float::from(n);
177 // We compute u = log2p1(x) with wprec + extra bits, since we lose some bits in 2^u.
178 let mut extra = 0u64;
179 let mut increment = Limb::WIDTH;
180 let mut nloop = 0u32;
181 let t = loop {
182 let precu = wprec + extra;
183 // We compute (1+x)^n as 2^(n*log2p1(x)), and we round toward 1, thus we round n*log2p1(x)
184 // toward 0, which implies we round log2p1(x) toward 0. lg is nonzero and cannot underflow:
185 // |log2(1+x)| > |x| >= 2^(MIN_EXPONENT-1), and toward-zero rounding cannot take it below
186 // the minimum positive Float.
187 let (lg, o_lg) = x.log_base_2_1_plus_x_prec_round_ref(precu, Down);
188 let mut inex = o_lg != Equal;
189 let mut e = i64::from(lg.get_exponent().unwrap());
190 // |lg - log2(1+x)| <= ulp(lg) = 2^(e-precu)
191 let (u, o_mul) = lg.mul_prec_round_val_ref(&nf, precu, Down);
192 inex |= o_mul != Equal;
193 // u is nonzero: |lg| >= 2^(MIN_EXPONENT-1) and |n| >= 1, and the toward-zero rounding of
194 // the product cannot reach below the minimum positive Float.
195 let e2 = i64::from(u.get_exponent().unwrap());
196 // ```
197 // |u - n*log2(1+x)| <= 2^(e2-precu) + |n|*2^(e-precu)
198 // <= 2^(e2-precu) + 2^(e+k-precu) <= 2^(e+k+1-precu)
199 // ``` where |n| <= 2^k, and e2 is the new exponent of u.
200 debug_assert!(e2 <= e + k);
201 e += k + 1;
202 let new_extra = if e2 > 0 { u64::exact_from(e2) } else { 0 };
203 // |u - n*log2(1+x)| <= 2^(e-precu) detect overflow: since we rounded n*log2p1(x) toward 0,
204 // if n*log2p1(x) >= MAX_EXPONENT, we are sure there is overflow.
205 if u >= Float::MAX_EXPONENT {
206 return exp_overflow(prec, rm);
207 }
208 // detect underflow: similarly, since we rounded n*log2p1(x) toward 0, if n*log2p1(x) <
209 // MIN_EXPONENT - 1, we are sure there is underflow.
210 if u < const { Float::MIN_EXPONENT - 1 } {
211 return exp_underflow(prec, if rm == Nearest { Down } else { rm });
212 }
213 // Detect cases where the result is 1 or 1+ulp(1) or 1-(1/2)ulp(1): |2^u - 1| =
214 // |exp(u*log(2)) - 1| <= |u|*log(2) < |u|
215 if nloop == 0 && e2 < -i64::exact_from(prec) {
216 // since ulp(1) = 2^(1-prec), we have |u| < (1/4)ulp(1)
217 return compound_near_one(prec, u.is_sign_positive(), rm);
218 }
219 // round 2^u toward 1
220 let rnd2 = if u.is_sign_positive() { Floor } else { Ceiling };
221 let (mut t, o_exp2) = Float::power_of_2_of_float_prec_round(u, wprec, rnd2);
222 inex |= o_exp2 != Equal;
223 // we had |u - n*log2(1+x)| < 2^(e-precu), thus u = n*log2(1+x) + delta with |delta| <
224 // 2^(e-precu), then 2^u = (1+x)^n * 2^delta. For |delta| < 0.5, |2^delta - 1| <= |delta|
225 // thus |t - (1+x)^n| <= ulp(t) + |t|*2^(e-precu) < 2^(EXP(t)-wprec) + 2^(EXP(t)+e-precu) If
226 // e >= precu, the rounding error on u is too large, and we have to loop again (though the
227 // escapes below may still exit the loop).
228 if e < i64::exact_from(precu) {
229 let extra_i = i64::exact_from(precu - wprec);
230 let err = if extra_i >= e { 1 } else { e + 1 - extra_i };
231 // now |t - (1+x)^n| < 2^(EXP(t)+err-wprec)
232 if !inex
233 || (rm != Exact
234 && i64::exact_from(wprec) > err
235 && float_can_round(
236 t.significand_ref().unwrap(),
237 wprec - u64::exact_from(err),
238 prec,
239 rm,
240 ))
241 {
242 break t;
243 }
244 // If t fits in the target precision (or with 1 more bit), then we can round, assuming
245 // the working precision is large enough, but the above float_can_round will fail
246 // because we cannot determine the ternary value. However, since we rounded t toward 1,
247 // we can determine it. Since the error in the approximation t is at most 2^err ulp(t),
248 // this error should be less than (1/2)ulp(y), thus we should have wprec - prec >= err +
249 // 1. (For Exact rounding we skip this escape, since nudging t would turn an
250 // exactly-representable result into a spurious panic; the exact-1+x escape below
251 // decides exactness instead.)
252 if rm != Exact
253 && t.get_min_prec().unwrap() <= prec + 1
254 && i64::exact_from(wprec - prec) > err
255 {
256 // we step t one place away from 1 to get the correct rounding
257 if rnd2 == Floor {
258 // t was rounded downwards. t cannot be the largest finite significand (its
259 // min_prec is at most prec + 1 < wprec), so this cannot overflow.
260 t.increment();
261 break t;
262 }
263 if t.get_min_prec() != Some(1) || t.get_exponent() != Some(Float::MIN_EXPONENT) {
264 t.decrement();
265 break t;
266 }
267 // Otherwise t is the minimum positive Float, and stepping below it would leave the
268 // representable exponent range. (In MPFR's extended exponent range the step and the
269 // final rounding happen normally, and mpfr_check_range then maps the result back;
270 // the following resolution is equivalent.) The true result lies strictly below t --
271 // t was rounded toward 1 and inex holds, so some rounding was strict -- but within
272 // half an ulp of the target precision, so the rounding resolves directly.
273 return match rm {
274 Floor | Down => (Float::ZERO, Less),
275 // Ceiling, Up, or Nearest; rm is not Exact here
276 _ => (Float::min_positive_value_prec(prec), Greater),
277 };
278 }
279 }
280 // Detect particular cases where Ziv's strategy may take too much memory and be too long.
281 // Since this does not depend on the working precision, we only check this at the first
282 // iteration.
283 debug_assert!(!(0..=1).contains(&n));
284 if nloop == 0
285 && n > 1
286 && let Some(result) = compound_x_n_fits(x, n, prec, rm, wprec)
287 {
288 return result;
289 }
290 // Exact cases like compound(0.5, 2) = 9/4 must be detected, since except for 1+x a power of
291 // 2, the log2p1 above will be inexact, so that in the Ziv test, inex != 0 and
292 // float_can_round will fail (even for Nearest, as the ternary value cannot be determined),
293 // yielding an infinite loop. For an exact case in precision prec, 1+x will necessarily be
294 // exact in precision prec, thus also in wprec, where wprec >= prec, and we can use pow_s
295 // under this condition (which will also evaluate some non-exact cases).
296 let (s, o_s) = x.add_prec_round_ref_val(Float::ONE, wprec, Down);
297 if o_s == Equal {
298 return s.pow_s_prec_round(n, prec, rm);
299 }
300 wprec += increment;
301 increment = wprec >> 1;
302 extra = new_extra;
303 nloop += 1;
304 };
305 Float::from_float_prec_round(t, prec, rm)
306}
307
308impl Float {
309 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
310 /// result to the specified precision and with the specified rounding mode. The [`Float`] is
311 /// taken by value. An [`Ordering`] is also returned, indicating whether the rounded value is
312 /// less than, equal to, or greater than the exact value. Although `NaN`s are not comparable to
313 /// any [`Float`], whenever this function returns a `NaN` it also returns `Equal`.
314 ///
315 /// The compound function is defined in IEEE 754 and is useful for computing compound interest:
316 /// if $x$ is an interest rate, then $(1+x)^n$ is the factor by which a principal grows after
317 /// $n$ compounding periods.
318 ///
319 /// See [`RoundingMode`] for a description of the possible rounding modes.
320 ///
321 /// $$
322 /// f(x,n,p,m) = (1+x)^n+\varepsilon.
323 /// $$
324 /// - If $(1+x)^n$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
325 /// - If $(1+x)^n$ is finite and nonzero, and $m$ is not `Nearest`, then $|\varepsilon| <
326 /// 2^{\lfloor\log_2 (1+x)^n\rfloor-p+1}$.
327 /// - If $(1+x)^n$ is finite and nonzero, and $m$ is `Nearest`, then $|\varepsilon| \leq
328 /// 2^{\lfloor\log_2 (1+x)^n\rfloor-p}$.
329 ///
330 /// Special cases:
331 /// - $f(\text{NaN},n)=\text{NaN}$ if $n\neq 0$, and $1.0$ if $n=0$
332 /// - $f(-\infty,n)=\text{NaN}$, even if $n=0$
333 /// - $f(\infty,0)=1.0$
334 /// - $f(\infty,n)=\infty$ if $n>0$, and $0.0$ if $n<0$
335 /// - $f(\pm 0.0,n)=1.0$
336 /// - $f(x,n)=\text{NaN}$ if $x<-1$, even if $n=0$
337 /// - $f(-1.0,n)=1.0$ if $n=0$, $0.0$ if $n>0$, and $\infty$ if $n<0$
338 /// - $f(x,0)=1.0$ if $x\geq -1$
339 ///
340 /// The result is never negative, and a zero result is always positive.
341 ///
342 /// Overflow and underflow:
343 /// - If $f(x,n,p,m)\geq 2^{2^{30}-1}$ and $m$ is `Ceiling`, `Up`, or `Nearest`, $\infty$ is
344 /// returned instead.
345 /// - If $f(x,n,p,m)\geq 2^{2^{30}-1}$ and $m$ is `Floor` or `Down`, $(1-(1/2)^p)2^{2^{30}-1}$
346 /// is returned instead.
347 /// - If $0<f(x,n,p,m)<2^{-2^{30}}$ and $m$ is `Floor` or `Down`, $0.0$ is returned instead.
348 /// - If $0<f(x,n,p,m)<2^{-2^{30}}$ and $m$ is `Ceiling` or `Up`, $2^{-2^{30}}$ is returned
349 /// instead.
350 /// - If $0<f(x,n,p,m)\leq 2^{-2^{30}-1}$ and $m$ is `Nearest`, $0.0$ is returned instead.
351 /// - If $2^{-2^{30}-1}<f(x,n,p,m)<2^{-2^{30}}$ and $m$ is `Nearest`, either $0.0$ or
352 /// $2^{-2^{30}}$ may be returned. This matches the behavior of MPFR's compound function,
353 /// whose underflow test rounds to nearest as if it were rounding toward zero, except for
354 /// inputs that it resolves by exact powering.
355 ///
356 /// If you know you'll be using `Nearest`, consider using [`Float::compound_prec`] instead. If
357 /// you know that your target precision is the precision of the input, consider using
358 /// [`Float::compound_round`] instead. If both of these things are true, consider using the
359 /// [`Compound`] trait instead.
360 ///
361 /// # Worst-case complexity
362 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
363 ///
364 /// $M(n) = O(n \log n)$
365 ///
366 /// where $T$ is time, $M$ is additional memory, $n$ is `max(prec, self.significant_bits())`,
367 /// and $m$ is the number of significant bits of the exponent `n`.
368 ///
369 /// # Panics
370 /// Panics if `prec` is zero, or if `rm` is `Exact` but the result cannot be represented exactly
371 /// with the given precision.
372 ///
373 /// # Examples
374 /// ```
375 /// use malachite_base::num::basic::traits::Two;
376 /// use malachite_base::rounding_modes::RoundingMode::*;
377 /// use malachite_float::Float;
378 /// use std::cmp::Ordering::*;
379 ///
380 /// let (c, o) = Float::from(3).compound_prec_round(2, 10, Nearest);
381 /// assert_eq!(c.to_string(), "16.000");
382 /// assert_eq!(o, Equal);
383 ///
384 /// let (c, o) = Float::TWO.compound_prec_round(-2, 10, Floor);
385 /// assert_eq!(c.to_string(), "0.11108");
386 /// assert_eq!(o, Less);
387 ///
388 /// let (c, o) = Float::TWO.compound_prec_round(-2, 10, Ceiling);
389 /// assert_eq!(c.to_string(), "0.11121");
390 /// assert_eq!(o, Greater);
391 /// ```
392 #[inline]
393 pub fn compound_prec_round(self, n: i64, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
394 compound_prec_round_helper(&self, n, prec, rm)
395 }
396
397 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
398 /// result to the specified precision and with the specified rounding mode. The [`Float`] is
399 /// taken by reference. An [`Ordering`] is also returned, indicating whether the rounded value
400 /// is less than, equal to, or greater than the exact value. Although `NaN`s are not comparable
401 /// to any [`Float`], whenever this function returns a `NaN` it also returns `Equal`.
402 ///
403 /// See [`RoundingMode`] for a description of the possible rounding modes.
404 ///
405 /// $$
406 /// f(x,n,p,m) = (1+x)^n+\varepsilon.
407 /// $$
408 /// - If $(1+x)^n$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
409 /// - If $(1+x)^n$ is finite and nonzero, and $m$ is not `Nearest`, then $|\varepsilon| <
410 /// 2^{\lfloor\log_2 (1+x)^n\rfloor-p+1}$.
411 /// - If $(1+x)^n$ is finite and nonzero, and $m$ is `Nearest`, then $|\varepsilon| \leq
412 /// 2^{\lfloor\log_2 (1+x)^n\rfloor-p}$.
413 ///
414 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
415 /// overflow, and underflow.
416 ///
417 /// If you know you'll be using `Nearest`, consider using [`Float::compound_prec_ref`] instead.
418 /// If you know that your target precision is the precision of the input, consider using
419 /// [`Float::compound_round_ref`] instead. If both of these things are true, consider using the
420 /// [`Compound`] trait instead.
421 ///
422 /// # Worst-case complexity
423 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
424 ///
425 /// $M(n) = O(n \log n)$
426 ///
427 /// where $T$ is time, $M$ is additional memory, $n$ is `max(prec, self.significant_bits())`,
428 /// and $m$ is the number of significant bits of the exponent `n`.
429 ///
430 /// # Panics
431 /// Panics if `prec` is zero, or if `rm` is `Exact` but the result cannot be represented exactly
432 /// with the given precision.
433 ///
434 /// # Examples
435 /// ```
436 /// use malachite_base::num::basic::traits::Two;
437 /// use malachite_base::rounding_modes::RoundingMode::*;
438 /// use malachite_float::Float;
439 /// use std::cmp::Ordering::*;
440 ///
441 /// let (c, o) = Float::from(3).compound_prec_round_ref(2, 10, Nearest);
442 /// assert_eq!(c.to_string(), "16.000");
443 /// assert_eq!(o, Equal);
444 ///
445 /// let (c, o) = Float::TWO.compound_prec_round_ref(-2, 10, Ceiling);
446 /// assert_eq!(c.to_string(), "0.11121");
447 /// assert_eq!(o, Greater);
448 /// ```
449 #[inline]
450 pub fn compound_prec_round_ref(&self, n: i64, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
451 compound_prec_round_helper(self, n, prec, rm)
452 }
453
454 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
455 /// result to the nearest value of the specified precision. The [`Float`] is taken by value. An
456 /// [`Ordering`] is also returned, indicating whether the rounded value is less than, equal to,
457 /// or greater than the exact value. Although `NaN`s are not comparable to any [`Float`],
458 /// whenever this function returns a `NaN` it also returns `Equal`.
459 ///
460 /// If the compound value is equidistant from two [`Float`]s with the specified precision, the
461 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
462 /// description of the `Nearest` rounding mode.
463 ///
464 /// $$
465 /// f(x,n,p) = (1+x)^n+\varepsilon.
466 /// $$
467 /// - If $(1+x)^n$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
468 /// - If $(1+x)^n$ is finite and nonzero, then $|\varepsilon| \leq 2^{\lfloor\log_2
469 /// (1+x)^n\rfloor-p}$.
470 ///
471 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
472 /// overflow, and underflow.
473 ///
474 /// If you want to use a rounding mode other than `Nearest`, consider using
475 /// [`Float::compound_prec_round`] instead. If you know that your target precision is the
476 /// precision of the input, consider using the [`Compound`] trait instead.
477 ///
478 /// # Worst-case complexity
479 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
480 ///
481 /// $M(n) = O(n \log n)$
482 ///
483 /// where $T$ is time, $M$ is additional memory, $n$ is `max(prec, self.significant_bits())`,
484 /// and $m$ is the number of significant bits of the exponent `n`.
485 ///
486 /// # Panics
487 /// Panics if `prec` is zero.
488 ///
489 /// # Examples
490 /// ```
491 /// use malachite_base::num::basic::traits::Two;
492 /// use malachite_float::Float;
493 /// use std::cmp::Ordering::*;
494 ///
495 /// let (c, o) = Float::from(3).compound_prec(2, 10);
496 /// assert_eq!(c.to_string(), "16.000");
497 /// assert_eq!(o, Equal);
498 ///
499 /// let (c, o) = Float::TWO.compound_prec(-2, 10);
500 /// assert_eq!(c.to_string(), "0.11108");
501 /// assert_eq!(o, Less);
502 /// ```
503 #[inline]
504 pub fn compound_prec(self, n: i64, prec: u64) -> (Self, Ordering) {
505 self.compound_prec_round(n, prec, Nearest)
506 }
507
508 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
509 /// result to the nearest value of the specified precision. The [`Float`] is taken by reference.
510 /// An [`Ordering`] is also returned, indicating whether the rounded value is less than, equal
511 /// to, or greater than the exact value. Although `NaN`s are not comparable to any [`Float`],
512 /// whenever this function returns a `NaN` it also returns `Equal`.
513 ///
514 /// If the compound value is equidistant from two [`Float`]s with the specified precision, the
515 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
516 /// description of the `Nearest` rounding mode.
517 ///
518 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
519 /// overflow, and underflow.
520 ///
521 /// If you want to use a rounding mode other than `Nearest`, consider using
522 /// [`Float::compound_prec_round_ref`] instead. If you know that your target precision is the
523 /// precision of the input, consider using the [`Compound`] trait instead.
524 ///
525 /// # Worst-case complexity
526 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
527 ///
528 /// $M(n) = O(n \log n)$
529 ///
530 /// where $T$ is time, $M$ is additional memory, $n$ is `max(prec, self.significant_bits())`,
531 /// and $m$ is the number of significant bits of the exponent `n`.
532 ///
533 /// # Panics
534 /// Panics if `prec` is zero.
535 ///
536 /// # Examples
537 /// ```
538 /// use malachite_base::num::basic::traits::Two;
539 /// use malachite_float::Float;
540 /// use std::cmp::Ordering::*;
541 ///
542 /// let (c, o) = Float::from(3).compound_prec_ref(2, 10);
543 /// assert_eq!(c.to_string(), "16.000");
544 /// assert_eq!(o, Equal);
545 ///
546 /// let (c, o) = Float::TWO.compound_prec_ref(-2, 10);
547 /// assert_eq!(c.to_string(), "0.11108");
548 /// assert_eq!(o, Less);
549 /// ```
550 #[inline]
551 pub fn compound_prec_ref(&self, n: i64, prec: u64) -> (Self, Ordering) {
552 self.compound_prec_round_ref(n, prec, Nearest)
553 }
554
555 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
556 /// result to the precision of the input with the specified rounding mode. The [`Float`] is
557 /// taken by value. An [`Ordering`] is also returned, indicating whether the rounded value is
558 /// less than, equal to, or greater than the exact value. Although `NaN`s are not comparable to
559 /// any [`Float`], whenever this function returns a `NaN` it also returns `Equal`.
560 ///
561 /// See [`RoundingMode`] for a description of the possible rounding modes.
562 ///
563 /// $$
564 /// f(x,n,m) = (1+x)^n+\varepsilon.
565 /// $$
566 /// - If $(1+x)^n$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
567 /// - If $(1+x)^n$ is finite and nonzero, and $m$ is not `Nearest`, then $|\varepsilon| <
568 /// 2^{\lfloor\log_2 (1+x)^n\rfloor-p+1}$, where $p$ is the precision of the input. Similarly,
569 /// $p$ is the precision of the input in the `Nearest` bullet below.
570 /// - If $(1+x)^n$ is finite and nonzero, and $m$ is `Nearest`, then $|\varepsilon| \leq
571 /// 2^{\lfloor\log_2 (1+x)^n\rfloor-p}$.
572 ///
573 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
574 /// overflow, and underflow.
575 ///
576 /// If you know you'll be using `Nearest`, consider using the [`Compound`] trait instead. If you
577 /// want to specify an output precision, consider using [`Float::compound_prec_round`] instead.
578 ///
579 /// # Worst-case complexity
580 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
581 ///
582 /// $M(n) = O(n \log n)$
583 ///
584 /// where $T$ is time, $M$ is additional memory, $n$ is `self.significant_bits()`, and $m$ is
585 /// the number of significant bits of the exponent `n`.
586 ///
587 /// # Panics
588 /// Panics if `rm` is `Exact` but the result cannot be represented exactly with the precision of
589 /// the input.
590 ///
591 /// # Examples
592 /// ```
593 /// use malachite_base::rounding_modes::RoundingMode::*;
594 /// use malachite_float::Float;
595 /// use std::cmp::Ordering::*;
596 ///
597 /// let (c, o) = Float::from(1.5).compound_round(2, Floor);
598 /// assert_eq!(c.to_string(), "6.0");
599 /// assert_eq!(o, Less);
600 ///
601 /// let (c, o) = Float::from(1.5).compound_round(2, Ceiling);
602 /// assert_eq!(c.to_string(), "8.0");
603 /// assert_eq!(o, Greater);
604 /// ```
605 #[inline]
606 pub fn compound_round(self, n: i64, rm: RoundingMode) -> (Self, Ordering) {
607 let prec = self.significant_bits();
608 self.compound_prec_round(n, prec, rm)
609 }
610
611 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
612 /// result to the precision of the input with the specified rounding mode. The [`Float`] is
613 /// taken by reference. An [`Ordering`] is also returned, indicating whether the rounded value
614 /// is less than, equal to, or greater than the exact value. Although `NaN`s are not comparable
615 /// to any [`Float`], whenever this function returns a `NaN` it also returns `Equal`.
616 ///
617 /// See [`RoundingMode`] for a description of the possible rounding modes.
618 ///
619 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
620 /// overflow, and underflow.
621 ///
622 /// If you know you'll be using `Nearest`, consider using the [`Compound`] trait instead. If you
623 /// want to specify an output precision, consider using [`Float::compound_prec_round_ref`]
624 /// instead.
625 ///
626 /// # Worst-case complexity
627 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
628 ///
629 /// $M(n) = O(n \log n)$
630 ///
631 /// where $T$ is time, $M$ is additional memory, $n$ is `self.significant_bits()`, and $m$ is
632 /// the number of significant bits of the exponent `n`.
633 ///
634 /// # Panics
635 /// Panics if `rm` is `Exact` but the result cannot be represented exactly with the precision of
636 /// the input.
637 ///
638 /// # Examples
639 /// ```
640 /// use malachite_base::rounding_modes::RoundingMode::*;
641 /// use malachite_float::Float;
642 /// use std::cmp::Ordering::*;
643 ///
644 /// let (c, o) = Float::from(1.5).compound_round_ref(2, Floor);
645 /// assert_eq!(c.to_string(), "6.0");
646 /// assert_eq!(o, Less);
647 ///
648 /// let (c, o) = Float::from(1.5).compound_round_ref(2, Ceiling);
649 /// assert_eq!(c.to_string(), "8.0");
650 /// assert_eq!(o, Greater);
651 /// ```
652 #[inline]
653 pub fn compound_round_ref(&self, n: i64, rm: RoundingMode) -> (Self, Ordering) {
654 let prec = self.significant_bits();
655 self.compound_prec_round_ref(n, prec, rm)
656 }
657
658 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$ in place,
659 /// rounding the result to the specified precision and with the specified rounding mode. An
660 /// [`Ordering`] is returned, indicating whether the rounded value is less than, equal to, or
661 /// greater than the exact value.
662 ///
663 /// See [`RoundingMode`] for a description of the possible rounding modes.
664 ///
665 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
666 /// overflow, and underflow.
667 ///
668 /// If you know you'll be using `Nearest`, consider using [`Float::compound_prec_assign`]
669 /// instead. If you know that your target precision is the precision of the input, consider
670 /// using [`Float::compound_round_assign`] instead. If both of these things are true, consider
671 /// using the [`CompoundAssign`] trait instead.
672 ///
673 /// # Worst-case complexity
674 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
675 ///
676 /// $M(n) = O(n \log n)$
677 ///
678 /// where $T$ is time, $M$ is additional memory, $n$ is `max(prec, self.significant_bits())`,
679 /// and $m$ is the number of significant bits of the exponent `n`.
680 ///
681 /// # Panics
682 /// Panics if `prec` is zero, or if `rm` is `Exact` but the result cannot be represented exactly
683 /// with the given precision.
684 ///
685 /// # Examples
686 /// ```
687 /// use malachite_base::rounding_modes::RoundingMode::*;
688 /// use malachite_float::Float;
689 /// use std::cmp::Ordering::*;
690 ///
691 /// let mut x = Float::from(3);
692 /// assert_eq!(x.compound_prec_round_assign(2, 10, Floor), Equal);
693 /// assert_eq!(x.to_string(), "16.000");
694 /// ```
695 pub fn compound_prec_round_assign(&mut self, n: i64, prec: u64, rm: RoundingMode) -> Ordering {
696 let (y, o) = self.compound_prec_round_ref(n, prec, rm);
697 *self = y;
698 o
699 }
700
701 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$ in place,
702 /// rounding the result to the nearest value of the specified precision. An [`Ordering`] is
703 /// returned, indicating whether the rounded value is less than, equal to, or greater than the
704 /// exact value.
705 ///
706 /// If the compound value is equidistant from two [`Float`]s with the specified precision, the
707 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
708 /// description of the `Nearest` rounding mode.
709 ///
710 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
711 /// overflow, and underflow.
712 ///
713 /// If you want to use a rounding mode other than `Nearest`, consider using
714 /// [`Float::compound_prec_round_assign`] instead. If you know that your target precision is the
715 /// precision of the input, consider using the [`CompoundAssign`] trait instead.
716 ///
717 /// # Worst-case complexity
718 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
719 ///
720 /// $M(n) = O(n \log n)$
721 ///
722 /// where $T$ is time, $M$ is additional memory, $n$ is `max(prec, self.significant_bits())`,
723 /// and $m$ is the number of significant bits of the exponent `n`.
724 ///
725 /// # Panics
726 /// Panics if `prec` is zero.
727 ///
728 /// # Examples
729 /// ```
730 /// use malachite_base::num::basic::traits::Two;
731 /// use malachite_float::Float;
732 /// use std::cmp::Ordering::*;
733 ///
734 /// let mut x = Float::TWO;
735 /// assert_eq!(x.compound_prec_assign(-2, 10), Less);
736 /// assert_eq!(x.to_string(), "0.11108");
737 /// ```
738 #[inline]
739 pub fn compound_prec_assign(&mut self, n: i64, prec: u64) -> Ordering {
740 self.compound_prec_round_assign(n, prec, Nearest)
741 }
742
743 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$ in place,
744 /// rounding the result to the precision of the input with the specified rounding mode. An
745 /// [`Ordering`] is returned, indicating whether the rounded value is less than, equal to, or
746 /// greater than the exact value.
747 ///
748 /// See [`RoundingMode`] for a description of the possible rounding modes.
749 ///
750 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
751 /// overflow, and underflow.
752 ///
753 /// If you know you'll be using `Nearest`, consider using the [`CompoundAssign`] trait instead.
754 /// If you want to specify an output precision, consider using
755 /// [`Float::compound_prec_round_assign`] instead.
756 ///
757 /// # Worst-case complexity
758 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
759 ///
760 /// $M(n) = O(n \log n)$
761 ///
762 /// where $T$ is time, $M$ is additional memory, $n$ is `self.significant_bits()`, and $m$ is
763 /// the number of significant bits of the exponent `n`.
764 ///
765 /// # Panics
766 /// Panics if `rm` is `Exact` but the result cannot be represented exactly with the precision of
767 /// the input.
768 ///
769 /// # Examples
770 /// ```
771 /// use malachite_base::rounding_modes::RoundingMode::*;
772 /// use malachite_float::Float;
773 /// use std::cmp::Ordering::*;
774 ///
775 /// let mut x = Float::from(1.5);
776 /// assert_eq!(x.compound_round_assign(2, Ceiling), Greater);
777 /// assert_eq!(x.to_string(), "8.0");
778 /// ```
779 #[inline]
780 pub fn compound_round_assign(&mut self, n: i64, rm: RoundingMode) -> Ordering {
781 let prec = self.significant_bits();
782 self.compound_prec_round_assign(n, prec, rm)
783 }
784}
785
786impl Compound<i64> for Float {
787 type Output = Self;
788
789 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
790 /// result to the nearest value with the precision of the input. The [`Float`] is taken by
791 /// value.
792 ///
793 /// The compound function is defined in IEEE 754 and is useful for computing compound interest:
794 /// if $x$ is an interest rate, then $(1+x)^n$ is the factor by which a principal grows after
795 /// $n$ compounding periods.
796 ///
797 /// If the compound value is equidistant from two [`Float`]s with the specified precision, the
798 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
799 /// description of the `Nearest` rounding mode.
800 ///
801 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
802 /// overflow, and underflow.
803 ///
804 /// # Worst-case complexity
805 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
806 ///
807 /// $M(n) = O(n \log n)$
808 ///
809 /// where $T$ is time, $M$ is additional memory, $n$ is `self.significant_bits()`, and $m$ is
810 /// the number of significant bits of the exponent `n`.
811 ///
812 /// # Examples
813 /// ```
814 /// use malachite_base::num::arithmetic::traits::Compound;
815 /// use malachite_base::num::basic::traits::Two;
816 /// use malachite_float::Float;
817 ///
818 /// assert_eq!(
819 /// Float::from(0.1).compound(10).to_string(),
820 /// "2.5937424601000005"
821 /// );
822 /// assert_eq!(Float::from(3).compound(2).to_string(), "16.0");
823 /// assert_eq!(Float::TWO.compound(-2).to_string(), "0.12");
824 /// ```
825 #[inline]
826 fn compound(self, n: i64) -> Self {
827 let prec = self.significant_bits();
828 self.compound_prec_round(n, prec, Nearest).0
829 }
830}
831
832impl Compound<i64> for &Float {
833 type Output = Float;
834
835 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$, rounding the
836 /// result to the nearest value with the precision of the input. The [`Float`] is taken by
837 /// reference.
838 ///
839 /// The compound function is defined in IEEE 754 and is useful for computing compound interest:
840 /// if $x$ is an interest rate, then $(1+x)^n$ is the factor by which a principal grows after
841 /// $n$ compounding periods.
842 ///
843 /// If the compound value is equidistant from two [`Float`]s with the specified precision, the
844 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
845 /// description of the `Nearest` rounding mode.
846 ///
847 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
848 /// overflow, and underflow.
849 ///
850 /// # Worst-case complexity
851 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
852 ///
853 /// $M(n) = O(n \log n)$
854 ///
855 /// where $T$ is time, $M$ is additional memory, $n$ is `self.significant_bits()`, and $m$ is
856 /// the number of significant bits of the exponent `n`.
857 ///
858 /// # Examples
859 /// ```
860 /// use malachite_base::num::arithmetic::traits::Compound;
861 /// use malachite_base::num::basic::traits::Two;
862 /// use malachite_float::Float;
863 ///
864 /// assert_eq!(
865 /// (&Float::from(0.1)).compound(10).to_string(),
866 /// "2.5937424601000005"
867 /// );
868 /// assert_eq!((&Float::from(3)).compound(2).to_string(), "16.0");
869 /// assert_eq!((&Float::TWO).compound(-2).to_string(), "0.12");
870 /// ```
871 #[inline]
872 fn compound(self, n: i64) -> Float {
873 let prec = self.significant_bits();
874 self.compound_prec_round_ref(n, prec, Nearest).0
875 }
876}
877
878impl CompoundAssign<i64> for Float {
879 /// Computes the compound function $(1+x)^n$ of a [`Float`] $x$ and an [`i64`] $n$ in place,
880 /// rounding the result to the nearest value with the precision of the input.
881 ///
882 /// The compound function is defined in IEEE 754 and is useful for computing compound interest:
883 /// if $x$ is an interest rate, then $(1+x)^n$ is the factor by which a principal grows after
884 /// $n$ compounding periods.
885 ///
886 /// If the compound value is equidistant from two [`Float`]s with the specified precision, the
887 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
888 /// description of the `Nearest` rounding mode.
889 ///
890 /// See the [`Float::compound_prec_round`] documentation for information on special cases,
891 /// overflow, and underflow.
892 ///
893 /// # Worst-case complexity
894 /// $T(n, m) = O(mn^{3/2} \log n \log\log n)$
895 ///
896 /// $M(n) = O(n \log n)$
897 ///
898 /// where $T$ is time, $M$ is additional memory, $n$ is `self.significant_bits()`, and $m$ is
899 /// the number of significant bits of the exponent `n`.
900 ///
901 /// # Examples
902 /// ```
903 /// use malachite_base::num::arithmetic::traits::CompoundAssign;
904 /// use malachite_float::Float;
905 ///
906 /// let mut x = Float::from(0.1);
907 /// x.compound_assign(10);
908 /// assert_eq!(x.to_string(), "2.5937424601000005");
909 /// ```
910 #[inline]
911 fn compound_assign(&mut self, n: i64) {
912 let prec = self.significant_bits();
913 self.compound_prec_round_assign(n, prec, Nearest);
914 }
915}
916
917/// Computes the compound function $(1+x)^n$ of a primitive float and an [`i64`], returning a
918/// primitive float.
919///
920/// The result is correctly rounded to the nearest value.
921///
922/// $$
923/// f(x,n) = (1+x)^n+\varepsilon.
924/// $$
925/// - If $(1+x)^n$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or assumed to be 0.
926/// - If $(1+x)^n$ is finite and nonzero, then $|\varepsilon| < 2^{\lfloor\log_2 (1+x)^n\rfloor-p}$,
927/// where $p$ is the precision of the output (typically 24 if `T` is a [`f32`] and 53 if `T` is a
928/// [`f64`], but less if the output is subnormal).
929///
930/// Special cases:
931/// - $f(\text{NaN},n)=\text{NaN}$ if $n\neq 0$, and $1.0$ if $n=0$
932/// - $f(-\infty,n)=\text{NaN}$, even if $n=0$
933/// - $f(\infty,0)=1.0$
934/// - $f(\infty,n)=\infty$ if $n>0$, and $0.0$ if $n<0$
935/// - $f(\pm 0.0,n)=1.0$
936/// - $f(x,n)=\text{NaN}$ if $x<-1$, even if $n=0$
937/// - $f(-1.0,n)=1.0$ if $n=0$, $0.0$ if $n>0$, and $\infty$ if $n<0$
938/// - $f(x,0)=1.0$ if $x\geq -1$
939///
940/// The result is never negative. If the result overflows, $\infty$ is returned, and if it
941/// underflows, $0.0$ is returned.
942///
943/// # Worst-case complexity
944/// Constant time and additional memory.
945///
946/// # Examples
947/// ```
948/// use malachite_base::num::float::NiceFloat;
949/// use malachite_float::float::arithmetic::compound::primitive_float_compound;
950///
951/// assert_eq!(NiceFloat(primitive_float_compound(0.5, 2)), NiceFloat(2.25));
952/// assert_eq!(
953/// NiceFloat(primitive_float_compound(0.1, 10)),
954/// NiceFloat(2.5937424601)
955/// );
956/// assert_eq!(
957/// NiceFloat(primitive_float_compound(-0.5, -2)),
958/// NiceFloat(4.0)
959/// );
960/// ```
961#[allow(clippy::type_repetition_in_bounds)]
962#[inline]
963pub fn primitive_float_compound<T: PrimitiveFloat>(x: T, n: i64) -> T
964where
965 Float: From<T> + PartialOrd<T>,
966 for<'a> T: ExactFrom<&'a Float> + RoundingFrom<&'a Float>,
967{
968 emulate_float_to_float_fn(|x, prec| x.compound_prec(n, prec), x)
969}