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}