Skip to main content

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}