malachite_float/float/constants/pi.rs
1// Copyright © 2026 Mikhail Hogrefe
2//
3// Uses code adopted from the GNU MPFR Library.
4//
5// Copyright 1999-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 core::cmp::Ordering;
17use malachite_base::num::arithmetic::traits::{Sqrt, Square};
18use malachite_base::num::basic::integers::PrimitiveInt;
19use malachite_base::num::conversion::traits::ExactFrom;
20use malachite_base::rounding_modes::RoundingMode::{self, *};
21use malachite_nz::natural::arithmetic::float::round::float_can_round;
22use malachite_nz::platform::Limb;
23
24impl Float {
25 /// Returns an approximation of $\pi$, with the given precision and rounded using the given
26 /// [`RoundingMode`]. An [`Ordering`] is also returned, indicating whether the rounded value is
27 /// less than or greater than the exact value of the constant. (Since the constant is
28 /// irrational, the rounded value is never equal to the exact value.)
29 ///
30 /// $$
31 /// x = \pi+\varepsilon.
32 /// $$
33 /// - If $m$ is not `Nearest`, then $|\varepsilon| < 2^{-p+2}$.
34 /// - If $m$ is `Nearest`, then $|\varepsilon| < 2^{-p+1}$.
35 ///
36 /// The constant is irrational and transcendental.
37 ///
38 /// The output has precision `prec`.
39 ///
40 /// # Worst-case complexity
41 /// $T(n) = O(n (\log n)^2 \log\log n)$
42 ///
43 /// $M(n) = O(n \log n)$
44 ///
45 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`: the iteration holds a
46 /// bounded number of working-precision values at a time (freeing each round's temporaries), so
47 /// the peak is the internal memory of a single full-precision multiplication.
48 ///
49 /// # Panics
50 /// Panics if `prec` is zero or if `rm` is `Exact`.
51 ///
52 /// # Examples
53 /// ```
54 /// use malachite_base::rounding_modes::RoundingMode::*;
55 /// use malachite_float::Float;
56 /// use std::cmp::Ordering::*;
57 ///
58 /// let (pi, o) = Float::pi_prec_round(100, Floor);
59 /// assert_eq!(pi.to_string(), "3.1415926535897932384626433832793");
60 /// assert_eq!(o, Less);
61 ///
62 /// let (pi, o) = Float::pi_prec_round(100, Ceiling);
63 /// assert_eq!(pi.to_string(), "3.1415926535897932384626433832825");
64 /// assert_eq!(o, Greater);
65 /// ```
66 ///
67 // This is mpfr_const_pi_internal from const_pi.c, MPFR 4.2.0.
68 #[inline]
69 pub fn pi_prec_round(prec: u64, rm: RoundingMode) -> (Self, Ordering) {
70 // we need 9 * 2 ^ kmax - 4 >= px + 2 * kmax + 8
71 let mut kmax = 2;
72 while ((prec + (kmax << 1) + 12) / 9) >> kmax != 0 {
73 kmax += 1;
74 }
75 // guarantees no recomputation for px <= 10000
76 let mut working_prec = prec + 3 * kmax + 14;
77 let mut increment = Limb::WIDTH;
78 loop {
79 let mut a = Self::one_prec(working_prec);
80 let mut big_a = a.clone();
81 let mut big_b = Self::one_half_prec(working_prec);
82 let mut big_d = Self::one_prec(working_prec) >> 2u32;
83 let mut k = 0;
84 loop {
85 let s = (&big_a + &big_b) >> 2u32;
86 a = (a + big_b.sqrt()) >> 1u32;
87 big_a = (&a).square();
88 big_b = (&big_a - s) << 1u32;
89 let mut s = &big_a - &big_b;
90 assert!(s < 1u32);
91 let ip = i64::exact_from(working_prec);
92 let cancel = if s == 0u32 {
93 ip
94 } else {
95 i64::from(-s.get_exponent().unwrap())
96 };
97 s <<= k;
98 big_d -= s;
99 // stop when |A_k - B_k| <= 2 ^ (k - p) i.e. cancel >= p - k
100 if cancel >= ip - i64::exact_from(k) {
101 break;
102 }
103 k += 1;
104 }
105 let pi: Self = big_b / big_d;
106 if float_can_round(
107 pi.significand_ref().unwrap(),
108 working_prec - (k << 1) - 8,
109 prec,
110 rm,
111 ) {
112 return Self::from_float_prec_round(pi, prec, rm);
113 }
114 working_prec += kmax + increment;
115 increment = working_prec >> 1;
116 }
117 }
118
119 /// Returns an approximation of $\pi$, with the given precision and rounded to the nearest
120 /// [`Float`] of that precision. An [`Ordering`] is also returned, indicating whether the
121 /// rounded value is less than or greater than the exact value of the constant. (Since the
122 /// constant is irrational, the rounded value is never equal to the exact value.)
123 ///
124 /// $$
125 /// x = \pi+\varepsilon.
126 /// $$
127 /// - $|\varepsilon| < 2^{-p+1}$.
128 ///
129 /// The constant is irrational and transcendental.
130 ///
131 /// The output has precision `prec`.
132 ///
133 /// # Worst-case complexity
134 /// $T(n) = O(n (\log n)^2 \log\log n)$
135 ///
136 /// $M(n) = O(n \log n)$
137 ///
138 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`: the iteration holds a
139 /// bounded number of working-precision values at a time (freeing each round's temporaries), so
140 /// the peak is the internal memory of a single full-precision multiplication.
141 ///
142 /// # Panics
143 /// Panics if `prec` is zero.
144 ///
145 /// # Examples
146 /// ```
147 /// use malachite_float::Float;
148 /// use std::cmp::Ordering::*;
149 ///
150 /// let (pi, o) = Float::pi_prec(1);
151 /// assert_eq!(pi.to_string(), "4.0");
152 /// assert_eq!(o, Greater);
153 ///
154 /// let (pi, o) = Float::pi_prec(10);
155 /// assert_eq!(pi.to_string(), "3.1406");
156 /// assert_eq!(o, Less);
157 ///
158 /// let (pi, o) = Float::pi_prec(100);
159 /// assert_eq!(pi.to_string(), "3.1415926535897932384626433832793");
160 /// assert_eq!(o, Less);
161 /// ```
162 #[inline]
163 pub fn pi_prec(prec: u64) -> (Self, Ordering) {
164 Self::pi_prec_round(prec, Nearest)
165 }
166}