Skip to main content

malachite_float/constants/
ln_2.rs

1// Copyright © 2026 Mikhail Hogrefe
2//
3// Uses code adopted from the GNU MPFR Library.
4//
5//      Copyright 1999, 2001-2024 Free Software Foundation, Inc.
6//
7//      Contributed by the AriC 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 alloc::vec;
17use core::cmp::Ordering;
18use core::mem::swap;
19use malachite_base::num::arithmetic::traits::CeilingLogBase2;
20use malachite_base::num::basic::integers::PrimitiveInt;
21use malachite_base::num::basic::traits::{One, Zero};
22use malachite_base::num::conversion::traits::WrappingFrom;
23use malachite_base::rounding_modes::RoundingMode::{self, *};
24use malachite_nz::integer::Integer;
25use malachite_nz::natural::arithmetic::float_extras::float_can_round;
26use malachite_nz::platform::Limb;
27
28// Auxiliary function: Compute the terms from n1 to n2 (excluded) 3 / 4 * sum((-1) ^ n * n! ^ 2 / 2
29// ^ n / (2 * n + 1)!, n = n1...n2 - 1).s
30//
31// Numerator is T[0], denominator is Q[0], Compute P[0] only when need_P is non-zero.
32//
33// Need 1 + ceil(log(n2 - n1) / log(2)) cells in T[], P[], Q[].
34//
35// This is S from const_log2.c, MPFR 4.2.0.
36fn sum(t: &mut [Integer], p: &mut [Integer], q: &mut [Integer], n1: u64, n2: u64, need_p: bool) {
37    if n2 == n1 + 1 {
38        p[0] = if n1 == 0 {
39            const { Integer::const_from_unsigned(3) }
40        } else {
41            -Integer::from(n1)
42        };
43        q[0] = ((Integer::from(n1) << 1u32) + Integer::ONE) << 2u32;
44        t[0].clone_from(&p[0]);
45    } else {
46        let m = (n1 >> 1) + (n2 >> 1) + (n1 & 1 & n2);
47        sum(t, p, q, n1, m, true);
48        let (t_head, t_tail) = t.split_first_mut().unwrap();
49        let (p_head, p_tail) = p.split_first_mut().unwrap();
50        let (q_head, q_tail) = q.split_first_mut().unwrap();
51        sum(t_tail, p_tail, q_tail, m, n2, need_p);
52        *t_head *= &q_tail[0];
53        t_tail[0] *= &*p_head;
54        *t_head += &t_tail[0];
55        if need_p {
56            *p_head *= &p_tail[0];
57        }
58        *q_head *= &q_tail[0];
59        // remove common trailing zeros if any
60        let mut tz = t_head.trailing_zeros().unwrap();
61        if tz != 0 {
62            let mut qz = q_head.trailing_zeros().unwrap();
63            if qz < tz {
64                tz = qz;
65            }
66            if need_p {
67                qz = p_head.trailing_zeros().unwrap();
68                if qz < tz {
69                    tz = qz;
70                }
71            }
72            // now tz = min(val(T), val(Q), val(P))
73            if tz != 0 {
74                *t_head >>= tz;
75                *q_head >>= tz;
76                if need_p {
77                    *p_head >>= tz;
78                }
79            }
80        }
81    }
82}
83
84impl Float {
85    /// Returns an approximation of the natural logarithm of 2, with the given precision and rounded
86    /// using the given [`RoundingMode`]. An [`Ordering`] is also returned, indicating whether the
87    /// rounded value is less than or greater than the exact value of the constant. (Since the
88    /// constant is irrational, the rounded value is never equal to the exact value.)
89    ///
90    /// $$
91    /// x = \ln 2+\varepsilon.
92    /// $$
93    /// - If $m$ is not `Nearest`, then $|\varepsilon| < 2^{-p}$.
94    /// - If $m$ is `Nearest`, then $|\varepsilon| < 2^{-p-1}$.
95    ///
96    /// The constant is irrational and transcendental.
97    ///
98    /// The output has precision `prec`.
99    ///
100    /// # Worst-case complexity
101    /// $T(n) = O(n (\log n)^2 \log\log n)$
102    ///
103    /// $M(n) = O(n (\log n)^2)$
104    ///
105    /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
106    ///
107    /// # Panics
108    /// Panics if `prec` is zero or if `rm` is `Exact`.
109    ///
110    /// # Examples
111    /// ```
112    /// use malachite_base::rounding_modes::RoundingMode::*;
113    /// use malachite_float::Float;
114    /// use std::cmp::Ordering::*;
115    ///
116    /// let (l2, o) = Float::ln_2_prec_round(100, Floor);
117    /// assert_eq!(l2.to_string(), "0.693147180559945309417232121458");
118    /// assert_eq!(o, Less);
119    ///
120    /// let (l2, o) = Float::ln_2_prec_round(100, Ceiling);
121    /// assert_eq!(l2.to_string(), "0.693147180559945309417232121459");
122    /// assert_eq!(o, Greater);
123    /// ```
124    ///
125    /// This is mpfr_const_log2_internal from const_log2.c, MPFR 4.2.0.
126    pub fn ln_2_prec_round(prec: u64, rm: RoundingMode) -> (Self, Ordering) {
127        let mut working_prec = prec + prec.ceiling_log_base_2() + 3;
128        let mut increment = Limb::WIDTH;
129        loop {
130            let big_n = working_prec / 3 + 1;
131            // the following are needed for error analysis (see algorithms.tex)
132            assert!(working_prec >= 3 && big_n >= 2);
133            let lg_big_n = usize::wrapping_from(big_n.ceiling_log_base_2()) + 1;
134            let mut scratch = vec![Integer::ZERO; 3 * lg_big_n];
135            split_into_chunks_mut!(scratch, lg_big_n, [t, p], q);
136            sum(t, p, q, 0, big_n, false);
137            let mut t0 = Integer::ZERO;
138            let mut q0 = Integer::ZERO;
139            swap(&mut t0, &mut t[0]);
140            swap(&mut q0, &mut q[0]);
141            let ln_2 = Self::from_integer_prec(t0, working_prec).0
142                / Self::from_integer_prec(q0, working_prec).0;
143            if float_can_round(ln_2.significand_ref().unwrap(), working_prec - 2, prec, rm) {
144                return Self::from_float_prec_round(ln_2, prec, rm);
145            }
146            working_prec += increment;
147            increment = working_prec >> 1;
148        }
149    }
150
151    /// Returns an approximation of the natural logarithm of 2, with the given precision and rounded
152    /// to the nearest [`Float`] of that precision. An [`Ordering`] is also returned, indicating
153    /// whether the rounded value is less than or greater than the exact value of the constant.
154    /// (Since the constant is irrational, the rounded value is never equal to the exact value.)
155    ///
156    /// $$
157    /// x = \ln 2+\varepsilon.
158    /// $$
159    /// - $|\varepsilon| < 2^{-p-1}$.
160    ///
161    /// The constant is irrational and transcendental.
162    ///
163    /// The output has precision `prec`.
164    ///
165    /// # Worst-case complexity
166    /// $T(n) = O(n (\log n)^2 \log\log n)$
167    ///
168    /// $M(n) = O(n (\log n)^2)$
169    ///
170    /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
171    ///
172    /// # Panics
173    /// Panics if `prec` is zero.
174    ///
175    /// # Examples
176    /// ```
177    /// use malachite_float::Float;
178    /// use std::cmp::Ordering::*;
179    ///
180    /// let (l2, o) = Float::ln_2_prec(1);
181    /// assert_eq!(l2.to_string(), "0.5");
182    /// assert_eq!(o, Less);
183    ///
184    /// let (l2, o) = Float::ln_2_prec(10);
185    /// assert_eq!(l2.to_string(), "0.693");
186    /// assert_eq!(o, Greater);
187    ///
188    /// let (l2, o) = Float::ln_2_prec(100);
189    /// assert_eq!(l2.to_string(), "0.693147180559945309417232121458");
190    /// assert_eq!(o, Less);
191    /// ```
192    #[inline]
193    pub fn ln_2_prec(prec: u64) -> (Self, Ordering) {
194        Self::ln_2_prec_round(prec, Nearest)
195    }
196}