malachite_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_extras::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)^2)$
44 ///
45 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
46 ///
47 /// # Panics
48 /// Panics if `prec` is zero or if `rm` is `Exact`.
49 ///
50 /// # Examples
51 /// ```
52 /// use malachite_base::rounding_modes::RoundingMode::*;
53 /// use malachite_float::Float;
54 /// use std::cmp::Ordering::*;
55 ///
56 /// let (pi, o) = Float::pi_prec_round(100, Floor);
57 /// assert_eq!(pi.to_string(), "3.141592653589793238462643383279");
58 /// assert_eq!(o, Less);
59 ///
60 /// let (pi, o) = Float::pi_prec_round(100, Ceiling);
61 /// assert_eq!(pi.to_string(), "3.141592653589793238462643383282");
62 /// assert_eq!(o, Greater);
63 /// ```
64 ///
65 // This is mpfr_const_pi_internal from const_pi.c, MPFR 4.2.0.
66 #[inline]
67 pub fn pi_prec_round(prec: u64, rm: RoundingMode) -> (Self, Ordering) {
68 // we need 9 * 2 ^ kmax - 4 >= px + 2 * kmax + 8
69 let mut kmax = 2;
70 while ((prec + 2 * kmax + 12) / 9) >> kmax != 0 {
71 kmax += 1;
72 }
73 // guarantees no recomputation for px <= 10000
74 let mut working_prec = prec + 3 * kmax + 14;
75 let mut increment = Limb::WIDTH;
76 loop {
77 let mut a = Self::one_prec(working_prec);
78 let mut big_a = a.clone();
79 let mut big_b = Self::one_half_prec(working_prec);
80 let mut big_d = Self::one_prec(working_prec) >> 2;
81 let mut k = 0;
82 loop {
83 let s = (&big_a + &big_b) >> 2;
84 a = (a + big_b.sqrt()) >> 1;
85 big_a = (&a).square();
86 big_b = (&big_a - s) << 1;
87 let mut s = &big_a - &big_b;
88 assert!(s < 1u32);
89 let ip = i64::exact_from(working_prec);
90 let cancel = if s == 0 {
91 ip
92 } else {
93 i64::from(-s.get_exponent().unwrap())
94 };
95 s <<= k;
96 big_d -= s;
97 // stop when |A_k - B_k| <= 2 ^ (k - p) i.e. cancel >= p - k
98 if cancel >= ip - i64::exact_from(k) {
99 break;
100 }
101 k += 1;
102 }
103 let pi: Self = big_b / big_d;
104 if float_can_round(
105 pi.significand_ref().unwrap(),
106 working_prec - (k << 1) - 8,
107 prec,
108 rm,
109 ) {
110 return Self::from_float_prec_round(pi, prec, rm);
111 }
112 working_prec += kmax + increment;
113 increment = working_prec >> 1;
114 }
115 }
116
117 /// Returns an approximation of $\pi$, with the given precision and rounded to the nearest
118 /// [`Float`] of that precision. An [`Ordering`] is also returned, indicating whether the
119 /// rounded value is less than or greater than the exact value of the constant. (Since the
120 /// constant is irrational, the rounded value is never equal to the exact value.)
121 ///
122 /// $$
123 /// x = \pi+\varepsilon.
124 /// $$
125 /// - $|\varepsilon| < 2^{-p+1}$.
126 ///
127 /// The constant is irrational and transcendental.
128 ///
129 /// The output has precision `prec`.
130 ///
131 /// # Worst-case complexity
132 /// $T(n) = O(n (\log n)^2 \log\log n)$
133 ///
134 /// $M(n) = O(n (\log n)^2)$
135 ///
136 /// where $T$ is time, $M$ is additional memory, and $n$ is `prec`.
137 ///
138 /// # Panics
139 /// Panics if `prec` is zero.
140 ///
141 /// # Examples
142 /// ```
143 /// use malachite_float::Float;
144 /// use std::cmp::Ordering::*;
145 ///
146 /// let (pi, o) = Float::pi_prec(1);
147 /// assert_eq!(pi.to_string(), "4.0");
148 /// assert_eq!(o, Greater);
149 ///
150 /// let (pi, o) = Float::pi_prec(10);
151 /// assert_eq!(pi.to_string(), "3.141");
152 /// assert_eq!(o, Less);
153 ///
154 /// let (pi, o) = Float::pi_prec(100);
155 /// assert_eq!(pi.to_string(), "3.141592653589793238462643383279");
156 /// assert_eq!(o, Less);
157 /// ```
158 #[inline]
159 pub fn pi_prec(prec: u64) -> (Self, Ordering) {
160 Self::pi_prec_round(prec, Nearest)
161 }
162}