malachite_float/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 crate::float::basic::extended::ExtendedFloat;
17use alloc::vec;
18use core::cmp::Ordering;
19use core::mem::swap;
20use malachite_base::num::arithmetic::traits::CeilingLogBase2;
21use malachite_base::num::basic::integers::PrimitiveInt;
22use malachite_base::num::basic::traits::{One, Zero};
23use malachite_base::num::conversion::traits::WrappingFrom;
24use malachite_base::num::logic::traits::SignificantBits;
25use malachite_base::rounding_modes::RoundingMode::{self, *};
26use malachite_nz::integer::Integer;
27use malachite_nz::natural::arithmetic::float_extras::float_can_round;
28use malachite_nz::platform::Limb;
29use malachite_q::Rational;
30
31// Auxiliary function: Compute the terms from n1 to n2 (excluded) 3 / 4 * sum((-1) ^ n * n! ^ 2 / 2
32// ^ n / (2 * n + 1)!, n = n1...n2 - 1).s
33//
34// Numerator is T[0], denominator is Q[0], Compute P[0] only when need_P is non-zero.
35//
36// Need 1 + ceil(log(n2 - n1) / log(2)) cells in T[], P[], Q[].
37//
38// This is S from const_log2.c, MPFR 4.2.0.
39fn sum(t: &mut [Integer], p: &mut [Integer], q: &mut [Integer], n1: u64, n2: u64, need_p: bool) {
40 if n2 == n1 + 1 {
41 p[0] = if n1 == 0 {
42 const { Integer::const_from_unsigned(3) }
43 } else {
44 -Integer::from(n1)
45 };
46 q[0] = ((Integer::from(n1) << 1u32) + Integer::ONE) << 2u32;
47 t[0].clone_from(&p[0]);
48 } else {
49 let m = (n1 >> 1) + (n2 >> 1) + (n1 & 1 & n2);
50 sum(t, p, q, n1, m, true);
51 let (t_head, t_tail) = t.split_first_mut().unwrap();
52 let (p_head, p_tail) = p.split_first_mut().unwrap();
53 let (q_head, q_tail) = q.split_first_mut().unwrap();
54 sum(t_tail, p_tail, q_tail, m, n2, need_p);
55 *t_head *= &q_tail[0];
56 t_tail[0] *= &*p_head;
57 *t_head += &t_tail[0];
58 if need_p {
59 *p_head *= &p_tail[0];
60 }
61 *q_head *= &q_tail[0];
62 // remove common trailing zeros if any
63 let mut tz = t_head.trailing_zeros().unwrap();
64 if tz != 0 {
65 let mut qz = q_head.trailing_zeros().unwrap();
66 if qz < tz {
67 tz = qz;
68 }
69 if need_p {
70 qz = p_head.trailing_zeros().unwrap();
71 if qz < tz {
72 tz = qz;
73 }
74 }
75 // now tz = min(val(T), val(Q), val(P))
76 if tz != 0 {
77 *t_head >>= tz;
78 *q_head >>= tz;
79 if need_p {
80 *p_head >>= tz;
81 }
82 }
83 }
84 }
85}
86
87impl Float {
88 /// Returns an approximation of the natural logarithm of 2, with the given precision and rounded
89 /// using the given [`RoundingMode`]. An [`Ordering`] is also returned, indicating whether the
90 /// rounded value is less than or greater than the exact value of the constant. (Since the
91 /// constant is irrational, the rounded value is never equal to the exact value.)
92 ///
93 /// $$
94 /// x = \ln 2+\varepsilon.
95 /// $$
96 /// - If $m$ is not `Nearest`, then $|\varepsilon| < 2^{-p}$.
97 /// - If $m$ is `Nearest`, then $|\varepsilon| < 2^{-p-1}$.
98 ///
99 /// The constant is irrational and transcendental.
100 ///
101 /// The output has precision `prec`.
102 ///
103 /// # Worst-case complexity
104 /// $T(n) = O(n (\log n)^2 \log\log n)$
105 ///
106 /// $M(n) = O(n (\log n)^2)$
107 ///
108 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
109 ///
110 /// # Panics
111 /// Panics if `prec` is zero or if `rm` is `Exact`.
112 ///
113 /// # Examples
114 /// ```
115 /// use malachite_base::rounding_modes::RoundingMode::*;
116 /// use malachite_float::Float;
117 /// use std::cmp::Ordering::*;
118 ///
119 /// let (l2, o) = Float::ln_2_prec_round(100, Floor);
120 /// assert_eq!(l2.to_string(), "0.69314718055994530941723212145798");
121 /// assert_eq!(o, Less);
122 ///
123 /// let (l2, o) = Float::ln_2_prec_round(100, Ceiling);
124 /// assert_eq!(l2.to_string(), "0.69314718055994530941723212145877");
125 /// assert_eq!(o, Greater);
126 /// ```
127 ///
128 /// This is mpfr_const_log2_internal from const_log2.c, MPFR 4.2.0.
129 pub fn ln_2_prec_round(prec: u64, rm: RoundingMode) -> (Self, Ordering) {
130 let mut working_prec = prec + prec.ceiling_log_base_2() + 3;
131 let mut increment = Limb::WIDTH;
132 loop {
133 let big_n = working_prec / 3 + 1;
134 // the following are needed for error analysis (see algorithms.tex)
135 assert!(working_prec >= 3 && big_n >= 2);
136 let lg_big_n = usize::wrapping_from(big_n.ceiling_log_base_2()) + 1;
137 let mut scratch = vec![Integer::ZERO; 3 * lg_big_n];
138 split_into_chunks_mut!(scratch, lg_big_n, [t, p], q);
139 sum(t, p, q, 0, big_n, false);
140 let mut t0 = Integer::ZERO;
141 let mut q0 = Integer::ZERO;
142 swap(&mut t0, &mut t[0]);
143 swap(&mut q0, &mut q[0]);
144 // ln(2) = t0 / q0 ~ 0.69. For large `working_prec`, t0 and q0 each have more than
145 // `MAX_EXPONENT` bits, so `Float::from_integer_prec` would overflow them to infinity.
146 // Round each into an `ExtendedFloat` (whose exponent is an `i64`) and divide there; the
147 // in-range quotient then converts back to a `Float`. First shift off the bits below
148 // `working_prec` (the same shift for both, so the ratio is preserved): without this,
149 // the `ExtendedFloat` conversion would build a power-of-2 denominator as wide as t0/q0
150 // themselves (up to ~1 GB), rather than ~`working_prec` bits.
151 let extra = t0
152 .significant_bits()
153 .max(q0.significant_bits())
154 .saturating_sub(working_prec + Limb::WIDTH);
155 let (t0, q0) = (t0 >> extra, q0 >> extra);
156 let ext_t =
157 ExtendedFloat::from_rational_prec_round(Rational::from(t0), working_prec, Nearest)
158 .0;
159 let ext_q =
160 ExtendedFloat::from_rational_prec_round(Rational::from(q0), working_prec, Nearest)
161 .0;
162 let ln_2 = Self::try_from(ext_t.div_prec_val_ref(&ext_q, working_prec).0).unwrap();
163 if float_can_round(ln_2.significand_ref().unwrap(), working_prec - 2, prec, rm) {
164 return Self::from_float_prec_round(ln_2, prec, rm);
165 }
166 working_prec += increment;
167 increment = working_prec >> 1;
168 }
169 }
170
171 /// Returns an approximation of the natural logarithm of 2, with the given precision and rounded
172 /// to the nearest [`Float`] of that precision. An [`Ordering`] is also returned, indicating
173 /// whether the rounded value is less than or greater than the exact value of the constant.
174 /// (Since the constant is irrational, the rounded value is never equal to the exact value.)
175 ///
176 /// $$
177 /// x = \ln 2+\varepsilon.
178 /// $$
179 /// - $|\varepsilon| < 2^{-p-1}$.
180 ///
181 /// The constant is irrational and transcendental.
182 ///
183 /// The output has precision `prec`.
184 ///
185 /// # Worst-case complexity
186 /// $T(n) = O(n (\log n)^2 \log\log n)$
187 ///
188 /// $M(n) = O(n (\log n)^2)$
189 ///
190 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
191 ///
192 /// # Panics
193 /// Panics if `prec` is zero.
194 ///
195 /// # Examples
196 /// ```
197 /// use malachite_float::Float;
198 /// use std::cmp::Ordering::*;
199 ///
200 /// let (l2, o) = Float::ln_2_prec(1);
201 /// assert_eq!(l2.to_string(), "0.50");
202 /// assert_eq!(o, Less);
203 ///
204 /// let (l2, o) = Float::ln_2_prec(10);
205 /// assert_eq!(l2.to_string(), "0.69336");
206 /// assert_eq!(o, Greater);
207 ///
208 /// let (l2, o) = Float::ln_2_prec(100);
209 /// assert_eq!(l2.to_string(), "0.69314718055994530941723212145798");
210 /// assert_eq!(o, Less);
211 /// ```
212 #[inline]
213 pub fn ln_2_prec(prec: u64) -> (Self, Ordering) {
214 Self::ln_2_prec_round(prec, Nearest)
215 }
216}