% Options for packages loaded elsewhere \PassOptionsToPackage{unicode}{hyperref} \PassOptionsToPackage{hyphens}{url} % \documentclass[ ]{report} \usepackage{amsmath,amssymb} \usepackage{iftex} \ifPDFTeX \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage{textcomp} % provide euro and other symbols \else % if luatex or xetex \usepackage{unicode-math} % this also loads fontspec \defaultfontfeatures{Scale=MatchLowercase} \defaultfontfeatures[\rmfamily]{Ligatures=TeX,Scale=1} \fi \usepackage{lmodern} \ifPDFTeX\else % xetex/luatex font selection \fi % Use upquote if available, for straight quotes in verbatim environments \IfFileExists{upquote.sty}{\usepackage{upquote}}{} \IfFileExists{microtype.sty}{% use microtype if available \usepackage[]{microtype} \UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts }{} \makeatletter \@ifundefined{KOMAClassName}{% if non-KOMA class \IfFileExists{parskip.sty}{% \usepackage{parskip} }{% else \setlength{\parindent}{0pt} \setlength{\parskip}{6pt plus 2pt minus 1pt}} }{% if KOMA class \KOMAoptions{parskip=half}} \makeatother \usepackage{xcolor} \usepackage[margin=2.0cm,a4paper]{geometry} \usepackage{color} \usepackage{fancyvrb} \newcommand{\VerbBar}{|} \newcommand{\VERB}{\Verb[commandchars=\\\{\}]} \DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}} % Add ',fontsize=\small' for more characters per line \newenvironment{Shaded}{}{} \newcommand{\AlertTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{#1}}} \newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}} \newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.49,0.56,0.16}{#1}} \newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}} \newcommand{\BuiltInTok}[1]{\textcolor[rgb]{0.00,0.50,0.00}{#1}} \newcommand{\CharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}} \newcommand{\CommentTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textit{#1}}} \newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}} \newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.53,0.00,0.00}{#1}} \newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{#1}}} \newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.56,0.13,0.00}{#1}} \newcommand{\DecValTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}} \newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.73,0.13,0.13}{\textit{#1}}} \newcommand{\ErrorTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{#1}}} \newcommand{\ExtensionTok}[1]{#1} \newcommand{\FloatTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{#1}} \newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.02,0.16,0.49}{#1}} \newcommand{\ImportTok}[1]{\textcolor[rgb]{0.00,0.50,0.00}{\textbf{#1}}} \newcommand{\InformationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}} \newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{#1}}} \newcommand{\NormalTok}[1]{#1} \newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.40,0.40,0.40}{#1}} \newcommand{\OtherTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{#1}} \newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.74,0.48,0.00}{#1}} \newcommand{\RegionMarkerTok}[1]{#1} \newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}} \newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.73,0.40,0.53}{#1}} \newcommand{\StringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}} \newcommand{\VariableTok}[1]{\textcolor[rgb]{0.10,0.09,0.49}{#1}} \newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{#1}} \newcommand{\WarningTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{#1}}}} \usepackage{longtable,booktabs,array} \usepackage{calc} % for calculating minipage widths % Correct order of tables after \paragraph or \subparagraph \usepackage{etoolbox} \makeatletter \patchcmd\longtable{\par}{\if@noskipsec\mbox{}\fi\par}{}{} \makeatother % Allow footnotes in longtable head/foot \IfFileExists{footnotehyper.sty}{\usepackage{footnotehyper}}{\usepackage{footnote}} \makesavenoteenv{longtable} \setlength{\emergencystretch}{3em} % prevent overfull lines \providecommand{\tightlist}{% \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}} \setcounter{secnumdepth}{-\maxdimen} % remove section numbering \usepackage{titlesec} \usepackage{fancyvrb} \usepackage{fvextra} \usepackage{enumitem} \usepackage{pdfpages} \usepackage{longtable} \usepackage{etoolbox} \usepackage{fontspec} \setmainfont{lmroman10-regular.otf}[ BoldFont = lmroman10-bold.otf, ItalicFont = lmroman10-italic.otf, BoldItalicFont = lmroman10-bolditalic.otf, OpticalSize = 0 ] \AtBeginEnvironment{longtable}{\fontsize{6}{8}\selectfont} \newcommand{\chapfnt}{\fontsize{19}{21}} \newcommand{\secfnt}{\fontsize{14}{17}} \newcommand{\ssecfnt}{\fontsize{12}{14}} \newcommand{\sectionbreak}{\clearpage} \newcommand{\llbracket}{[\![} \newcommand{\rrbracket}{]\!]} \titleformat{\chapter}[display] {\normalfont\chapfnt\bfseries}{\chaptertitlename\ \thechapter}{20pt}{\chapfnt} \titleformat{\section} {\normalfont\secfnt\bfseries}{\thesection}{1em}{} \titleformat{\subsection} {\normalfont\ssecfnt\bfseries}{\thesubsection}{1em}{} \titlespacing*{\chapter} {0pt}{50pt}{40pt} \titlespacing*{\section} {0pt}{3.5ex plus 1ex minus .2ex}{2.3ex plus .2ex} \titlespacing*{\subsection} {0pt}{3.25ex plus 1ex minus .2ex}{1.5ex plus .2ex} \DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\},fontsize=\scriptsize,frame=single,rulecolor=\color{lightgray},breaklines,samepage,label=\tiny{Code},labelposition=topline} \DefineVerbatimEnvironment{verbatim}{Verbatim}{commandchars=\\\{\},fontsize=\scriptsize,frame=single,rulecolor=\color{lightgray},breaklines,samepage,label=\tiny{Output},labelposition=topline,fontshape=it} \setlist{after=\bigskip} \let\OldRule\rule \renewcommand{\rule}[2]{\OldRule{0.0\linewidth}{#2}} \ifLuaTeX \usepackage{selnolig} % disable illegal ligatures \fi \usepackage{bookmark} \IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available \urlstyle{same} \hypersetup{ hidelinks, pdfcreator={LaTeX via pandoc}} \title{Numerical Calculation of Pi} \author{The Publicator using Qwen/Qwen3.8-27B-FP8} \date{} \begin{document} \maketitle { \setcounter{tocdepth}{2} \tableofcontents } \chapter{Numerical Calculation of Pi}\label{numerical-calculation-of-pi} \textbf{Abstract:} This paper presents a comprehensive study of numerical methods for calculating the mathematical constant pi. It begins by outlining the importance of pi in mathematics, science, and computing, and reviews the historical development of approximation techniques from classical geometric methods to modern computational algorithms. The paper then examines the mathematical foundations underlying pi calculation, including geometric, trigonometric, series-based, and integral representations. Both classical methods, such as Archimedes' polygon approach and the Gregory-Leibniz series, and modern high-precision algorithms, including Machin-like formulas and the Chudnovsky algorithm, are discussed. The paper also addresses practical implementation issues, convergence behavior, error analysis, and computational performance, with attention to precision handling, algorithmic efficiency, and optimization strategies. Finally, it explores applications of pi computation in numerical analysis, benchmarking, cryptography, and scientific computing, and concludes by comparing the strengths and limitations of different methods while suggesting directions for future research and implementation improvements. \section{1. Introduction}\label{introduction} \subsection{1.1 The Significance of Pi}\label{the-significance-of-pi} The constant \(\pi\) is one of the most fundamental quantities in mathematics. Defined as the ratio of a circle's circumference to its diameter, it first arises in elementary geometry, but its influence extends far beyond the study of circles. It appears in trigonometry, complex analysis, number theory, probability, statistics, and mathematical physics. In analysis, \(\pi\) is central to Fourier series, the Gaussian integral, and the theory of periodic functions. In physics and engineering, it occurs in wave equations, orbital mechanics, electromagnetism, and fluid dynamics. A particularly important property of \(\pi\) is that it is irrational and transcendental. This means that it cannot be expressed as a ratio of two integers, and it is not a root of any non-zero polynomial with rational coefficients. Consequently, its decimal expansion is infinite and non-repeating. While \(\pi\) can be defined exactly through mathematical identities, its value cannot be written down completely in finite decimal form. This makes numerical approximation both necessary and meaningful. \subsection{1.2 Why Numerical Approximation Is Important}\label{why-numerical-approximation-is-important} Although \(\pi\) has a precise mathematical definition, practical work in science, engineering, and computing usually requires a finite numerical representation. Calculations involving circles, waves, rotations, probabilities, or special functions often depend on a sufficiently accurate value of \(\pi\). The required precision varies widely: a few digits may be enough for everyday engineering estimates, while scientific simulations, cryptographic tests, or high-precision numerical experiments may require hundreds, thousands, or even millions of digits. Numerical approximation is also important because it connects theoretical mathematics with computation. Computing \(\pi\) provides a clear and well-understood problem for studying algorithms, floating-point arithmetic, arbitrary-precision arithmetic, convergence, and computational efficiency. It serves as a benchmark for numerical software and hardware, and it illustrates broader principles in numerical analysis, such as error control, stability, and the trade-off between accuracy and computational cost. \subsection{1.3 Scope of This Publication}\label{scope-of-this-publication} This publication examines the numerical calculation of \(\pi\) from both mathematical and computational perspectives. It begins with the historical development of approximation methods, then presents the mathematical foundations that make such calculations possible. The discussion proceeds to classical numerical methods, modern high-performance algorithms, practical implementation issues, convergence and error analysis, and performance optimization. Finally, it considers applications and extensions, showing how the study of \(\pi\) relates to wider problems in numerical computation and scientific computing. \section{2. Historical Background}\label{historical-background} \subsection{2.1 Ancient Geometric Approaches}\label{ancient-geometric-approaches} The earliest known approximations of \(\pi\) were obtained through geometric reasoning, especially by comparing the circumference or area of a circle with the perimeter or area of inscribed and circumscribed polygons. These methods were practical and intuitive, but their accuracy was limited by the number of polygon sides that could be constructed and computed by hand. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3000}} >{\raggedleft\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.4000}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3000}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Period / Culture \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Approximation \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Historical Significance \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Babylonian & \(25/8 = 3.125\) & One of the earliest recorded numerical approximations. \\ Egyptian, Rhind Papyrus & \((16/9)^2 = 256/81 \approx 3.1605\) & Derived from an area formula for a circle. \\ Greek, Archimedes & \(223/71 < \pi < 22/7\) & First rigorous bounding method using 96-gons. \\ \end{longtable} Archimedes' polygon method is especially important because it established a clear numerical strategy: increase the number of sides of a regular polygon to obtain tighter upper and lower bounds for \(\pi\). His result, \(\frac{223}{71} < \pi < \frac{22}{7},\) was a major achievement in ancient mathematics and remains a standard example of geometric approximation. \subsection{2.2 Medieval and Early Modern Refinements}\label{medieval-and-early-modern-refinements} After the Greek period, mathematicians in China, India, and the Islamic world continued to improve geometric approximations and developed more sophisticated computational techniques. In China, Liu Hui refined polygon-based methods and obtained approximations close to \(3.14159\). His student Zu Chongzhi later produced the famous rational approximation \(\frac{355}{113} \approx 3.14159292035,\) which is remarkably accurate for its size and remained one of the best rational approximations for centuries. In India, Aryabhata gave the approximation \(3.1416\), while Madhava of Sangamagrama developed early infinite series for trigonometric functions, including a series for \(\pi/4\). These Indian contributions were among the first to move beyond purely geometric methods toward analytical representations. In the Islamic world, Al-Kashi computed \(\pi\) to 16 decimal places in the 15th century using a polygon method with a very large number of sides. This work demonstrated that hand computation could achieve high precision, but it also highlighted the increasing difficulty of extending geometric methods further. \subsection{2.3 Early Analytical Series and Products}\label{early-analytical-series-and-products} The development of calculus and infinite series transformed the calculation of \(\pi\). Instead of relying on polygons, mathematicians began to express \(\pi\) through infinite sums, products, and arctangent identities. A central example is the Gregory-Leibniz series, \(\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots,\) which is simple but converges slowly. Its historical importance lies in showing that \(\pi\) could be represented analytically, not only geometrically. Other early analytical developments include: \begin{itemize} \tightlist \item Viète's infinite product for \(2/\pi\), \item Wallis's product for \(\pi/2\), \item Euler's transformations of series and products, \item Machin's formula, \end{itemize} \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right),\) which greatly accelerated numerical computation by using rapidly converging arctangent series. Machin-like formulas became especially important because they allowed mathematicians to compute many more digits of \(\pi\) than was practical with the Gregory-Leibniz series alone. This period marks the transition from classical geometric approximation to analytical numerical calculation. \subsection{2.4 Modern Computational Techniques}\label{modern-computational-techniques} The modern era of \(\pi\) calculation began with the development of electronic computers. Once high-precision arithmetic, fast multiplication, and efficient memory management became available, the calculation of \(\pi\) shifted from a mathematical exercise to a benchmark for computational performance. Several developments were particularly influential: \begin{itemize} \tightlist \item The Brent-Salamin algorithm, based on the arithmetic-geometric mean, provided quadratic convergence and became a standard method for high-precision computation. \item The Chudnovsky algorithm introduced a rapidly converging series that is especially efficient for large-scale calculations. \item The BBP formula enabled the extraction of individual hexadecimal digits of \(\pi\) without computing all preceding digits. \item Fast multiplication algorithms, such as those based on the fast Fourier transform, reduced the cost of handling very large numbers. \end{itemize} Modern \(\pi\) computations now involve trillions of digits and are used not only to study the constant itself, but also to test numerical algorithms, hardware reliability, parallelization strategies, and high-precision arithmetic. \subsection{2.5 Historical Lessons for Numerical Calculation}\label{historical-lessons-for-numerical-calculation} The historical development of \(\pi\) approximation shows a clear progression: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item \textbf{Geometric bounding} gave the first reliable numerical estimates. \item \textbf{Analytical series} made it possible to compute \(\pi\) to arbitrary precision in principle. \item \textbf{Modern algorithms} made high-precision computation practical and efficient. \end{enumerate} This progression is directly relevant to the methods examined in \textbf{Classical Numerical Methods} and \textbf{Modern Computational Algorithms}. It also motivates the later discussion of convergence, error bounds, and performance in \textbf{Convergence and Error Analysis} and \textbf{Performance and Optimization}. In short, the history of \(\pi\) calculation illustrates how mathematical insight and computational capability have evolved together. \section{3. Mathematical Foundations}\label{mathematical-foundations} \subsection{3.1 Geometric Definition and Circle Geometry}\label{geometric-definition-and-circle-geometry} The constant \(\pi\) is most fundamentally defined through the geometry of the circle. For a circle of radius \(r\), the circumference \(C\) and area \(A\) are given by \(C = 2\pi r, \qquad A = \pi r^2.\) Equivalently, \(\pi\) is the ratio of the circumference of a circle to its diameter: \(\pi = \frac{C}{d} = \frac{C}{2r}.\) For the unit circle, \(x^2 + y^2 = 1,\) the circumference is \(2\pi\), and the area is \(\pi\). In angular measure, a full rotation corresponds to \(2\pi\) radians, while a half-turn corresponds to \(\pi\) radians. Thus \(\pi\) also appears naturally as the half-period of the trigonometric functions. A classical geometric approach to approximating \(\pi\) is based on regular polygons inscribed in and circumscribed about a circle. For a regular \(n\)-gon inscribed in the unit circle, the side length is \(s_n = 2\sin\left(\frac{\pi}{n}\right),\) so its perimeter is \(P_n = 2n\sin\left(\frac{\pi}{n}\right).\) For a regular \(n\)-gon circumscribed about the unit circle, the side length is \(t_n = 2\tan\left(\frac{\pi}{n}\right),\) so its perimeter is \(Q_n = 2n\tan\left(\frac{\pi}{n}\right).\) Because the inscribed polygon lies inside the circle and the circumscribed polygon lies outside it, \(P_n < 2\pi < Q_n.\) Dividing by \(2\), one obtains the rigorous bounds \(n\sin\left(\frac{\pi}{n}\right) < \pi < n\tan\left(\frac{\pi}{n}\right).\) These inequalities form the mathematical basis of the polygon method discussed in \textbf{4. Classical Numerical Methods}. Historically, this approach was central to the work of Archimedes and is reviewed in \textbf{2. Historical Background}. Although geometrically intuitive, polygon-based approximations converge relatively slowly: the error decreases only algebraically with \(n\), which makes this method less suitable for computing very high-precision values of \(\pi\) compared with later analytical methods. \subsection{3.2 Trigonometric Identities and the Arctangent}\label{trigonometric-identities-and-the-arctangent} Trigonometric functions provide a direct bridge between circle geometry and analysis. The unit circle parametrization \((\cos \theta, \sin \theta)\) shows that \(\pi\) is the angle corresponding to a half-turn. The fundamental identity \(\sin^2 \theta + \cos^2 \theta = 1\) encodes the equation of the unit circle, while the periodicity \(\sin(\theta + 2\pi) = \sin \theta, \qquad \cos(\theta + 2\pi) = \cos \theta\) shows that \(2\pi\) is the natural period of circular motion. A particularly important identity is \(\tan\left(\frac{\pi}{4}\right) = 1.\) Since the arctangent function is the inverse of the tangent function on the principal branch, \(\arctan(1) = \frac{\pi}{4}.\) Therefore, \(\pi = 4\arctan(1).\) This identity is one of the most important starting points for numerical computation of \(\pi\), because it reduces the problem of approximating \(\pi\) to the problem of approximating an inverse trigonometric function. The addition formula for the tangent function, \(\tan(a+b) = \frac{\tan a + \tan b}{1 - \tan a \tan b},\) leads to the arctangent addition identity \(\arctan x + \arctan y = \arctan\left(\frac{x+y}{1-xy}\right),\) up to the usual adjustment by integer multiples of \(\pi\) when the argument crosses branch boundaries. This identity allows one to rewrite \(\pi/4\) as a combination of arctangents of smaller arguments. For example, \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right).\) Formulas of this type are known as Machin-like formulas. They are mathematically significant because replacing \(\arctan(1)\) by arctangents of smaller numbers greatly improves the convergence of the associated power series. Machin-like formulas are developed further in \textbf{5. Modern Computational Algorithms}. Other trigonometric identities are also important. The double-angle formulas \(\sin(2\theta) = 2\sin\theta\cos\theta,\) \(\cos(2\theta) = \cos^2\theta - \sin^2\theta,\) and the half-angle formulas \(\sin\left(\frac{\theta}{2}\right) = \sqrt{\frac{1-\cos\theta}{2}},\) \(\cos\left(\frac{\theta}{2}\right) = \sqrt{\frac{1+\cos\theta}{2}},\) are useful both in geometric polygon constructions and in iterative numerical algorithms. They allow one to compute trigonometric values for successively smaller angles, a technique that appears in both historical and modern methods for approximating \(\pi\). \subsection{3.3 Infinite Series and Products}\label{infinite-series-and-products} The transition from geometric approximations to analytical computation of \(\pi\) is made possible by infinite series and infinite products. These representations express \(\pi\) as the limit of a sequence of rational or algebraic quantities, which is essential for numerical calculation. The Taylor series for the sine and cosine functions are \(\sin x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!},\) \(\cos x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n}}{(2n)!}.\) These series converge for all real \(x\), and they provide a rigorous analytic definition of the trigonometric functions. They also show that \(\pi\) is deeply connected to the factorial function and to the structure of analytic functions. A more directly useful series for computing \(\pi\) is the Taylor series for the arctangent: \(\arctan x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{2n+1}, \qquad |x| \le 1.\) Setting \(x=1\) gives the Gregory-Leibniz series, \(\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \cdots.\) Equivalently, \(\pi = 4\left( 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \cdots \right).\) This series is historically important and mathematically simple, but it converges slowly. If \(S_N = \sum_{n=0}^{N} (-1)^n \frac{1}{2n+1},\) then the alternating-series remainder satisfies \(\left| \frac{\pi}{4} - S_N \right| \le \frac{1}{2N+3}.\) Thus, to obtain \(d\) correct decimal digits using the Gregory-Leibniz series, roughly \(10^d\) terms are required. This slow convergence is one of the main reasons why later methods use arctangents of smaller arguments or more rapidly convergent series. Machin-like formulas exploit the arctangent series by using identities such as \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right).\) Because the arguments \(1/5\) and \(1/239\) are much smaller than \(1\), the corresponding series converge much faster. This principle is central to many efficient classical and modern algorithms for computing \(\pi\). Infinite products also provide representations of \(\pi\). The Wallis product is \(\frac{\pi}{2} = \prod_{n=1}^{\infty} \frac{(2n)^2}{(2n-1)(2n+1)}.\) The Viète product is \(\frac{2}{\pi} = \prod_{n=2}^{\infty} \cos\left(\frac{\pi}{2^n}\right).\) These products are historically significant and are closely related to trigonometric identities and polygonal approximations. However, for high-precision numerical computation, infinite products are generally less convenient than rapidly convergent series or iterative algorithms. \subsection{3.4 Integral Representations}\label{integral-representations} Integral representations of \(\pi\) connect the constant to calculus, probability, and special functions. They are especially useful because many numerical methods approximate definite integrals by quadrature, series expansion, or transformation into other integral forms. A basic integral representation follows from the arctangent: \(\int_0^1 \frac{dx}{1+x^2} = \arctan(1) = \frac{\pi}{4}.\) Therefore, \(\pi = 4\int_0^1 \frac{dx}{1+x^2}.\) This integral is closely related to the Gregory-Leibniz series, since \(\frac{1}{1+x^2} = \sum_{n=0}^{\infty} (-1)^n x^{2n}, \qquad |x|<1,\) and termwise integration gives \(\int_0^1 \frac{dx}{1+x^2} = \sum_{n=0}^{\infty} (-1)^n \frac{1}{2n+1}.\) Another important integral is \(\pi = \int_{-1}^{1} \frac{dx}{\sqrt{1-x^2}}.\) This follows from the antiderivative \(\int \frac{dx}{\sqrt{1-x^2}} = \arcsin x.\) A geometrically direct integral representation is the area of the unit circle: \(\pi = \int_{-1}^{1} 2\sqrt{1-x^2}\,dx.\) This expresses \(\pi\) as the area under the upper and lower semicircles of the unit circle. A related improper integral is \(\pi = 2\int_0^{\infty} \frac{dx}{1+x^2}.\) This representation is useful in complex analysis and in the study of rational functions. A particularly important integral representation is the Gaussian integral: \(\int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi}.\) Squaring both sides gives \(\pi = \left( \int_{-\infty}^{\infty} e^{-x^2}\,dx \right)^2.\) This identity connects \(\pi\) to probability theory, Fourier analysis, and numerical quadrature. It also illustrates that \(\pi\) is not merely a geometric constant but a central object in analysis. Finally, \(\pi\) can be expressed using the beta and gamma functions: \(B\left(\frac{1}{2},\frac{1}{2}\right) = \int_0^1 t^{-1/2}(1-t)^{-1/2}\,dt = \pi.\) Since \(B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)},\) and \(\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi},\) one obtains \(B\left(\frac{1}{2},\frac{1}{2}\right) = \frac{\Gamma(1/2)\Gamma(1/2)}{\Gamma(1)} = \pi.\) These special-function representations are less commonly used for elementary computation of \(\pi\), but they are important in advanced numerical analysis and in the derivation of rapidly convergent series. \subsection{3.5 Role of the Mathematical Foundations in Numerical Calculation}\label{role-of-the-mathematical-foundations-in-numerical-calculation} The mathematical foundations of \(\pi\) calculation can be summarized as follows. Circle geometry provides the original definition of \(\pi\) and yields rigorous bounds through polygonal approximations. Trigonometric identities connect \(\pi\) to inverse trigonometric functions, especially the arctangent. Infinite series and products convert these identities into computable limits. Integral representations provide additional analytical forms that can be approximated by quadrature or transformed into series. Each representation has different computational properties. Geometric polygon methods are simple but converge slowly. The Gregory-Leibniz series is easy to state but inefficient for high precision. Machin-like arctangent formulas improve convergence by using smaller arguments. Integral representations offer flexibility but require careful numerical treatment. Rapidly convergent series and iterative algorithms, discussed in \textbf{5. Modern Computational Algorithms}, build on these foundations to achieve high-precision values of \(\pi\) with far fewer operations. Because \(\pi\) is irrational and transcendental, as noted in \textbf{1. Introduction}, no finite decimal representation is exact. Therefore, every numerical method produces an approximation whose accuracy depends on the number of terms, iterations, or quadrature points used. The convergence behavior, error bounds, and numerical stability of these methods are analyzed in \textbf{7. Convergence and Error Analysis}. Their practical realization in software is treated in \textbf{6. Algorithm Implementation}, while their computational cost and optimization are examined in \textbf{8. Performance and Optimization}. \section{4. Classical Numerical Methods}\label{classical-numerical-methods} \subsection{4.1 Archimedes' Polygon Method}\label{archimedes-polygon-method} The oldest systematic numerical method for approximating \(\pi\) is geometric. It is based on comparing the circumference of a circle with the perimeters of regular polygons inscribed in and circumscribed about the circle. This method is historically important because it gives rigorous upper and lower bounds for \(\pi\) without using infinite series or limits. Consider a circle of radius \(1\). Its circumference is \(2\pi\). If a regular \(n\)-gon is inscribed in the circle, its perimeter is smaller than the circumference. If a regular \(n\)-gon is circumscribed about the circle, its perimeter is larger. For the unit circle, the side length of an inscribed regular \(n\)-gon is \(2\sin\left(\frac{\pi}{n}\right),\) and the side length of a circumscribed regular \(n\)-gon is \(2\tan\left(\frac{\pi}{n}\right).\) Therefore, the perimeters are \(P_{\text{in}}(n)=2n\sin\left(\frac{\pi}{n}\right)\) and \(P_{\text{out}}(n)=2n\tan\left(\frac{\pi}{n}\right).\) Since \(P_{\text{in}}(n) < 2\pi < P_{\text{out}}(n),\) we obtain the classical bounds \(n\sin\left(\frac{\pi}{n}\right) < \pi < n\tan\left(\frac{\pi}{n}\right).\) These inequalities provide a lower and an upper approximation to \(\pi\) for every integer \(n \ge 3\). Archimedes used this idea by starting with a hexagon and repeatedly doubling the number of sides until he reached a 96-gon. Using only elementary geometry and careful rational approximations to square roots, he proved the famous bounds \(\frac{223}{71} < \pi < \frac{22}{7}.\) In decimal form, \(3.140845\ldots < \pi < 3.142857\ldots\) This was a remarkable achievement for its time and remained one of the best known approximations for many centuries. The polygon method can also be written in a recursive form that avoids direct use of trigonometric functions. Let \(a_n = 2\sin\left(\frac{\pi}{n}\right)\) be the side length of the inscribed regular \(n\)-gon, and let \(b_n = 2\tan\left(\frac{\pi}{n}\right)\) be the side length of the circumscribed regular \(n\)-gon. Starting from a hexagon, \(a_6 = 1, \qquad b_6 = \frac{2}{\sqrt{3}}.\) When the number of sides is doubled, the side lengths satisfy \(a_{2n} = \sqrt{2-\sqrt{4-a_n^2}},\) and \(b_{2n} = \frac{b_n}{1+\sqrt{1+\left(\frac{b_n}{2}\right)^2}}.\) The corresponding bounds for \(\pi\) are then \(\frac{n a_n}{2} < \pi < \frac{n b_n}{2}.\) The following table shows how the bounds improve as the number of sides increases. \begin{longtable}[]{@{}rrr@{}} \toprule\noalign{} \(n\) & Lower bound \(n\sin(\pi/n)\) & Upper bound \(n\tan(\pi/n)\) \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot 6 & 3.000000 & 3.464102 \\ 12 & 3.105829 & 3.215390 \\ 24 & 3.132629 & 3.159660 \\ 48 & 3.139350 & 3.143332 \\ 96 & 3.141032 & 3.142714 \\ \end{longtable} The convergence of the polygon method is relatively slow compared with modern algorithms, but it is much faster than the simplest arctangent series. Using the Taylor expansions \(\sin x = x - \frac{x^3}{6} + \frac{x^5}{120} - \cdots\) and \(\tan x = x + \frac{x^3}{3} + \frac{x^5}{15} + \cdots,\) with \(x=\pi/n\), we obtain \(n\sin\left(\frac{\pi}{n}\right) = \pi - \frac{\pi^3}{6n^2} + O\left(\frac{1}{n^4}\right),\) and \(n\tan\left(\frac{\pi}{n}\right) = \pi + \frac{\pi^3}{3n^2} + O\left(\frac{1}{n^4}\right).\) Thus the lower and upper errors are both of order \(O\left(\frac{1}{n^2}\right).\) If the number of sides is doubled at each step, the error is reduced by approximately a factor of \(4\). In terms of the number of decimal digits \(d\), the required number of sides grows roughly like \(n \sim 10^{d/2}.\) Equivalently, if one doubles the number of sides repeatedly, the number of doublings needed to gain \(d\) decimal digits is proportional to \(d\). This is still far less efficient than the fastest modern algorithms, but it is a clear and rigorous classical method. \subsection{4.2 The Leibniz Series}\label{the-leibniz-series} The development of calculus shifted the computation of \(\pi\) from geometry to analysis. One of the earliest and simplest analytical methods is the Leibniz series, which comes from the Taylor series for the arctangent function. For \(|x|<1\), the arctangent function has the power series expansion \(\arctan x = \sum_{k=0}^{\infty} (-1)^k \frac{x^{2k+1}}{2k+1}.\) At \(x=1\), the series converges conditionally by the alternating series test: \(\arctan(1) = 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots.\) Since \(\arctan(1)=\frac{\pi}{4},\) we obtain \(\frac{\pi}{4} = 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots,\) or equivalently, \(\pi = 4\left( 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots \right).\) This is the Leibniz series for \(\pi\). Let \(S_N = \sum_{k=0}^{N} (-1)^k \frac{1}{2k+1}.\) Then \(\pi \approx 4S_N.\) Because the terms decrease in magnitude to zero, the alternating series error estimate gives \(\left| \frac{\pi}{4}-S_N \right| \le \frac{1}{2N+3}.\) Multiplying by \(4\), we obtain the error bound \(|\pi-4S_N| \le \frac{4}{2N+3}.\) The partial sums also alternate around the true value. If \(N\) is even, then \(4S_N\) is an overestimate of \(\pi\). If \(N\) is odd, then \(4S_N\) is an underestimate. Thus the sequence of partial sums provides alternating upper and lower bounds for \(\pi\). The convergence of the Leibniz series is slow. The error is of order \(O\left(\frac{1}{N}\right).\) This means that gaining one additional decimal digit generally requires about ten times as many terms. For example: \begin{longtable}[]{@{} >{\raggedleft\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedleft\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedleft\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedleft Number of terms \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Error bound \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Approximate accuracy \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot 100 & \(4/201 \approx 0.0199\) & about 1 decimal digit \\ 10,000 & \(4/19999 \approx 0.000200\) & about 3 decimal digits \\ 1,000,000 & \(4/1999999 \approx 2.0\times 10^{-6}\) & about 5 to 6 decimal digits \\ \end{longtable} To obtain \(d\) correctly rounded decimal digits, the alternating-series bound suggests that roughly \(N \gtrsim 4\times 10^d\) terms are needed. For example, computing ten decimal digits of \(\pi\) using the Leibniz series would require on the order of tens of billions of terms. This makes the method impractical for high-precision computation, although it is extremely useful for illustrating the basic ideas of series convergence and error estimation. \subsection{4.3 The Gregory-Leibniz Formula}\label{the-gregory-leibniz-formula} The Gregory-Leibniz formula is the specific identity that connects the arctangent series to \(\pi\). Historically, James Gregory derived the arctangent series, and Gottfried Wilhelm Leibniz evaluated it at \(x=1\) to obtain a formula for \(\pi\). The Gregory-Leibniz formula is \(\frac{\pi}{4} = \arctan(1) = \sum_{k=0}^{\infty} (-1)^k \frac{1}{2k+1}.\) Equivalently, \(\pi = 4\arctan(1) = 4\sum_{k=0}^{\infty} (-1)^k \frac{1}{2k+1}.\) Although the names ``Leibniz series'' and ``Gregory-Leibniz formula'' are often used interchangeably, the distinction made here is that the Leibniz series refers to the general arctangent series, while the Gregory-Leibniz formula refers to the particular identity used to compute \(\pi\). The formula can be derived by integrating the geometric series \(\frac{1}{1+t^2} = 1-t^2+t^4-t^6+\cdots\) from \(0\) to \(1\). Since \(\int_0^1 \frac{dt}{1+t^2} = \arctan(1) = \frac{\pi}{4},\) term-by-term integration gives \(\frac{\pi}{4} = 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots.\) The convergence properties are the same as for the Leibniz series. The series is conditionally convergent at \(x=1\), not absolutely convergent. The error after summing through the term \(1/(2N+1)\) satisfies \(|\pi-4S_N| \le \frac{4}{2N+3}.\) The slow convergence is caused by the fact that the argument \(x=1\) lies at the boundary of the interval of convergence of the arctangent Taylor series. For \(|x|<1\), the general arctangent series \(\arctan x = x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\cdots\) converges faster because the terms contain powers of \(x^{2k+1}\). For \(0{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Method \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Basic approximation \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Error estimate \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Convergence behavior \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Archimedes' polygon method & \(n\sin(\pi/n) < \pi < n\tan(\pi/n)\) & Lower error \(\sim \pi^3/(6n^2)\), upper error \(\sim \pi^3/(3n^2)\) & \(O(n^{-2})\) in the number of sides; error reduced by about a factor of \(4\) when \(n\) is doubled \\ Leibniz series / Gregory-Leibniz formula & \(\pi \approx 4\sum_{k=0}^{N}(-1)^k/(2k+1)\) & \(|\pi-4S_N| \le 4/(2N+3)\) & \(O(N^{-1})\) in the number of terms; roughly ten times more terms are needed for one additional decimal digit \\ \end{longtable} The polygon method has the advantage of producing rigorous upper and lower bounds. It is also relatively efficient among classical methods because the error decreases quadratically with the number of sides. However, it requires increasingly accurate square roots as the number of sides grows, and naive implementations can suffer from numerical cancellation when the side lengths become very small. The Leibniz and Gregory-Leibniz series are simpler to state and implement, requiring only addition, subtraction, and division. They also provide alternating bounds for \(\pi\). Their major limitation is slow convergence. Because the error decreases only linearly with the number of terms, the method becomes impractical for high-precision computation. These classical methods are therefore best viewed as foundational tools. They illustrate how geometric and analytical ideas lead to numerical approximations of \(\pi\), and they provide simple error bounds that are useful for later discussion in \textbf{7. Convergence and Error Analysis}. Their limitations also motivate the faster arctangent-based formulas, iterative methods, and high-precision algorithms considered in \textbf{5. Modern Computational Algorithms}. \section{5. Modern Computational Algorithms}\label{modern-computational-algorithms} \subsection{5.1 Machin-like Formulas}\label{machin-like-formulas} A major step beyond the slowly converging series discussed in \textbf{4. Classical Numerical Methods} is the use of Machin-like formulas. These identities exploit the arctangent addition formula, introduced in \textbf{3. Mathematical Foundations}, to express \(\pi\) as a linear combination of arctangent series with small arguments. Because the Taylor series for \(\arctan x\) converges much faster when \(|x|\) is small, Machin-like formulas are far more efficient than the Gregory-Leibniz series. The arctangent addition formula is \(\arctan x + \arctan y = \arctan\left(\frac{x+y}{1-xy}\right)\) up to an additive multiple of \(\pi/2\), depending on the quadrant. By choosing suitable rational values of \(x\) and \(y\), one obtains identities of the form \(\frac{\pi}{4} = \sum_{j=1}^{m} c_j \arctan\left(\frac{1}{b_j}\right),\) where \(c_j\) are integers and \(b_j > 1\). The classical example, due to John Machin, is \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right).\) Using the Taylor expansion \(\arctan z = \sum_{k=0}^{\infty} (-1)^k \frac{z^{2k+1}}{2k+1}, \qquad |z|\le 1,\) Machin's formula becomes \(\frac{\pi}{4} = 4\sum_{k=0}^{\infty} (-1)^k \frac{1}{(2k+1)5^{2k+1}} - \sum_{k=0}^{\infty} (-1)^k \frac{1}{(2k+1)239^{2k+1}}.\) The first series dominates the convergence because \(1/5\) is larger than \(1/239\). Each additional term in the \(1/5\) series reduces the error by a factor of approximately \(25\), giving roughly \(\log_{10}(25) \approx 1.4\) decimal digits per term. This is a dramatic improvement over the Gregory-Leibniz series, where the number of terms grows exponentially with the number of desired decimal digits. More generally, if \(\frac{\pi}{4} = \sum_{j=1}^{m} c_j \arctan\left(\frac{1}{b_j}\right),\) and each arctangent series is truncated after \(N\) terms, then the truncation error satisfies \(\left| \frac{\pi}{4} - \sum_{j=1}^{m} c_j \sum_{k=0}^{N} (-1)^k \frac{1}{(2k+1)b_j^{2k+1}} \right| \le \sum_{j=1}^{m} \frac{|c_j|}{(2N+3)b_j^{2N+3}}.\) This bound is useful for estimating how many terms are required to reach a target precision. In implementation, Machin-like formulas are attractive because they can be evaluated using rational arithmetic or fixed-point arithmetic with guard digits. A common optimization is to update terms recursively rather than recomputing powers. For a fixed \(b\), define \(u_k = (-1)^k \frac{1}{(2k+1)b^{2k+1}}.\) Then \(u_{k+1} = -u_k \frac{2k+1}{2k+3} \frac{1}{b^2}.\) This avoids repeated exponentiation and makes the method efficient for moderate precision calculations. Machin-like formulas were historically important because they made high-precision computation of \(\pi\) feasible before electronic computers. They remain useful today for educational implementations, small-scale calculations, and as a baseline for comparing more advanced algorithms. However, for computations involving millions or billions of digits, their linear-in-terms convergence is usually outperformed by algorithms such as Chudnovsky's. \subsection{5.2 The Chudnovsky Algorithm}\label{the-chudnovsky-algorithm} One of the most important modern algorithms for high-precision computation of \(\pi\) is the Chudnovsky algorithm, developed by the Chudnovsky brothers in 1988. It is based on a rapidly convergent Ramanujan-type series and is one of the principal methods used in record-setting \(\pi\) calculations. The Chudnovsky formula is \(\frac{1}{\pi} = 12 \sum_{k=0}^{\infty} (-1)^k \frac{(6k)! \left(545140134k + 13591409\right)} {(3k)! (k!)^3 640320^{3k+3/2}}.\) Equivalently, \(\pi = \left[ 12 \sum_{k=0}^{\infty} (-1)^k \frac{(6k)! \left(545140134k + 13591409\right)} {(3k)! (k!)^3 640320^{3k+3/2}} \right]^{-1}.\) The series converges extremely quickly. Successive terms decrease by a factor of roughly \(10^{14.15}\), meaning that each term contributes approximately 14 decimal digits of accuracy. Thus, to compute \(d\) decimal digits of \(\pi\), one needs on the order of \(\frac{d}{14.15}\) terms, rather than the much larger number required by Machin-like formulas. A practical implementation avoids computing factorials directly. Instead, the terms are updated recursively. Let \(A_k = \frac{(6k)! \left(545140134k + 13591409\right)} {(3k)! (k!)^3 640320^{3k}}.\) Then the Chudnovsky sum can be written as \(\frac{1}{\pi} = \frac{12}{640320^{3/2}} \sum_{k=0}^{\infty} (-1)^k A_k.\) The terms \(A_k\) satisfy a rational recurrence relation of the form \(A_{k+1} = -A_k \frac{(6k+1)(6k+2)(6k+3)(6k+4)(6k+5)(6k+6)} {(3k+1)(3k+2)(3k+3)(k+1)^3 640320^3} \cdot \frac{545140134(k+1)+13591409} {545140134k+13591409}.\) This recurrence allows the summation to be performed using only multiplications, divisions, and additions of high-precision numbers. For very high precision, the Chudnovsky algorithm is usually combined with \textbf{binary splitting}. Binary splitting is a divide-and-conquer technique for summing a sequence of rational terms. Instead of adding terms one by one, the summation range is split into two halves, each half is evaluated recursively, and the partial results are combined. If a block of terms is represented as a rational number \(N/D\), then two blocks \(N_1/D_1\) and \(N_2/D_2\) can be combined as \(\frac{N_1}{D_1} + \frac{N_2}{D_2} = \frac{N_1 D_2 + N_2 D_1}{D_1 D_2}.\) This approach keeps intermediate numbers well balanced and allows the use of fast multiplication algorithms, such as those based on the fast Fourier transform. For record-scale computations, the cost of multiplication dominates the total running time, so the efficiency of the multiplication algorithm is crucial. The Chudnovsky algorithm is especially effective for large-scale computations because: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item It requires relatively few terms for a given number of digits. \item The terms can be updated recursively. \item Binary splitting reduces the cost of summation. \item Fast multiplication can be applied to the large integers that arise. \end{enumerate} Because of these properties, the Chudnovsky algorithm is one of the standard choices for computing \(\pi\) to millions, billions, or more decimal digits. \subsection{5.3 Iterative Methods: Brent-Salamin and the AGM}\label{iterative-methods-brent-salamin-and-the-agm} Another important class of modern algorithms for computing \(\pi\) is based on iterative methods, particularly the arithmetic-geometric mean, or AGM. The most widely known AGM-based method is the Brent-Salamin algorithm, also called the Gauss-Legendre algorithm. The AGM of two positive numbers \(a\) and \(b\) is the common limit of the sequences \(a_{n+1} = \frac{a_n + b_n}{2},\) \(b_{n+1} = \sqrt{a_n b_n}.\) Both sequences converge rapidly to the same value. The Brent-Salamin algorithm uses this convergence to compute \(\pi\). The algorithm is initialized with \(a_0 = 1, \qquad b_0 = \frac{1}{\sqrt{2}}, \qquad t_0 = \frac{1}{4}, \qquad p_0 = 1.\) The iteration is then \(a_{n+1} = \frac{a_n + b_n}{2},\) \(b_{n+1} = \sqrt{a_n b_n},\) \(t_{n+1} = t_n - p_n (a_n - a_{n+1})^2,\) \(p_{n+1} = 2p_n.\) After \(n\) iterations, \(\pi\) is approximated by \(\pi \approx \frac{(a_n + b_n)^2}{4t_n}.\) The convergence is quadratic: the number of correct digits roughly doubles with each iteration. For example, if an iteration produces 10 correct digits, the next iteration may produce about 20, the next about 40, and so on. This makes the Brent-Salamin algorithm very efficient for moderate to high precision. The main computational cost in each iteration is the high-precision square root. The square root can itself be computed using Newton's method: \(x_{k+1} = \frac{1}{2} \left( x_k + \frac{y}{x_k} \right),\) which converges quadratically to \(\sqrt{y}\). Thus, the Brent-Salamin algorithm is a nested iterative method: the outer AGM iteration converges quadratically, and the inner square-root iteration also converges quadratically. The Brent-Salamin algorithm has several practical advantages: \begin{itemize} \tightlist \item It is mathematically simple. \item It is numerically stable. \item It requires only addition, multiplication, and square roots. \item It is well suited to arbitrary-precision arithmetic. \item It is often easier to implement correctly than the Chudnovsky algorithm. \end{itemize} However, for extremely large precision, the Chudnovsky algorithm is often faster in practice because it can exploit fast multiplication more effectively and requires fewer high-precision operations per digit. The Brent-Salamin algorithm remains an excellent choice for many applications, especially when implementation simplicity and reliability are important. Other iterative methods, such as Borwein-type algorithms, also exist and can have higher-order convergence. Nevertheless, the Brent-Salamin algorithm is the most commonly cited iterative method for high-precision computation of \(\pi\). \subsection{5.4 Practical Comparison and Algorithm Selection}\label{practical-comparison-and-algorithm-selection} The choice of algorithm depends on the required precision, the available arithmetic, and the implementation environment. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Algorithm \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Convergence behavior \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Main strengths \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Main limitations \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Machin-like formulas & Exponential in the number of terms; roughly 1.4 digits per term for Machin's formula & Simple, rational arithmetic, easy to implement & Requires many terms for very high precision \\ Chudnovsky algorithm & Approximately 14 decimal digits per term & Very fast for large-scale computations; compatible with binary splitting and fast multiplication & More complex; requires careful high-precision integer arithmetic \\ Brent-Salamin / AGM & Quadratic convergence; digits roughly double each iteration & Simple, stable, efficient for moderate precision & Square-root cost; usually less efficient than Chudnovsky at extreme precision \\ \end{longtable} For a few decimal digits, any of these methods is more than sufficient. For educational purposes, Machin-like formulas are often the clearest because they connect directly to the arctangent series and require only basic arithmetic. For moderate precision, the Brent-Salamin algorithm is attractive because of its quadratic convergence and relative simplicity. For record-scale computations, the Chudnovsky algorithm is typically preferred, especially when combined with binary splitting and fast multiplication. In all cases, high-precision computation requires careful handling of rounding, guard digits, and error accumulation. The convergence rates described here give the ideal mathematical behavior, but practical performance also depends on the efficiency of the underlying arithmetic operations. These issues are examined more formally in \textbf{7. Convergence and Error Analysis}, while computational cost, memory usage, and optimization strategies are discussed in \textbf{8. Performance and Optimization}. \section{6. Algorithm Implementation}\label{algorithm-implementation} \subsection{6.1 Choosing a Numeric Representation}\label{choosing-a-numeric-representation} The first implementation decision is the numeric type used to represent intermediate and final values. The appropriate choice depends on the target precision, the algorithm, and the performance requirements. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Target precision \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Typical numeric type \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Practical notes \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot A few digits to about 15 decimal digits & IEEE 754 \texttt{double} or \texttt{float64} & Simple and fast, but limited to roughly 15-16 correct decimal digits. Suitable for demonstrations and low-precision applications. \\ About 16 to 30 decimal digits & \texttt{long\ double}, extended precision, or arbitrary-precision decimal & Useful when slightly more than double precision is needed, but support and behavior vary by platform. \\ Hundreds to millions of digits & Arbitrary-precision integers, rationals, or multiprecision floating-point types & Required for high-precision computation. Libraries such as MPFR, GMP, mpmath, Boost.Multiprecision, or Java \texttt{BigInteger}/\texttt{BigDecimal} are commonly used. \\ \end{longtable} For high-precision work, it is usually better to use arbitrary-precision arithmetic than to rely on native floating-point types. In many implementations, the algorithm is written in terms of a precision parameter, for example: \begin{Shaded} \begin{Highlighting}[] \NormalTok{mp.dps }\OperatorTok{=}\NormalTok{ target\_digits }\OperatorTok{+}\NormalTok{ guard\_digits} \end{Highlighting} \end{Shaded} where \texttt{target\_digits} is the number of decimal digits required and \texttt{guard\_digits} is an additional safety margin. Guard digits protect against accumulated rounding error. A common choice is 10-20 extra decimal digits for moderate precision, and more for very large computations. For series-based methods, integer or rational arithmetic can be preferable to floating-point arithmetic. Instead of computing each term as a floating-point number, one can compute scaled integer terms and perform the final division only once. This reduces rounding error and can improve performance when fast integer multiplication is available. \subsection{6.2 Implementing Classical Methods}\label{implementing-classical-methods} Classical methods are useful for teaching and for low-precision computation. As discussed in \textbf{4. Classical Numerical Methods}, they are mathematically simple but generally converge more slowly than modern algorithms. \subsubsection{6.2.1 Archimedes' Polygon Method}\label{archimedes-polygon-method-1} The polygon method approximates \(\pi\) by comparing the perimeters of inscribed and circumscribed regular polygons. For a unit circle, \(n\sin\left(\frac{\pi}{n}\right) < \pi < n\tan\left(\frac{\pi}{n}\right).\) A direct implementation can evaluate the sine and tangent at high precision: \begin{Shaded} \begin{Highlighting}[] \ImportTok{import}\NormalTok{ mpmath }\ImportTok{as}\NormalTok{ mp} \KeywordTok{def}\NormalTok{ pi\_polygon(n, digits):} \NormalTok{ mp.dps }\OperatorTok{=}\NormalTok{ digits }\OperatorTok{+} \DecValTok{10} \NormalTok{ lower }\OperatorTok{=}\NormalTok{ n }\OperatorTok{*}\NormalTok{ mp.sin(mp.pi }\OperatorTok{/}\NormalTok{ n)} \NormalTok{ upper }\OperatorTok{=}\NormalTok{ n }\OperatorTok{*}\NormalTok{ mp.tan(mp.pi }\OperatorTok{/}\NormalTok{ n)} \ControlFlowTok{return}\NormalTok{ (lower }\OperatorTok{+}\NormalTok{ upper) }\OperatorTok{/} \DecValTok{2} \end{Highlighting} \end{Shaded} However, this implementation has a circularity problem: it uses \texttt{mp.pi} to compute the angle. A production implementation should either use a known initial approximation of \(\pi\) or use a recurrence for polygon side lengths that avoids explicit trigonometric evaluation. A more robust approach is to start from a hexagon and repeatedly double the number of sides using stable recurrences for the side lengths. This avoids repeated calls to high-precision trigonometric functions. The polygon method is rarely used for high-precision computation because its convergence is only \(O(n^{-2})\), where \(n\) is the number of sides. \subsubsection{6.2.2 Gregory-Leibniz Series}\label{gregory-leibniz-series} The Gregory-Leibniz series is \(\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots\) or \(\pi = 4\sum_{k=0}^{\infty} \frac{(-1)^k}{2k+1}.\) Because the series is alternating and the terms decrease monotonically, the error after adding the \(k\)-th term is bounded by the magnitude of the next term. Therefore, a safe stopping condition is to stop when the current term is smaller than the desired tolerance divided by 4. \begin{Shaded} \begin{Highlighting}[] \ImportTok{import}\NormalTok{ mpmath }\ImportTok{as}\NormalTok{ mp} \KeywordTok{def}\NormalTok{ pi\_leibniz(digits):} \NormalTok{ mp.dps }\OperatorTok{=}\NormalTok{ digits }\OperatorTok{+} \DecValTok{10} \NormalTok{ tol }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{10}\NormalTok{) }\OperatorTok{**}\NormalTok{ (}\OperatorTok{{-}}\NormalTok{digits) }\OperatorTok{/} \DecValTok{4} \NormalTok{ s }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{0}\NormalTok{)} \NormalTok{ sign }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{1}\NormalTok{)} \NormalTok{ k }\OperatorTok{=} \DecValTok{0} \ControlFlowTok{while} \VariableTok{True}\NormalTok{:} \NormalTok{ term }\OperatorTok{=}\NormalTok{ sign }\OperatorTok{/}\NormalTok{ (}\DecValTok{2}\OperatorTok{*}\NormalTok{k }\OperatorTok{+} \DecValTok{1}\NormalTok{)} \NormalTok{ s }\OperatorTok{+=}\NormalTok{ term} \ControlFlowTok{if} \BuiltInTok{abs}\NormalTok{(term) }\OperatorTok{\textless{}}\NormalTok{ tol:} \ControlFlowTok{break} \NormalTok{ sign }\OperatorTok{=} \OperatorTok{{-}}\NormalTok{sign} \NormalTok{ k }\OperatorTok{+=} \DecValTok{1} \ControlFlowTok{return} \DecValTok{4} \OperatorTok{*}\NormalTok{ s} \end{Highlighting} \end{Shaded} This implementation is simple and numerically stable, but it is inefficient for high precision. Since the convergence is only \(O(N^{-1})\), obtaining \(p\) decimal digits requires on the order of \(10^p\) terms. For example, 100 correct digits would require an impractically large number of iterations. \subsection{6.3 Implementing Machin-like Formulas}\label{implementing-machin-like-formulas} Machin-like formulas improve on the Gregory-Leibniz series by using arctangent arguments smaller than 1. A standard identity is \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right).\) Thus, \(\pi = 16\arctan\left(\frac{1}{5}\right) - 4\arctan\left(\frac{1}{239}\right).\) The arctangent Taylor series is \(\arctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \cdots = \sum_{k=0}^{\infty} \frac{(-1)^k x^{2k+1}}{2k+1}, \qquad |x| \le 1.\) For high-precision computation, it is usually better to implement this series directly rather than calling a library \texttt{atan} function, because the series gives explicit control over precision and stopping criteria. A stable implementation uses a recurrence for the terms: \(t_k = \frac{(-1)^k x^{2k+1}}{2k+1},\) \(t_{k+1} = -t_k x^2 \frac{2k+1}{2k+3}.\) \begin{Shaded} \begin{Highlighting}[] \ImportTok{import}\NormalTok{ mpmath }\ImportTok{as}\NormalTok{ mp} \KeywordTok{def}\NormalTok{ atan\_series(x, tol):} \NormalTok{ s }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{0}\NormalTok{)} \NormalTok{ term }\OperatorTok{=}\NormalTok{ x} \NormalTok{ k }\OperatorTok{=} \DecValTok{0} \ControlFlowTok{while} \BuiltInTok{abs}\NormalTok{(term) }\OperatorTok{\textgreater{}=}\NormalTok{ tol:} \NormalTok{ s }\OperatorTok{+=}\NormalTok{ term} \NormalTok{ k }\OperatorTok{+=} \DecValTok{1} \NormalTok{ term }\OperatorTok{=} \OperatorTok{{-}}\NormalTok{term }\OperatorTok{*}\NormalTok{ x }\OperatorTok{*}\NormalTok{ x }\OperatorTok{*}\NormalTok{ (}\DecValTok{2}\OperatorTok{*}\NormalTok{k }\OperatorTok{{-}} \DecValTok{1}\NormalTok{) }\OperatorTok{/}\NormalTok{ (}\DecValTok{2}\OperatorTok{*}\NormalTok{k }\OperatorTok{+} \DecValTok{1}\NormalTok{)} \ControlFlowTok{return}\NormalTok{ s} \KeywordTok{def}\NormalTok{ pi\_machin(digits):} \NormalTok{ mp.dps }\OperatorTok{=}\NormalTok{ digits }\OperatorTok{+} \DecValTok{10} \NormalTok{ tol }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{10}\NormalTok{) }\OperatorTok{**}\NormalTok{ (}\OperatorTok{{-}}\NormalTok{digits) }\OperatorTok{/} \DecValTok{20} \NormalTok{ a }\OperatorTok{=}\NormalTok{ atan\_series(mp.mpf(}\DecValTok{1}\NormalTok{) }\OperatorTok{/} \DecValTok{5}\NormalTok{, tol)} \NormalTok{ b }\OperatorTok{=}\NormalTok{ atan\_series(mp.mpf(}\DecValTok{1}\NormalTok{) }\OperatorTok{/} \DecValTok{239}\NormalTok{, tol)} \ControlFlowTok{return} \DecValTok{16} \OperatorTok{*}\NormalTok{ a }\OperatorTok{{-}} \DecValTok{4} \OperatorTok{*}\NormalTok{ b} \end{Highlighting} \end{Shaded} The tolerance is divided by 20 because the final result combines the two arctangent approximations with coefficients 16 and 4. Machin-like formulas are much faster than the Gregory-Leibniz series and are relatively easy to implement. They are a good choice for moderate precision and for educational implementations. \subsection{6.4 Implementing the Brent-Salamin / AGM Method}\label{implementing-the-brent-salamin-agm-method} The Brent-Salamin method, based on the arithmetic-geometric mean, is one of the simplest high-precision algorithms. As noted in \textbf{5. Modern Computational Algorithms}, it has quadratic convergence, roughly doubling the number of correct digits at each iteration. The iteration is: \(a_0 = 1, \qquad b_0 = \frac{1}{\sqrt{2}}, \qquad t_0 = \frac{1}{4}, \qquad p_0 = 1,\) and for \(n \ge 0\), \(a_{n+1} = \frac{a_n + b_n}{2},\) \(b_{n+1} = \sqrt{a_n b_n},\) \(t_{n+1} = t_n - p_n (a_n - a_{n+1})^2,\) \(p_{n+1} = 2p_n.\) The approximation to \(\pi\) is \(\pi \approx \frac{(a_n + b_n)^2}{4t_n}.\) A straightforward implementation is: \begin{Shaded} \begin{Highlighting}[] \ImportTok{import}\NormalTok{ mpmath }\ImportTok{as}\NormalTok{ mp} \KeywordTok{def}\NormalTok{ pi\_agm(digits):} \NormalTok{ mp.dps }\OperatorTok{=}\NormalTok{ digits }\OperatorTok{+} \DecValTok{10} \NormalTok{ tol }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{10}\NormalTok{) }\OperatorTok{**}\NormalTok{ (}\OperatorTok{{-}}\NormalTok{digits)} \NormalTok{ a }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{1}\NormalTok{)} \NormalTok{ b }\OperatorTok{=} \DecValTok{1} \OperatorTok{/}\NormalTok{ mp.sqrt(}\DecValTok{2}\NormalTok{)} \NormalTok{ t }\OperatorTok{=}\NormalTok{ mp.mpf(}\StringTok{\textquotesingle{}0.25\textquotesingle{}}\NormalTok{)} \NormalTok{ p }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{1}\NormalTok{)} \ControlFlowTok{while} \VariableTok{True}\NormalTok{:} \NormalTok{ a\_new }\OperatorTok{=}\NormalTok{ (a }\OperatorTok{+}\NormalTok{ b) }\OperatorTok{/} \DecValTok{2} \NormalTok{ b\_new }\OperatorTok{=}\NormalTok{ mp.sqrt(a }\OperatorTok{*}\NormalTok{ b)} \NormalTok{ t }\OperatorTok{=}\NormalTok{ t }\OperatorTok{{-}}\NormalTok{ p }\OperatorTok{*}\NormalTok{ (a }\OperatorTok{{-}}\NormalTok{ a\_new) }\OperatorTok{**} \DecValTok{2} \NormalTok{ p }\OperatorTok{=} \DecValTok{2} \OperatorTok{*}\NormalTok{ p} \NormalTok{ a, b }\OperatorTok{=}\NormalTok{ a\_new, b\_new} \NormalTok{ pi }\OperatorTok{=}\NormalTok{ (a }\OperatorTok{+}\NormalTok{ b) }\OperatorTok{**} \DecValTok{2} \OperatorTok{/}\NormalTok{ (}\DecValTok{4} \OperatorTok{*}\NormalTok{ t)} \ControlFlowTok{if} \BuiltInTok{abs}\NormalTok{(a }\OperatorTok{{-}}\NormalTok{ b) }\OperatorTok{\textless{}}\NormalTok{ tol:} \ControlFlowTok{break} \ControlFlowTok{return}\NormalTok{ pi} \end{Highlighting} \end{Shaded} The stopping test \texttt{abs(a\ -\ b)\ \textless{}\ tol} is conservative. Because the method converges quadratically, a tighter test can be derived from the error analysis, but the simple test is easy to implement and reliable. Implementation considerations for the AGM method include: \begin{itemize} \tightlist \item Use high-precision square roots. \item Keep the working precision slightly above the target precision. \item Avoid unnecessary recomputation of \(\sqrt{2}\); compute it once at the required precision. \item The iteration is sequential, so parallelization is less natural than for series-based methods. \end{itemize} The Brent-Salamin method is attractive for implementation simplicity and numerical stability. It is often a good default for moderate to high precision when extremely large record-scale computations are not required. \subsection{6.5 Implementing the Chudnovsky Algorithm}\label{implementing-the-chudnovsky-algorithm} The Chudnovsky algorithm is a rapidly convergent series for \(1/\pi\). As discussed in \textbf{5. Modern Computational Algorithms}, it yields about 14 decimal digits per term and is especially suitable for high-precision computation. A standard form is \(\frac{1}{\pi} = 12 \sum_{k=0}^{\infty} (-1)^k \frac{(6k)! \left(13591409 + 545140134k\right)} {(3k)! (k!)^3 640320^{3k+3/2}}.\) Therefore, \(\pi = \frac{1}{12S},\) where \(S = \sum_{k=0}^{\infty} (-1)^k \frac{(6k)! \left(13591409 + 545140134k\right)} {(3k)! (k!)^3 640320^{3k+3/2}}.\) A naive implementation is conceptually simple: \begin{Shaded} \begin{Highlighting}[] \ImportTok{import}\NormalTok{ mpmath }\ImportTok{as}\NormalTok{ mp} \KeywordTok{def}\NormalTok{ pi\_chudnovsky\_naive(digits):} \NormalTok{ mp.dps }\OperatorTok{=}\NormalTok{ digits }\OperatorTok{+} \DecValTok{20} \NormalTok{ tol }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{10}\NormalTok{) }\OperatorTok{**}\NormalTok{ (}\OperatorTok{{-}}\NormalTok{digits) }\OperatorTok{/} \DecValTok{20} \NormalTok{ K }\OperatorTok{=} \BuiltInTok{int}\NormalTok{(digits }\OperatorTok{/} \DecValTok{14}\NormalTok{) }\OperatorTok{+} \DecValTok{2} \NormalTok{ S }\OperatorTok{=}\NormalTok{ mp.mpf(}\DecValTok{0}\NormalTok{)} \ControlFlowTok{for}\NormalTok{ k }\KeywordTok{in} \BuiltInTok{range}\NormalTok{(K }\OperatorTok{+} \DecValTok{1}\NormalTok{):} \NormalTok{ term }\OperatorTok{=}\NormalTok{ mp.factorial(}\DecValTok{6}\OperatorTok{*}\NormalTok{k) }\OperatorTok{*}\NormalTok{ (}\DecValTok{13591409} \OperatorTok{+} \DecValTok{545140134}\OperatorTok{*}\NormalTok{k)} \NormalTok{ term }\OperatorTok{/=}\NormalTok{ mp.factorial(}\DecValTok{3}\OperatorTok{*}\NormalTok{k) }\OperatorTok{*}\NormalTok{ (mp.factorial(k) }\OperatorTok{**} \DecValTok{3}\NormalTok{)} \NormalTok{ term }\OperatorTok{/=}\NormalTok{ mp.power(}\DecValTok{640320}\NormalTok{, }\DecValTok{3}\OperatorTok{*}\NormalTok{k }\OperatorTok{+}\NormalTok{ mp.mpf(}\StringTok{\textquotesingle{}1.5\textquotesingle{}}\NormalTok{))} \ControlFlowTok{if}\NormalTok{ k }\OperatorTok{\%} \DecValTok{2}\NormalTok{:} \NormalTok{ term }\OperatorTok{=} \OperatorTok{{-}}\NormalTok{term} \NormalTok{ S }\OperatorTok{+=}\NormalTok{ term} \ControlFlowTok{if} \BuiltInTok{abs}\NormalTok{(term) }\OperatorTok{\textless{}}\NormalTok{ tol:} \ControlFlowTok{break} \ControlFlowTok{return} \DecValTok{1} \OperatorTok{/}\NormalTok{ (}\DecValTok{12} \OperatorTok{*}\NormalTok{ S)} \end{Highlighting} \end{Shaded} This version is useful for small precision but is not efficient for large-scale computation. The factorials grow rapidly, and repeated factorial evaluation is expensive. Production implementations usually avoid direct factorial computation by using one of the following strategies: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \item \textbf{Term recurrence}\\ Compute successive terms using a rational recurrence instead of recomputing factorials. \item \textbf{Binary splitting}\\ Evaluate the sum using a divide-and-conquer strategy that combines partial products and partial sums. \item \textbf{Fast integer multiplication}\\ Use algorithms such as FFT-based or NTT-based multiplication for large integers. \end{enumerate} Binary splitting is especially important for record-scale computations. It reduces the number of large-number operations and allows the computation to be organized in a way that is more cache-friendly and more amenable to parallelization. The exact implementation depends on the arithmetic library, but the general idea is to compute the sum as a rational number \(S = \frac{A}{B}\) using large integers \(A\) and \(B\), and then perform the final division only once. For Chudnovsky, the number of terms needed for \(p\) decimal digits is approximately \(K \approx \frac{p}{14} + 1.\) This makes the algorithm far more efficient than classical series for high precision. \subsection{6.6 Loop Structures and Stopping Criteria}\label{loop-structures-and-stopping-criteria} Different algorithms require different loop structures and stopping criteria. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Method \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Typical loop structure \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Practical stopping criterion \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Archimedes polygon method & \texttt{for} over polygon side counts, often doubling \(n\) & Width of the lower/upper bound is below tolerance \\ Gregory-Leibniz series & \texttt{while} over terms & Next term is small enough that the alternating-series error bound is satisfied \\ Machin-like formulas & \texttt{while} for each arctangent series & Term magnitude is below a tolerance scaled by the final coefficient \\ Brent-Salamin / AGM & \texttt{while} over iterations & \(|a - b|\) is below tolerance, or a fixed number of quadratic-convergence iterations is used \\ Chudnovsky series & \texttt{for} or \texttt{while} over terms, often with binary splitting & Term magnitude is below tolerance, or a precomputed number of terms is used \\ \end{longtable} A robust implementation should include a maximum iteration count to prevent infinite loops in the event of a bug or an inappropriate tolerance. For alternating series, such as Gregory-Leibniz and Machin-like arctangent series, the error can be bounded by the first omitted term. For non-alternating or more complex series, such as Chudnovsky, the stopping criterion should be based on the observed term size and the desired relative or absolute error. When using arbitrary-precision floating-point arithmetic, it is important to distinguish between: \begin{itemize} \tightlist \item \textbf{Absolute tolerance}: useful when the final answer is expected to be near a known scale, such as \(\pi \approx 3.14\). \item \textbf{Relative tolerance}: useful when the quantity being computed may vary in magnitude. \item \textbf{Digit-based stopping}: stop when the number of verified correct digits reaches the target. \end{itemize} For \(\pi\), absolute tolerance is usually sufficient because the value is known to be near 3.14. \subsection{6.7 Practical Coding Considerations}\label{practical-coding-considerations} Several practical issues affect the reliability and performance of a \(\pi\)-computation program. \subsubsection{6.7.1 Precision and Guard Digits}\label{precision-and-guard-digits} Always use a working precision larger than the requested output precision. For example, if 100 decimal digits are required, compute with 110 or 120 digits and round the final result. \begin{Shaded} \begin{Highlighting}[] \NormalTok{target\_digits }\OperatorTok{=} \DecValTok{100} \NormalTok{guard\_digits }\OperatorTok{=} \DecValTok{20} \NormalTok{mp.dps }\OperatorTok{=}\NormalTok{ target\_digits }\OperatorTok{+}\NormalTok{ guard\_digits} \end{Highlighting} \end{Shaded} The final result should be rounded to the requested number of digits: \begin{Shaded} \begin{Highlighting}[] \NormalTok{result }\OperatorTok{=}\NormalTok{ mp.nstr(pi\_value, target\_digits)} \end{Highlighting} \end{Shaded} \subsubsection{6.7.2 Avoiding Overflow and Underflow}\label{avoiding-overflow-and-underflow} In native floating-point arithmetic, large factorials or powers can overflow. In arbitrary-precision integer arithmetic, overflow is not a problem, but memory usage can become large. For series with rapidly growing terms, use recurrences or binary splitting rather than direct evaluation of large factorials. \subsubsection{6.7.3 Rounding and Final Output}\label{rounding-and-final-output} The final value should be rounded consistently. If the program claims to produce \(p\) correct decimal digits, it should not print more digits than have been verified. A common validation string is: \begin{Shaded} \begin{Highlighting}[] \NormalTok{3.14159265358979323846264338327950288419716939937510} \end{Highlighting} \end{Shaded} This can be used as a unit test for the first 50 digits. \subsubsection{6.7.4 Library Selection}\label{library-selection} Common choices include: \begin{itemize} \tightlist \item \textbf{C/C++}: MPFR, GMP, Boost.Multiprecision. \item \textbf{Python}: mpmath, decimal, gmpy2. \item \textbf{Java}: \texttt{BigInteger}, \texttt{BigDecimal}, or external arbitrary-precision libraries. \item \textbf{Rust}: \texttt{num-bigint}, \texttt{num-integer}, or custom fixed-point arithmetic. \end{itemize} For high-performance implementations, the multiplication algorithm is often the dominant cost. Fast multiplication, as discussed in \textbf{8. Performance and Optimization}, is essential for computing millions of digits. \subsubsection{6.7.5 Numerical Stability}\label{numerical-stability} Most of the methods discussed here are numerically stable when implemented with sufficient guard digits. However, care is still needed: \begin{itemize} \tightlist \item The Gregory-Leibniz series is stable but slow. \item Machin-like formulas are stable and much faster. \item The AGM method is stable and simple. \item Chudnovsky is stable when implemented with exact integer arithmetic or sufficient precision. \end{itemize} Catastrophic cancellation is not a major issue for these particular \(\pi\) algorithms, but it can appear in auxiliary computations, such as high-precision square roots or trigonometric evaluations. \subsubsection{6.7.6 Testing and Validation}\label{testing-and-validation} A good implementation should include tests such as: \begin{itemize} \tightlist \item Compare the result against known digits of \(\pi\). \item Test at several precision levels. \item Verify that increasing the precision does not change already verified digits. \item Check that the error decreases as expected for each method. \item Test edge cases, such as very small precision and very large precision. \item Validate that the stopping criterion is not triggered prematurely. \end{itemize} For example: \begin{Shaded} \begin{Highlighting}[] \NormalTok{known }\OperatorTok{=} \StringTok{"3.14159265358979323846264338327950288419716939937510"} \KeywordTok{def}\NormalTok{ test\_pi(digits):} \NormalTok{ value }\OperatorTok{=}\NormalTok{ pi\_agm(digits)} \NormalTok{ text }\OperatorTok{=}\NormalTok{ mp.nstr(value, digits)} \ControlFlowTok{assert}\NormalTok{ text.startswith(known[:digits])} \end{Highlighting} \end{Shaded} \subsubsection{6.7.7 Code Organization}\label{code-organization} A clean implementation separates concerns: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item \textbf{Arithmetic layer}: precision setup, integer or multiprecision operations. \item \textbf{Algorithm layer}: series, iteration, or binary splitting logic. \item \textbf{Validation layer}: known-digit checks and error estimates. \item \textbf{I/O layer}: formatting and output. \end{enumerate} This separation makes it easier to switch between algorithms, change precision, or port the code to another language. \subsection{6.8 Summary of Implementation Choices}\label{summary-of-implementation-choices} For low-precision or educational implementations, the Gregory-Leibniz series or a Machin-like formula is often sufficient. The Gregory-Leibniz series is the simplest to understand, while Machin-like formulas provide much better convergence with only modest additional complexity. For moderate to high precision, the Brent-Salamin / AGM method is an excellent choice because it is simple, stable, and converges quadratically. It is often easier to implement correctly than record-scale series algorithms. For very high precision, the Chudnovsky algorithm is usually preferred, especially when combined with binary splitting and fast multiplication. It requires more implementation effort but is far more efficient for large digit counts. In all cases, the key implementation principles are: \begin{itemize} \tightlist \item Use a numeric type appropriate to the target precision. \item Add guard digits to the working precision. \item Use stable recurrences rather than direct evaluation of large factorials or powers when possible. \item Choose stopping criteria based on the convergence properties of the method. \item Validate the result against known digits of \(\pi\). \item Optimize only after correctness has been established. \end{itemize} \section{7. Convergence and Error Analysis}\label{convergence-and-error-analysis} \subsection{7.1 Convergence Rates}\label{convergence-rates} A central question in numerical pi calculation is how quickly an approximation approaches the true value as the number of terms, iterations, or polygon sides increases. Let \(p\) denote the number of requested decimal digits. A method is useful for high-precision computation if the number of required operations grows slowly with \(p\), ideally logarithmically or linearly rather than exponentially. The main methods discussed in this publication have very different convergence behavior: \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.2727}} >{\raggedleft\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3636}} >{\raggedleft\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3636}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Method \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Convergence behavior \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Approximate work to obtain \(p\) digits \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Archimedes polygon method & \(O(n^{-2})\) in number of sides \(n\) & \(n \sim 10^{p/2}\) \\ Gregory-Leibniz series & \(O(N^{-1})\) in number of terms \(N\) & \(N \sim 4\cdot 10^{p}\) \\ Machin-like arctangent formulas & Exponential in number of terms & \(N \sim 0.72p\) for Machin's formula \\ Brent-Salamin / AGM iteration & Quadratic & \(n \sim \log_2 p\) iterations \\ Chudnovsky series & Exponential, about 14.18 digits per term & \(N \sim 0.071p\) terms \\ \end{longtable} \subsubsection{Polygonal bounds}\label{polygonal-bounds} For a regular \(n\)-gon inscribed in and circumscribed about the unit circle, \(L_n = n\sin\left(\frac{\pi}{n}\right), \qquad U_n = n\tan\left(\frac{\pi}{n}\right),\) with \(L_n < \pi < U_n.\) Using Taylor expansions, \(L_n = \pi - \frac{\pi^3}{6n^2} + O(n^{-4}),\) and \(U_n = \pi + \frac{\pi^3}{3n^2} + O(n^{-4}).\) Thus both the lower and upper polygonal bounds have error of order \(O(n^{-2})\). Doubling the number of sides reduces the error by roughly a factor of four. If one uses the midpoint \(M_n = \frac{L_n+U_n}{2},\) the leading \(n^{-2}\) error cancels, giving an approximation with error \(O(n^{-4})\). Even so, the polygon method is far less efficient than modern series or iterative methods for high precision. To obtain \(p\) decimal digits, one needs roughly \(n \sim 10^{p/2}\) sides, or \(n \sim 10^{p/4}\) if using the midpoint. For \(p=100\), this is already on the order of \(10^{50}\) sides, which is impractical. \subsubsection{Gregory-Leibniz series}\label{gregory-leibniz-series-1} The Gregory-Leibniz series is \(\frac{\pi}{4} = 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots.\) If \(S_N = \sum_{k=0}^{N} \frac{(-1)^k}{2k+1},\) then the alternating-series remainder gives \(\left|\frac{\pi}{4}-S_N\right| \le \frac{1}{2N+3},\) so \(|\pi - 4S_N| \le \frac{4}{2N+3}.\) The convergence is only linear in \(N\), with error \(O(N^{-1})\). To obtain \(p\) correct decimal digits, one needs approximately \(N \approx 4\cdot 10^p\) terms. For example, 10 digits require about \(4\times 10^9\) terms, while 100 digits would require about \(4\times 10^{100}\) terms. This makes the Gregory-Leibniz series mathematically simple but computationally inefficient for high precision. \subsubsection{Machin-like formulas}\label{machin-like-formulas-1} Machin-like formulas improve the convergence of arctangent series by using small rational arguments. The classical Machin formula is \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right).\) For a general arctangent series, \(\arctan\left(\frac{1}{q}\right) = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)q^{2k+1}}, \qquad q\ge 1.\) If the sum is truncated after the term \(k=N\), the remainder satisfies \(\left| \arctan\left(\frac{1}{q}\right) - \sum_{k=0}^{N} \frac{(-1)^k}{(2k+1)q^{2k+1}} \right| \le \frac{1}{(2N+3)q^{2N+3}}.\) For Machin's formula, the error is dominated by the \(\arctan(1/5)\) term: \(E_N \le \frac{4}{(2N+3)5^{2N+3}} + \frac{1}{(2N+3)239^{2N+3}}.\) Asymptotically, the error decreases by a factor of about \(25\) per term, giving approximately \(\log_{10}(25) \approx 1.398\) decimal digits per term. Thus Machin-like formulas require only on the order of \(0.72p\) terms to obtain \(p\) digits. More aggressive Machin-like identities with smaller arctangent arguments can improve this rate further, although the coefficients may become larger. \subsubsection{Brent-Salamin / AGM iteration}\label{brent-salamin-agm-iteration} The Brent-Salamin algorithm, based on the arithmetic-geometric mean, has quadratic convergence. With \(a_{n+1} = \frac{a_n+b_n}{2}, \qquad b_{n+1} = \sqrt{a_nb_n}, \qquad c_n = a_n-b_n,\) and initial values \(a_0=1, \qquad b_0=\frac{1}{\sqrt{2}},\) one obtains a pi estimate of the form \(\pi_n = \frac{2a_n^2} {1-\sum_{k=0}^{n-1}2^{k+1}c_k^2}.\) The error satisfies, asymptotically, \(e_{n+1} \approx K e_n^2,\) where \(e_n = |\pi-\pi_n|\) and \(K\) is a method-dependent constant. Consequently, the number of correct digits roughly doubles at each iteration. To obtain \(p\) digits, only about \(n \approx \log_2 p\) iterations are needed. For example, 100 digits require only about 7 iterations. \subsubsection{Chudnovsky series}\label{chudnovsky-series} The Chudnovsky series is one of the fastest known series for high-precision pi calculation. It computes \(1/\pi\) via \(\frac{1}{\pi} = 12 \sum_{k=0}^{\infty} (-1)^k \frac{(6k)!\left(13591409+545140134k\right)} {(3k)!(k!)^3\,640320^{3k+3/2}}.\) The magnitude of successive terms decreases asymptotically by the factor \(\frac{1728}{640320^3}.\) Therefore the number of decimal digits gained per term is \(-\log_{10}\left(\frac{1728}{640320^3}\right) \approx 14.18.\) Thus only about \(N \approx \frac{p}{14.18}\) terms are needed for \(p\) digits. For 100 digits, this is only about 8 terms. The main practical difficulty is not the number of terms but the efficient evaluation of the large factorial and power factors, which is normally handled using term recurrences, binary splitting, and fast multiplication. \subsection{7.2 Error Bounds and Stopping Criteria}\label{error-bounds-and-stopping-criteria} A reliable pi algorithm should provide a way to know when the desired precision has been reached. For many of the methods considered here, the error can be bounded rigorously or controlled with high confidence using guard digits. Let the requested precision be \(p\) decimal digits. A common absolute tolerance is \(\tau = \frac{1}{2}10^{-p}.\) In practice, one usually works with \(p+g\) guard digits, where \(g\) is a small positive integer, and uses \(\tau_g = \frac{1}{2}10^{-(p+g)}\) as the internal stopping tolerance. \subsubsection{Polygonal bounds}\label{polygonal-bounds-1} For the polygon method, the Taylor expansions give \(0 < \pi - L_n = \frac{\pi^3}{6n^2} + O(n^{-4}),\) and \(0 < U_n - \pi = \frac{\pi^3}{3n^2} + O(n^{-4}).\) Thus the interval \([L_n,U_n]\) has width \(U_n-L_n = O(n^{-2}).\) A rigorous stopping criterion can be obtained by requiring \(U_n-L_n < \tau_g.\) If the midpoint \(M_n=(L_n+U_n)/2\) is used, the error is \(O(n^{-4})\), but a simple rigorous bound is less immediate. In any case, the polygon method is rarely competitive for high precision because the required number of sides grows rapidly with \(p\). \subsubsection{Gregory-Leibniz series}\label{gregory-leibniz-series-2} For the Gregory-Leibniz series, the alternating-series test gives a simple rigorous bound. If \(S_N = \sum_{k=0}^{N} \frac{(-1)^k}{2k+1},\) then \(|\pi - 4S_N| \le \frac{4}{2N+3}.\) Therefore one may stop when \(\frac{4}{2N+3} < \tau_g.\) This bound is easy to use, but the required \(N\) grows exponentially in \(p\). The method is therefore suitable mainly for low-precision or educational use. \subsubsection{Machin-like formulas}\label{machin-like-formulas-2} For a Machin-like formula of the form \(\frac{\pi}{4} = \sum_{j=1}^{m} c_j \arctan\left(\frac{1}{q_j}\right),\) where the \(c_j\) are integers and the \(q_j\) are positive integers, the error is bounded by the sum of the individual arctangent remainders. If each arctangent series is truncated after the term \(k=N_j\), then \(\left| \frac{\pi}{4} - \sum_{j=1}^{m} c_j \sum_{k=0}^{N_j} \frac{(-1)^k}{(2k+1)q_j^{2k+1}} \right| \le \sum_{j=1}^{m} \frac{|c_j|}{(2N_j+3)q_j^{2N_j+3}}.\) For Machin's formula, \(\frac{\pi}{4} = 4\arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right),\) using the same truncation index \(N\) in both arctangent series gives \(E_N \le \frac{4}{(2N+3)5^{2N+3}} + \frac{1}{(2N+3)239^{2N+3}}.\) A safe stopping criterion is \(E_N < \tau_g.\) Because the remainder is dominated by the smallest denominator \(q_j\), the convergence rate is controlled primarily by the smallest \(q_j\) appearing in the formula. \subsubsection{Brent-Salamin / AGM iteration}\label{brent-salamin-agm-iteration-1} For the AGM-based Brent-Salamin method, the correction sum \(\sum_{k=0}^{\infty} 2^{k+1}c_k^2\) converges extremely rapidly because \(c_k\) decreases quadratically. If the sum is truncated after \(n-1\), the next omitted correction term is \(2^{n+1}c_n^2.\) A practical stopping criterion is therefore \(2^{n+1}c_n^2 < \tau_g.\) In addition, one should ensure that \(a_n\) and \(b_n\) are computed with enough guard digits to make rounding error negligible. Because the method is quadratically convergent, once the correction term is below the tolerance, the remaining error is usually far smaller than the requested precision. \subsubsection{Chudnovsky series}\label{chudnovsky-series-1} Let \(U_k\) denote the \(k\)-th term in the Chudnovsky series for \(1/\pi\), including the factor 12. If \(S_N = \sum_{k=0}^{N} U_k,\) then, because the series is alternating with decreasing term magnitudes, \(\left| \frac{1}{\pi} - S_N \right| \le |U_{N+1}|.\) If one defines \(\pi_N = \frac{1}{S_N},\) then the error in \(\pi_N\) is \(\left| \pi - \pi_N \right| = \left| \frac{1}{S_N+R_N} - \frac{1}{S_N} \right| = \frac{|R_N|}{S_N(S_N+R_N)},\) where \(R_N\) is the remainder. A conservative bound is \(\left| \pi - \pi_N \right| \le \frac{|U_{N+1}|}{(S_N-|U_{N+1}|)^2},\) provided \(S_N > |U_{N+1}|\). Since \(S_N \approx 1/\pi\), a practical rule is to require \(|U_{N+1}| < \frac{\tau_g}{\pi^2}.\) This gives a rigorous or near-rigorous stopping criterion, depending on how the intermediate arithmetic is rounded. \subsection{7.3 Numerical Stability}\label{numerical-stability-1} Convergence rate alone does not determine practical accuracy. A method may converge rapidly in exact arithmetic but become unstable when implemented with finite precision. The main stability issues are rounding error, cancellation, overflow, and the accumulation of many small operations. \subsubsection{Floating-point limitations}\label{floating-point-limitations} Standard double-precision floating-point arithmetic provides about 15 to 16 decimal digits of precision. Therefore, no method can reliably produce more than about 16 correct digits using double precision alone, regardless of its theoretical convergence rate. For higher precision, arbitrary-precision integer, rational, or multiprecision floating-point arithmetic is required. \subsubsection{Gregory-Leibniz series}\label{gregory-leibniz-series-3} The Gregory-Leibniz series is stable in the sense that it does not involve severe cancellation, but it requires a very large number of additions. In finite precision, each addition introduces a small rounding error. When the number of terms is enormous, these errors can accumulate. Moreover, once the terms become smaller than the machine epsilon relative to the current partial sum, further additions no longer change the result. Thus, in double precision, the series cannot be used to obtain more than about 16 digits, and in arbitrary precision it remains inefficient because of its slow convergence. \subsubsection{Machin-like formulas}\label{machin-like-formulas-3} Machin-like formulas are numerically stable and much more efficient than the Gregory-Leibniz series. The arctangent terms decrease rapidly, especially when the arguments \(1/q_j\) are small. A stable implementation computes the terms recursively, for example \(t_{k+1} = -\frac{t_k}{q^2},\) rather than repeatedly evaluating powers and factorials. This avoids unnecessary overflow and reduces rounding error. Because the number of terms is only proportional to the desired precision, a modest number of guard digits is usually sufficient. \subsubsection{Polygonal methods}\label{polygonal-methods} The polygon method can be implemented using trigonometric functions or recurrence relations. If trigonometric functions are used, the small-angle argument \(\pi/n\) must be handled accurately, especially for large \(n\). If recurrence relations are used, rounding error can accumulate over many doublings of the polygon side count. The method is stable for moderate \(n\), but it is not competitive for high precision because the number of sides grows rapidly with the desired number of digits. \subsubsection{Brent-Salamin / AGM iteration}\label{brent-salamin-agm-iteration-2} The AGM iteration is generally very stable. The sequences \(a_n\) and \(b_n\) remain positive and converge monotonically to the arithmetic-geometric mean. The main delicate operation is \(c_n = a_n - b_n,\) because \(a_n\) and \(b_n\) become very close as the iteration proceeds. This subtraction can lose relative precision in \(c_n\). However, what matters for the final pi estimate is the absolute error in \(c_n^2\), and the correction terms decrease so rapidly that a few guard digits are usually enough to make the rounding error negligible. A practical implementation should: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item use \(p+g\) working digits, where \(g\) is at least 5 and often larger for high precision; \item compute square roots with accuracy comparable to the working precision; \item stop when the next correction term \(2^{n+1}c_n^2\) is below the internal tolerance; \item validate the final result against known digits of \(\pi\). \end{enumerate} \subsubsection{Chudnovsky series}\label{chudnovsky-series-2} The Chudnovsky series is extremely fast in exact arithmetic, but naive floating-point evaluation can be problematic. The terms involve large factorials and powers, and direct evaluation may cause overflow, underflow, or loss of precision. Production implementations therefore avoid naive factorial computation. Instead, they use: \begin{itemize} \tightlist \item term recurrences to update successive terms; \item binary splitting to evaluate the sum efficiently; \item arbitrary-precision integer arithmetic for the numerator and denominator; \item fast multiplication for large integers; \item a final high-precision division to obtain \(\pi\). \end{itemize} With these techniques, the Chudnovsky algorithm is both accurate and stable. The main source of error is the final division and the chosen working precision, not the convergence of the series itself. \subsection{7.4 Accuracy versus Computational Cost}\label{accuracy-versus-computational-cost} To compare methods fairly, it is useful to express computational cost in terms of the cost of multiplying two \(p\)-digit numbers. Let \(M(p)\) denote this cost. For example, with schoolbook multiplication \(M(p)=O(p^2)\), while with fast Fourier-transform-based multiplication \(M(p)=O(p\log p)\) up to logarithmic factors. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.1667}} >{\raggedleft\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.2222}} >{\raggedleft\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.2222}} >{\raggedleft\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.2222}} >{\raggedright\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.1667}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Method \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Number of terms/iterations for \(p\) digits \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Naive cost \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Optimized cost \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Practical role \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Polygon method & \(n \sim 10^{p/2}\) & \(O(10^{p/2}M(p))\) & Not competitive & Historical, rigorous bounds \\ Gregory-Leibniz & \(N \sim 4\cdot 10^p\) & \(O(10^p M(p))\) & Not practical for high \(p\) & Educational, low precision \\ Machin-like & \(N \sim 0.72p\) & \(O(pM(p))\) & \(O(M(p)\log p)\) with binary splitting & Simple, stable, moderate precision \\ Brent-Salamin / AGM & \(n \sim \log_2 p\) & \(O(\log p\, M(p))\) & \(O(M(p)\log p)\) with fast square root & General-purpose, stable \\ Chudnovsky & \(N \sim 0.071p\) & \(O(pM(p))\) naive & \(O(M(p)\log p)\) with binary splitting & Very high precision, record-scale \\ \end{longtable} The Gregory-Leibniz series and polygon method have simple error bounds but poor computational scaling. Their convergence is too slow for high-precision work. Machin-like formulas are a major improvement: they are simple to implement, have rigorous alternating-series error bounds, and require only a linear number of terms in the desired precision. The Brent-Salamin / AGM method is attractive because it requires only logarithmically many iterations and is stable. It is often a good general-purpose choice, especially when implementation simplicity and reliability are important. The Chudnovsky algorithm is usually preferred for very large-scale computations because it produces about 14 decimal digits per term. When combined with binary splitting and fast multiplication, its asymptotic cost is favorable and its constant factor is small. For record-breaking pi calculations, Chudnovsky-type series are typically the method of choice. \subsection{7.5 Practical Recommendations}\label{practical-recommendations} The choice of method depends on the target precision, the available arithmetic environment, and the desired balance between simplicity and efficiency. \begin{itemize} \item For \(p \le 15\), double-precision implementations of Machin-like formulas or the Brent-Salamin method are usually sufficient. The Brent-Salamin method is often preferable because it requires only a few iterations. \item For moderate precision, say \(16 \le p \le 10^3\), Machin-like formulas and the Brent-Salamin method are both practical. Machin-like formulas are simpler to code, while Brent-Salamin is more stable and requires fewer iterations. \item For high precision, \(p \ge 10^4\), the Chudnovsky algorithm with binary splitting and fast multiplication is usually the most efficient. The Brent-Salamin method remains a strong alternative when implementation simplicity is prioritized. \item For all methods, use guard digits. A working precision of \(p+g\) digits, with \(g\) between 5 and 20 depending on the method and precision, is a safe practice. \item Use rigorous or conservative stopping criteria. For alternating series, the first omitted term provides a natural bound. For the AGM method, the next correction term \(2^{n+1}c_n^2\) provides a practical stopping test. For Chudnovsky, the first omitted term in the \(1/\pi\) series, scaled by approximately \(\pi^2\), gives a reliable error estimate. \item Validate the final result against known digits of \(\pi\). This is especially important in high-precision implementations, where subtle bugs in multiplication, division, square root, or binary splitting can produce plausible but incorrect results. \end{itemize} In summary, classical methods such as the polygon method and the Gregory-Leibniz series are valuable for understanding the mathematical foundations of pi approximation, but their slow convergence makes them unsuitable for high-precision computation. Machin-like formulas provide a simple and stable improvement. The Brent-Salamin / AGM method offers quadratic convergence and excellent stability. The Chudnovsky algorithm provides the highest asymptotic efficiency for very large precision and is the preferred method for record-scale pi calculations. \section{8. Performance and Optimization}\label{performance-and-optimization} \subsection{8.1 Time Complexity and Cost per Digit}\label{time-complexity-and-cost-per-digit} The performance of a large-scale calculation of \(\pi\) depends on two coupled factors: the convergence rate of the chosen algorithm and the cost of the high-precision arithmetic required to evaluate it. As summarized in \textbf{7. Convergence and Error Analysis}, the convergence behavior differs strongly among methods. Classical methods are mathematically simple but computationally inefficient, whereas modern algorithms reduce the number of required terms or iterations dramatically. Let \(D\) denote the number of requested decimal digits and let \(p = D + g\) be the working precision, where \(g\) is the number of guard digits. Let \(M(p)\) denote the cost of multiplying two \(p\)-digit integers. The value of \(M(p)\) depends on the multiplication algorithm: \(M(p) = \begin{cases} \Theta(p^2), & \text{schoolbook multiplication}, \\ \Theta(p^{\log_2 3}), & \text{Karatsuba multiplication}, \\ \Theta(p \log p \log \log p), & \text{FFT/NTT-based multiplication, asymptotically}. \end{cases}\) For high-precision \(\pi\) calculations, the total running time is usually dominated by a small number of very large multiplications, divisions, or square roots rather than by the number of loop iterations alone. \subsubsection{Classical methods}\label{classical-methods} The polygon method of Archimedes converges as \(O(n^{-2})\), where \(n\) is the number of polygon sides. To obtain \(D\) correct decimal digits, one needs roughly \(n \sim 10^{D/2},\) which grows exponentially in \(D\). This makes the method impractical for high precision. The Gregory-Leibniz series converges only as \(O(N^{-1})\), where \(N\) is the number of terms. To obtain \(D\) digits, one needs approximately \(N \sim 10^D\) terms. Even if each term were cheap, the total number of operations grows exponentially with the requested precision. Therefore, the Gregory-Leibniz series is useful mainly for educational purposes or very low-precision demonstrations. \subsubsection{Machin-like formulas}\label{machin-like-formulas-4} Machin-like formulas converge exponentially because they use arctangent series with small rational arguments. As noted in \textbf{7. Convergence and Error Analysis}, Machin's formula gives roughly \(1.4\) decimal digits per term. Thus the number of terms required for \(D\) digits is approximately \(N_{\text{Machin}} \approx \frac{D}{1.4}.\) With a straightforward fixed-precision implementation, each term update costs roughly \(O(p)\) arithmetic operations, giving a total cost on the order of \(O(Dp) \approx O(D^2)\) when \(p \sim D\). This is far better than the Gregory-Leibniz series, but still less efficient than the best modern algorithms for very large \(D\). \subsubsection{Brent-Salamin / AGM method}\label{brent-salamin-agm-method} The Brent-Salamin or arithmetic-geometric mean method has quadratic convergence: the number of correct digits roughly doubles at each iteration. Therefore, the number of iterations required to reach \(D\) digits is approximately \(O(\log D).\) Each iteration involves a small number of high-precision multiplications, additions, subtractions, and square roots. If the working precision is \(p\), the total cost is often modeled as \(O(M(p)\log D).\) This makes the Brent-Salamin method highly efficient for moderate to high precision and attractive for implementation because it is stable and relatively simple, as discussed in \textbf{6. Algorithm Implementation}. \subsubsection{Chudnovsky algorithm}\label{chudnovsky-algorithm} The Chudnovsky algorithm is especially efficient for very high precision because it yields about \(14.18\) decimal digits per term. The number of terms required for \(D\) digits is therefore approximately \(N_{\text{Chud}} \approx \frac{D}{14.18}.\) A naive implementation that evaluates factorials directly is inefficient and may overflow or require enormous intermediate integers. Production implementations instead use term recurrences, binary splitting, and fast multiplication, as emphasized in \textbf{5. Modern Computational Algorithms} and \textbf{6. Algorithm Implementation}. With binary splitting, the summation can be organized as a divide-and-conquer computation. The total cost is commonly modeled as \(O(M(p)\log D),\) up to constants and implementation-dependent overhead. When combined with FFT- or NTT-based multiplication, this gives a very favorable asymptotic cost for record-scale calculations. \subsection{8.2 Memory Usage}\label{memory-usage} Memory usage in high-precision \(\pi\) calculation is governed primarily by the size of the arbitrary-precision integers or multiprecision floating-point numbers that must be stored. If \(p\) decimal digits of working precision are used, a single large integer requires approximately \(p \log_2 10 \approx 3.322p\) bits, or about \(0.415p\) bytes, before accounting for library overhead, alignment, and temporary buffers. For example: \begin{longtable}[]{@{}rr@{}} \toprule\noalign{} Requested digits \(D\) & Approx. size of one \(D\)-digit integer \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot \(10^6\) & \(\approx 0.4\) MB \\ \(10^7\) & \(\approx 4\) MB \\ \(10^8\) & \(\approx 41\) MB \\ \(10^9\) & \(\approx 415\) MB \\ \end{longtable} In practice, a calculation may require several such integers simultaneously, so the total memory footprint can be several times larger than the size of a single number. \subsubsection{Brent-Salamin / AGM method}\label{brent-salamin-agm-method-1} The Brent-Salamin method requires only a small number of large values at each iteration, typically \(a_n\), \(b_n\), and \(c_n\), together with temporary values for multiplication and square root. Its memory usage is therefore roughly \(O(p),\) which is modest compared with more complex summation algorithms. \subsubsection{Chudnovsky algorithm with binary splitting}\label{chudnovsky-algorithm-with-binary-splitting} Binary splitting is memory-efficient if implemented with a depth-first traversal. The recursion stack has depth \(O(\log D)\), and the largest intermediate integers have size \(O(p)\). However, temporary buffers for multiplication, division, or FFT-based arithmetic may require additional memory. If all subproblem results are stored simultaneously, memory usage can increase substantially. A well-designed implementation therefore computes and combines subresults incrementally, keeping only the necessary intermediate values. \subsubsection{Classical and series-based methods}\label{classical-and-series-based-methods} The Gregory-Leibniz series and simple Machin-like implementations require relatively little memory, often only \(O(p)\), because they maintain a running sum and a current term. Their limitation is not memory but time: they require too many terms for high precision. \subsubsection{Large-scale considerations}\label{large-scale-considerations} For record-scale calculations, memory management becomes a first-class performance concern. Practical strategies include: \begin{itemize} \tightlist \item using optimized arbitrary-precision libraries with efficient memory allocation; \item reusing buffers for large integers; \item avoiding unnecessary copies of large numbers; \item checkpointing intermediate results to disk when memory is insufficient; \item monitoring peak memory usage, not just final result size; \item choosing multiplication algorithms whose temporary memory requirements fit the available hardware. \end{itemize} \subsection{8.3 Parallelization}\label{parallelization} Parallelization can improve the performance of \(\pi\) calculations at two levels: algorithmic parallelism and arithmetic-level parallelism. \subsubsection{Algorithmic parallelism}\label{algorithmic-parallelism} The Chudnovsky algorithm is well suited to parallelization when implemented with binary splitting. The summation range can be divided into independent subranges, each of which can be evaluated on a separate processor or node. The partial results are then combined in a tree-like reduction. This approach is effective because the subproblems are largely independent and can be balanced according to available compute resources. However, the final combination step is sequential and may become a bottleneck for very large calculations. Machin-like formulas can also be parallelized by summing blocks of arctangent terms independently. Each block can be computed in parallel and then combined. Care must be taken to preserve the error bounds and stopping criteria described in \textbf{7. Convergence and Error Analysis}, especially when using alternating series. The Brent-Salamin method is less amenable to algorithmic parallelization because each iteration depends on the previous one. The recurrence is inherently sequential: \(a_{n+1}, b_{n+1}, c_{n+1}\) depend on \(a_n\), \(b_n\), and \(c_n\). Therefore, parallel speedup for Brent-Salamin usually comes from parallelizing the internal high-precision operations rather than from parallelizing the iteration itself. \subsubsection{Arithmetic-level parallelism}\label{arithmetic-level-parallelism} Even when the algorithm is sequential, the underlying high-precision arithmetic can be parallelized. Fast multiplication algorithms such as Karatsuba, Toom-Cook, and FFT/NTT-based methods can exploit multiple cores. FFT-based multiplication in particular can be highly parallelizable because it involves many independent transform operations. Other operations that may benefit from parallelism include: \begin{itemize} \tightlist \item high-precision division; \item square root computation; \item Newton-Raphson refinement; \item modular reduction; \item remainder computation; \item final conversion from \(1/\pi\) to \(\pi\). \end{itemize} \subsubsection{Distributed computation}\label{distributed-computation} For extremely large calculations, distributed systems may be used. In such settings, the main challenges are: \begin{itemize} \tightlist \item communication overhead between nodes; \item load balancing; \item checkpointing and fault tolerance; \item deterministic reproducibility; \item synchronization of partial results. \end{itemize} Distributed computation is most useful when the problem size exceeds the memory or compute capacity of a single machine. \subsection{8.4 Optimization Strategies}\label{optimization-strategies} Optimizing a large-scale \(\pi\) calculation requires choosing the right algorithm, arithmetic representation, and implementation techniques. \subsubsection{Algorithm selection}\label{algorithm-selection} The choice of algorithm should match the target precision: \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}} >{\raggedright\arraybackslash}p{(\columnwidth - 4\tabcolsep) * \real{0.3333}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Precision range \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Recommended approach \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Reason \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Low precision & Double-precision floating point, simple series & Simplicity and speed \\ Low to moderate precision & Machin-like formulas & Fast convergence, simple implementation \\ Moderate to high precision & Brent-Salamin / AGM & Quadratic convergence, stability, simplicity \\ Very high precision & Chudnovsky with binary splitting & High digits per term, efficient large-integer arithmetic \\ \end{longtable} This hierarchy is consistent with the implementation guidance in \textbf{6. Algorithm Implementation} and the performance comparison in \textbf{7. Convergence and Error Analysis}. \subsubsection{Use arbitrary-precision arithmetic}\label{use-arbitrary-precision-arithmetic} For high precision, native floating-point types are insufficient. Arbitrary-precision integers, rationals, or multiprecision floating-point libraries should be used. The working precision should include guard digits above the requested output precision to reduce accumulated rounding error. For the Chudnovsky algorithm, exact integer arithmetic is often preferred. The series is evaluated for \(1/\pi\), and the final division produces \(\pi\). This avoids some of the rounding complications associated with floating-point summation. \subsubsection{Avoid naive factorial evaluation}\label{avoid-naive-factorial-evaluation} A direct implementation of the Chudnovsky series that computes factorials explicitly is inefficient. Instead, production implementations use term recurrences that update each term from the previous one using multiplications and divisions by smaller integers. This reduces both time and memory usage. \subsubsection{Binary splitting}\label{binary-splitting} Binary splitting is one of the most important optimizations for high-precision series evaluation. It computes partial sums and products recursively, reducing repeated work and improving numerical stability. For the Chudnovsky algorithm, binary splitting allows the computation to be organized as a balanced tree. This improves both cache behavior and parallelizability. It also makes it easier to control intermediate sizes and avoid unnecessary growth of integers. \subsubsection{Fast multiplication}\label{fast-multiplication} Because high-precision multiplication dominates the cost of many \(\pi\) algorithms, choosing an efficient multiplication method is critical. Common choices include: \begin{itemize} \tightlist \item schoolbook multiplication for small operands; \item Karatsuba multiplication for medium-sized operands; \item Toom-Cook multiplication for larger operands; \item FFT- or NTT-based multiplication for very large operands. \end{itemize} In practice, optimized libraries often use a hybrid strategy, switching between algorithms depending on operand size. \subsubsection{Precision management}\label{precision-management} The working precision should be chosen carefully. Too little precision leads to incorrect results, while too much precision wastes time and memory. A typical strategy is: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item determine the requested output precision \(D\); \item add a safety margin \(g\) of guard digits; \item compute at precision \(p = D + g\); \item round or truncate the final result to \(D\) digits; \item validate the result against known digits of \(\pi\). \end{enumerate} The value of \(g\) depends on the algorithm and implementation. For stable methods such as Brent-Salamin, a modest guard margin may suffice. For series-based methods, a larger margin may be needed to control accumulated error. \subsubsection{Stopping criteria}\label{stopping-criteria} Stopping criteria should be matched to the convergence behavior of the algorithm. \begin{itemize} \tightlist \item For the Gregory-Leibniz series, an alternating-series bound can be used. \item For Machin-like formulas, the remainders of the arctangent series can be bounded. \item For Brent-Salamin, the correction term \(2^{n+1}c_n^2\) provides a practical stopping criterion. \item For Chudnovsky, the first omitted term in the \(1/\pi\) series, scaled appropriately, gives a reliable error estimate. \end{itemize} These criteria are discussed in \textbf{7. Convergence and Error Analysis}. \subsubsection{Use optimized libraries}\label{use-optimized-libraries} For serious high-precision computation, using a well-tested arbitrary-precision library is usually preferable to writing all low-level arithmetic from scratch. Libraries such as GMP, MPFR, FLINT, or similar systems provide optimized multiplication, division, square root, and memory management routines. Using such libraries allows the implementation to focus on the mathematical algorithm while relying on highly optimized arithmetic kernels. \subsubsection{Checkpointing and validation}\label{checkpointing-and-validation} Large calculations may run for hours, days, or longer. Checkpointing intermediate results allows the computation to be resumed after interruption. Validation is also essential. Common validation strategies include: \begin{itemize} \tightlist \item comparing the result against known digits of \(\pi\); \item recomputing with a different algorithm; \item recomputing at higher precision; \item checking error bounds; \item verifying intermediate checksums or modular residues. \end{itemize} \subsection{8.5 Practical Performance Summary}\label{practical-performance-summary} For small-scale or educational calculations, simplicity is usually more important than asymptotic efficiency. A double-precision implementation of a Machin-like formula or the Brent-Salamin method is often sufficient. For high-precision calculations, performance is dominated by the cost of large-integer arithmetic. The Brent-Salamin method is a strong general-purpose choice because of its quadratic convergence, stability, and relatively simple implementation. For very large-scale calculations, the Chudnovsky algorithm is often preferred. Its high digits-per-term rate, combined with binary splitting and fast multiplication, makes it one of the most efficient known approaches for record-scale \(\pi\) computation. A practical optimization strategy for large-scale \(\pi\) calculation is therefore: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item choose an algorithm with fast convergence; \item use arbitrary-precision arithmetic with guard digits; \item avoid naive factorial or term evaluation; \item use binary splitting for series-based methods; \item use fast multiplication and division algorithms; \item parallelize independent subproblems and arithmetic kernels where possible; \item manage memory carefully and checkpoint when necessary; \item validate the final result using independent checks and known digits. \end{enumerate} With these strategies, the computational cost of \(\pi\) calculation can be reduced from an intractable problem for classical methods to a manageable high-performance computing task for modern algorithms. \section{9. Applications and Extensions}\label{applications-and-extensions} \subsection{9.1 Numerical Analysis and Validation of High-Precision Arithmetic}\label{numerical-analysis-and-validation-of-high-precision-arithmetic} The numerical calculation of \(\pi\) is a natural test problem in numerical analysis because it combines a simple mathematical definition with a wide range of numerical behaviors: convergence, rounding, cancellation, guard-digit management, and the cost of high-precision arithmetic. As noted in \textbf{1. Introduction}, the required precision depends strongly on the application, ranging from a few digits in elementary engineering calculations to millions of digits in benchmarking, verification, and symbolic-numeric hybrid computation. A central use of \(\pi\) in numerical analysis is as a reference constant for validating algorithms and libraries. Because \(\pi\) is irrational and transcendental, it cannot be represented exactly in finite decimal form, so every computed value is an approximation. This makes \(\pi\) useful for checking whether a numerical method produces the expected number of correct digits, whether its error bounds are respected, and whether its stopping criterion is reliable. For example, the integral representation \(\pi = 4\int_0^1 \frac{dx}{1+x^2}\) provides a simple quadrature test: a numerical integration routine can be checked by comparing its result for \(\pi/4\) against known digits of \(\pi\). Similarly, the identity \(\arctan(1)=\frac{\pi}{4}\) connects \(\pi\) to inverse trigonometric functions and makes \(\pi\) a useful test case for arctangent series, Machin-like formulas, and related special-function implementations. High-precision libraries also benefit from \(\pi\) computation as a validation workload. As discussed in \textbf{6. Algorithm Implementation}, high-precision computation requires arbitrary-precision integers, rationals, or multiprecision floating-point types, together with guard digits above the requested output precision. Computing \(\pi\) to a specified number of digits exercises these components in a controlled setting. A library can be tested by requiring, for example, 100, 1000, or 100000 correct decimal digits and then comparing the result against known digits of \(\pi\). This type of test can reveal problems in rounding, overflow, underflow, term recurrence, division, square-root computation, and final normalization. The convergence and error analysis of \(\pi\) approximation methods also has direct numerical-analysis value. As summarized in \textbf{7. Convergence and Error Analysis}, different methods have very different convergence rates: \begin{itemize} \tightlist \item Archimedes' polygon method converges as \(O(n^{-2})\) in the number of sides \(n\). \item The Gregory-Leibniz series converges only as \(O(N^{-1})\) in the number of terms \(N\). \item Machin-like formulas converge exponentially, with Machin's formula giving roughly \(1.4\) decimal digits per term. \item The Brent-Salamin / AGM method has quadratic convergence, roughly doubling the number of correct digits each iteration. \item The Chudnovsky series gives about \(14.18\) decimal digits per term. \end{itemize} These differences make \(\pi\) a useful object for comparing numerical methods under a common target precision. A numerical analyst can study how many operations are required, how large the intermediate integers become, how sensitive the result is to rounding, and how the error bound behaves in practice. Rigorous or interval-based computation is another important application. The polygonal bounds \(n\sin\left(\frac{\pi}{n}\right) < \pi < n\tan\left(\frac{\pi}{n}\right)\) provide explicit intervals containing \(\pi\). Although these bounds are not efficient for very high precision, they illustrate a broader principle: numerical computation can produce not only an approximate value but also a certified interval. Similar ideas apply to series-based methods, where remainders can be bounded. For the Gregory-Leibniz series, for example, the alternating-series bound \(|\pi - 4S_N| \le \frac{4}{2N+3}\) gives a simple rigorous error estimate. For the Chudnovsky series, the first omitted term in the \(1/\pi\) series, scaled appropriately, gives a practical error estimate. These techniques are directly relevant to verified numerical computation, where the goal is not merely to obtain a plausible decimal expansion but to prove that the result lies within a specified interval. Finally, \(\pi\) appears in many standard numerical formulas, including Gaussian integrals, Fourier transforms, probability densities, and special functions. The Gaussian integral \(\int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi}\) connects \(\pi\) to probability, statistics, and numerical quadrature. The normalization constant \(1/\sqrt{2\pi}\) appears in the standard normal density, and factors of \(2\pi\) appear in Fourier analysis. Thus, accurate computation of \(\pi\) supports the validation of a broad class of numerical software, not only programs whose explicit purpose is to compute \(\pi\). \subsection{9.2 Cryptography and Secure Computation}\label{cryptography-and-secure-computation} The role of \(\pi\) in cryptography is indirect but practically important. \(\pi\) itself is not a cryptographic primitive, and its digits should not be used as secret keys, initialization vectors, or random numbers. Although the decimal digits of \(\pi\) are widely believed to behave pseudo-randomly, this property is not proven, and cryptographic applications require rigorously justified randomness. Nevertheless, \(\pi\) computation is useful in cryptographic software development because it exercises the same high-precision arithmetic infrastructure that appears in many cryptographic systems. Modern cryptographic protocols rely heavily on large-integer arithmetic. RSA key generation, modular exponentiation, elliptic-curve point multiplication, hash-function implementations, and many lattice-based or post-quantum schemes all require efficient multiplication, division, reduction, and sometimes square-root or related operations. As discussed in \textbf{8. Performance and Optimization}, the dominant cost in high-precision \(\pi\) computation is often large-integer multiplication, division, or square-root computation rather than the number of loop iterations alone. Therefore, a \(\pi\) computation can serve as a deterministic stress test for arbitrary-precision arithmetic libraries that may also be used in cryptographic software. A practical cryptographic application of \(\pi\) computation is the generation of deterministic test vectors. Because the digits of \(\pi\) are known to very high precision, they provide a reproducible reference for regression testing. A cryptographic library that uses arbitrary-precision integers can be tested by computing \(\pi\) to a specified number of digits and comparing the result against a known value. This does not test the cryptographic security of the library, but it can detect implementation errors in the underlying arithmetic kernels. Such errors are important because a bug in multiplication, division, or modular reduction can compromise cryptographic correctness even if the high-level protocol is sound. \(\pi\) computation can also be used to benchmark the performance of arithmetic kernels that are relevant to cryptography. For example, fast multiplication algorithms such as Karatsuba, Toom-Cook, FFT-based multiplication, or NTT-based multiplication may be evaluated by measuring how quickly they support a large \(\pi\) computation. In lattice-based cryptography, NTT-based multiplication is especially relevant, and a \(\pi\) benchmark can help compare the performance of different polynomial or integer multiplication backends. However, a \(\pi\) benchmark measures arithmetic throughput, not side-channel resistance, constant-time behavior, or protocol security. In secure computation, high-precision arithmetic may appear in homomorphic encryption, secure multiparty computation, or privacy-preserving numerical evaluation. In such settings, the ability to compute constants such as \(\pi\) accurately can be useful for normalization, calibration, or verification. Again, the value of \(\pi\) is not that it is a cryptographic constant, but that it provides a well-understood, deterministic, and scalable numerical workload. A useful distinction is therefore: \begin{itemize} \tightlist \item \textbf{Not appropriate:} using \(\pi\) digits as cryptographic randomness or key material. \item \textbf{Appropriate:} using \(\pi\) computation to test, validate, and benchmark arbitrary-precision arithmetic libraries that support cryptographic protocols. \item \textbf{Appropriate:} using known digits of \(\pi\) as deterministic test vectors for regression testing. \item \textbf{Appropriate:} using \(\pi\) computation to profile large-integer multiplication, division, and square-root routines that may also be used in cryptographic software. \end{itemize} This indirect role is important because cryptographic software often depends on general-purpose arbitrary-precision libraries. A library that is correct and efficient for \(\pi\) computation is more likely to be reliable for other large-integer workloads, although cryptographic use still requires additional security-specific validation. \subsection{9.3 Benchmarking and Performance Evaluation}\label{benchmarking-and-performance-evaluation} One of the most visible applications of \(\pi\) computation is benchmarking. Because \(\pi\) has a simple definition, known digits to very high precision, and a scalable precision parameter, it is a canonical benchmark for numerical software, arbitrary-precision libraries, and high-performance computing systems. A \(\pi\) benchmark can measure several quantities: \begin{itemize} \tightlist \item time to compute \(D\) decimal digits; \item memory usage as a function of working precision; \item throughput in digits per second; \item energy consumption; \item parallel scaling efficiency; \item performance of multiplication, division, and square-root algorithms; \item correctness of high-precision arithmetic. \end{itemize} As discussed in \textbf{8. Performance and Optimization}, the performance of large-scale \(\pi\) calculations depends on both the convergence rate of the algorithm and the cost of high-precision arithmetic. This makes \(\pi\) a useful benchmark for comparing algorithms at different levels. At the algorithmic level, one can compare: \begin{itemize} \tightlist \item the Gregory-Leibniz series, which is simple but inefficient; \item Machin-like formulas, which are faster and suitable for low to moderate precision; \item the Brent-Salamin / AGM method, which has quadratic convergence and is attractive for moderate to high precision; \item the Chudnovsky algorithm, which is especially effective for very high-precision and record-scale calculations. \end{itemize} At the arithmetic level, one can compare: \begin{itemize} \tightlist \item schoolbook multiplication, with cost \(O(p^2)\); \item Karatsuba multiplication; \item FFT-based or NTT-based multiplication, which is asymptotically more efficient for very large operands; \item different division and square-root algorithms; \item different memory layouts and cache behavior. \end{itemize} The Brent-Salamin / AGM method is particularly useful as a general-purpose benchmark because it requires only \(O(\log D)\) iterations for \(D\) digits and uses relatively few large arithmetic operations. Its iterations are sequential, so it is less amenable to algorithmic parallelization than binary-splitting methods, but its internal arithmetic can still be parallelized. The Chudnovsky algorithm, by contrast, is well suited to parallel and distributed computation because its binary-splitting structure allows independent subproblems to be evaluated separately. This makes it a strong benchmark for parallel systems, distributed memory architectures, and high-performance arbitrary-precision libraries. A well-designed \(\pi\) benchmark should specify: \begin{itemize} \tightlist \item the target number of decimal digits; \item the number of guard digits used; \item the algorithm, such as Machin-like, Brent-Salamin, or Chudnovsky; \item the multiplication algorithm, such as schoolbook, Karatsuba, FFT, or NTT; \item the precision model, such as arbitrary-precision integers or multiprecision floating-point; \item the output format, such as decimal digits or hexadecimal digits; \item the validation method, such as comparison against known digits or independent recomputation. \end{itemize} Without these details, benchmark results can be misleading. For example, a computation that produces 100000 digits using a fast multiplication algorithm may be much faster than one using schoolbook multiplication, even if the same \(\pi\) formula is used. Similarly, a method that uses many guard digits may be more accurate but slower and more memory-intensive. \(\pi\) computation is also useful for hardware validation. It can be used to test CPUs, GPUs, FPGAs, distributed clusters, and specialized accelerators. Because the workload is deterministic and the expected result is known, it is easier to validate than many scientific simulations. A system that computes the first \(N\) digits of \(\pi\) correctly can be checked against published values, providing a clear pass/fail criterion. Record-scale \(\pi\) computations have historically served as public benchmarks for computational power. They demonstrate the ability of a system to perform sustained high-precision arithmetic, manage large memory footprints, and coordinate parallel tasks. As discussed in \textbf{5. Modern Computational Algorithms}, the Chudnovsky algorithm is often preferred for such record-scale calculations when combined with binary splitting and fast multiplication. This combination allows the computation to scale to millions or billions of digits while keeping the arithmetic cost manageable. In addition to raw performance, \(\pi\) benchmarks can evaluate software engineering properties such as: \begin{itemize} \tightlist \item reproducibility; \item checkpointing and restart; \item memory management; \item error handling; \item validation against known digits; \item portability across platforms. \end{itemize} As noted in \textbf{8. Performance and Optimization}, practical optimization strategies include using arbitrary-precision libraries, avoiding naive factorial evaluation, using term recurrences, applying binary splitting, choosing fast multiplication algorithms, managing guard digits, and checkpointing long computations. A \(\pi\) benchmark can therefore serve not only as a performance test but also as a software-quality test. \subsection{9.4 Scientific Computing and Engineering}\label{scientific-computing-and-engineering} The constant \(\pi\) appears throughout scientific computing and engineering because it is tied to circles, spheres, waves, probability, and many physical laws. As discussed in \textbf{3. Mathematical Foundations}, \(\pi\) is fundamentally defined by circle geometry as the ratio of a circle's circumference to its diameter, and for the unit circle it equals both the area and half the circumference. This geometric origin makes \(\pi\) unavoidable in any computation involving circular or spherical objects. In engineering, \(\pi\) appears in formulas for area, volume, surface area, and moment of inertia. For example: \begin{itemize} \tightlist \item area of a circle: \(A = \pi r^2\); \item circumference of a circle: \(C = 2\pi r\); \item volume of a sphere: \(V = \frac{4}{3}\pi r^3\); \item surface area of a sphere: \(S = 4\pi r^2\); \item volume of a cylinder: \(V = \pi r^2 h\). \end{itemize} These formulas are used in mechanical design, fluid mechanics, structural engineering, thermal analysis, and manufacturing. In many such applications, standard double-precision floating-point arithmetic is sufficient, because the required accuracy is limited by measurement error, material variability, or model uncertainty. However, in metrology, high-precision simulation, or ill-conditioned numerical problems, higher precision may be needed. In physics, \(\pi\) appears in wave phenomena, angular measurements, and normalization constants. Common examples include: \begin{itemize} \tightlist \item angular frequency: \(\omega = 2\pi f\); \item wave number: \(k = 2\pi/\lambda\); \item phase: \(\phi = 2\pi x/\lambda\); \item Fourier series and Fourier transforms; \item normalization of wavefunctions in quantum mechanics; \item Gaussian integrals in statistical mechanics and quantum field theory. \end{itemize} The Gaussian integral \(\int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi}\) is especially important because it connects \(\pi\) to probability, statistics, and many areas of mathematical physics. The standard normal density \(\frac{1}{\sqrt{2\pi}}e^{-x^2/2}\) appears in error analysis, confidence intervals, signal processing, and machine learning. Accurate values of \(\pi\) and related constants are therefore important for correctly normalizing probability distributions and evaluating integrals. In numerical methods, \(\pi\) appears in spectral methods, finite-element methods, boundary integral methods, and PDE solvers. For example, problems on circular or spherical domains require trigonometric functions, polar or spherical coordinates, and Fourier expansions. Heat conduction in cylindrical coordinates, fluid flow in pipes, electromagnetic waveguides, and acoustic resonators all involve \(\pi\) in their analytical or numerical formulations. \(\pi\) is also used in Monte Carlo methods, both as a test problem and as a physical constant. The classical Monte Carlo estimate of \(\pi\) by random sampling in a unit square is a standard pedagogical example. It illustrates the basic idea of stochastic integration, but it converges slowly, typically as \(O(N^{-1/2})\), and is not suitable for high-precision computation. Its value is educational rather than practical for large-scale numerical work. In scientific computing, the choice of precision is usually application-driven. For most engineering simulations, double precision provides about 15-16 decimal digits of accuracy, which is often more than enough. However, high-precision computation may be needed when: \begin{itemize} \tightlist \item the problem is ill-conditioned; \item cancellation is severe; \item the computation is part of a verification or validation study; \item the result is used to generate reference data; \item the simulation involves long-time integration or sensitive dependence on parameters; \item the computation is used to test numerical libraries or hardware. \end{itemize} In such cases, the methods discussed in \textbf{5. Modern Computational Algorithms} and \textbf{6. Algorithm Implementation} become relevant. A scientific code may not need to compute \(\pi\) to millions of digits, but it may need a reliable high-precision value of \(\pi\) or a high-precision library that can compute it. The same infrastructure used for \(\pi\) computation - arbitrary-precision arithmetic, guard digits, stable recurrences, and validated stopping criteria - can be applied to other scientific constants and special functions. \subsection{9.5 Extensions to Other Mathematical Constants}\label{extensions-to-other-mathematical-constants} The techniques developed for computing \(\pi\) extend naturally to many other mathematical constants. The core ideas are the same: represent the constant by a rapidly convergent series, product, integral, or iterative method; analyze the convergence; control rounding error with guard digits; use efficient high-precision arithmetic; and validate the result against known values. A useful way to view these extensions is to separate the mathematical representation from the computational infrastructure. The representation determines the convergence rate and error bounds, while the infrastructure determines the practical cost. As discussed in \textbf{7. Convergence and Error Analysis} and \textbf{8. Performance and Optimization}, the dominant cost in high-precision computation is often large-integer multiplication, division, or square-root computation. Therefore, a constant that can be computed with a rapidly convergent series and efficient binary splitting may be far more practical than one requiring a slowly convergent series, even if the latter is mathematically simpler. Several important constants can be computed using methods closely related to those used for \(\pi\). \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Constant \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Common representation \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Practical high-precision method \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Relation to \(\pi\) computation \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot \(e\) & \(\displaystyle e = \sum_{n=0}^{\infty} \frac{1}{n!}\) & Binary splitting, fast factorial evaluation & Similar to series-based \(\pi\) methods; requires stable term recurrences \\ \(\sqrt{2}\) & Newton iteration \(x_{n+1} = \frac{1}{2}\left(x_n + \frac{2}{x_n}\right)\) & Quadratic iteration & Analogous to the quadratic convergence of the AGM method \\ \(\log 2\) & \(\displaystyle \log 2 = 2\operatorname{artanh}\left(\frac{1}{3}\right)\) & Accelerated arctangent-like series, binary splitting & Uses the same idea of small arguments for faster convergence \\ \(\gamma\) & Stirling expansion for the gamma function & Asymptotic expansion with error bounds & Requires careful error analysis, similar to series remainders for \(\pi\) \\ \(\zeta(2)\) & \(\displaystyle \zeta(2)=\frac{\pi^2}{6}\) & Compute from \(\pi\), or use Euler-Maclaurin summation & Directly related to \(\pi\) \\ \(\zeta(3)\) & Apéry's constant \(\displaystyle \zeta(3)=\sum_{n=1}^{\infty}\frac{1}{n^3}\) & Euler-Maclaurin summation or Apéry-type rapidly convergent series & Uses high-precision summation and error control \\ Catalan's constant \(G\) & \(\displaystyle G=\sum_{n=0}^{\infty}\frac{(-1)^n}{(2n+1)^2}\) & Accelerated series or binary splitting & Related to inverse trigonometric and Dirichlet beta functions \\ Elliptic integrals & Complete elliptic integral \(K(k)\) & AGM method & The AGM method is closely connected to \(\pi\) and elliptic integrals \\ \end{longtable} The exponential constant \(e\) is a particularly close analogue of \(\pi\) computation. The series \(e = \sum_{n=0}^{\infty} \frac{1}{n!}\) is simple, but naive factorial evaluation is inefficient for high precision. As with the Chudnovsky algorithm, production implementations should avoid naive factorial evaluation and instead use term recurrences, binary splitting, and fast multiplication. The error after \(N\) terms can be bounded using the remainder of the exponential series, giving a reliable stopping criterion. The constant \(\sqrt{2}\) can be computed by Newton's method: \(x_{n+1} = \frac{1}{2}\left(x_n + \frac{2}{x_n}\right).\) This iteration has quadratic convergence, similar in spirit to the Brent-Salamin / AGM method for \(\pi\). It is simple, stable, and well suited to high-precision computation. The main practical issues are the same as in \(\pi\) computation: choosing an initial approximation, controlling rounding error, and deciding when the desired number of digits has been reached. The constant \(\log 2\) can be computed using the slowly convergent series \(\log 2 = \sum_{n=1}^{\infty} \frac{(-1)^{n+1}}{n},\) but this is inefficient for high precision. A better approach uses an accelerated form such as \(\log 2 = 2\operatorname{artanh}\left(\frac{1}{3}\right) = 2\sum_{n=0}^{\infty} \frac{1}{(2n+1)3^{2n+1}}.\) This is analogous to the use of Machin-like formulas for \(\pi\), where smaller arguments lead to faster convergence. The Euler-Mascheroni constant \(\gamma\) is more challenging because it does not have a simple rapidly convergent series as straightforward as the Chudnovsky series for \(1/\pi\). High-precision computation often uses the Stirling expansion for the gamma function, together with careful error bounds. This illustrates an important point: extending \(\pi\) computation techniques to other constants requires new mathematical analysis, not merely copying the same code. The Riemann zeta values provide another important extension. Since \(\zeta(2) = \frac{\pi^2}{6},\) a high-precision value of \(\pi\) immediately gives a high-precision value of \(\zeta(2)\). Conversely, \(\zeta(2)\) can be computed independently using Euler-Maclaurin summation or other accelerated methods. For \(\zeta(3)\), Apéry's constant, rapidly convergent Apéry-type series exist, and high-precision computation again relies on stable recurrences, binary splitting, and fast multiplication. The AGM method used for \(\pi\) also extends to elliptic integrals. The complete elliptic integral of the first kind satisfies \(K(k) = \frac{\pi}{2\,\operatorname{AGM}(1,\sqrt{1-k^2})}.\) Thus, the same iterative machinery that computes \(\pi\) can be used to compute elliptic integrals, which appear in geometry, potential theory, and mathematical physics. This is a direct extension of the Brent-Salamin / AGM method from a single constant to a family of related functions. BBP-type formulas provide another extension. The BBP formula for \(\pi\), \(\pi = \sum_{k=0}^{\infty} \frac{1}{16^k} \left( \frac{4}{8k+1} - \frac{2}{8k+4} - \frac{1}{8k+5} - \frac{1}{8k+6} \right),\) allows extraction of base-16 digits of \(\pi\) without necessarily computing all preceding digits. Related digit-extraction ideas have been developed for other constants, including variants for \(\log 2\) and certain zeta values. These formulas are important because they change the computational problem from ``compute all digits up to \(N\)'' to ``extract a block of digits near position \(N\).'' The main lesson from these extensions is that \(\pi\) computation is not an isolated problem. It is a representative case of high-precision constant computation. The same principles apply: \begin{itemize} \tightlist \item choose a rapidly convergent representation; \item derive or use reliable error bounds; \item use guard digits; \item avoid unstable or inefficient recurrences; \item use fast multiplication and division; \item apply binary splitting when the computation is a large sum; \item parallelize independent subproblems when possible; \item validate the final result against known digits or independent methods. \end{itemize} \subsection{9.6 Practical Guidance}\label{practical-guidance} The appropriate use of \(\pi\) computation depends on the application. The following table summarizes common scenarios and practical recommendations. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2308}} >{\raggedleft\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.3077}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2308}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2308}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Application \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedleft Typical precision \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Suggested approach \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Validation \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Elementary geometry or engineering & 6-15 digits & Standard floating-point constant or Machin-like formula & Compare with standard library value \\ Numerical analysis teaching & 10-100 digits & Gregory-Leibniz for simplicity, Machin-like for efficiency & Known digits, error bounds \\ High-precision library testing & 100-100000 digits & Brent-Salamin / AGM or Chudnovsky & Known digits, independent recomputation \\ Benchmarking arbitrary-precision arithmetic & 10000-millions of digits & Chudnovsky with binary splitting and fast multiplication & Known digits, checksums, independent runs \\ Parallel or distributed benchmarking & Large scale & Chudnovsky with binary splitting & Known digits, task-level validation \\ Scientific simulation & Usually double precision & Use a trusted constant; high precision only if needed & Compare with reference values \\ Cryptographic software testing & Arbitrary precision & Use \(\pi\) as deterministic test vector for big-integer arithmetic & Known digits, regression tests \\ Extension to other constants & Problem-dependent & Use analogous series, AGM, binary splitting, or asymptotic methods & Known digits, error bounds, independent methods \\ \end{longtable} For low to moderate precision, Machin-like formulas or the Brent-Salamin / AGM method are often sufficient. As discussed in \textbf{7. Convergence and Error Analysis} and \textbf{8. Performance and Optimization}, Machin-like formulas are simple, stable, and suitable for low to moderate precision, while Brent-Salamin is a strong general-purpose choice because of its quadratic convergence and implementation simplicity. For very high-precision or record-scale calculations, the Chudnovsky algorithm is usually preferred when combined with binary splitting and fast multiplication. In all cases, the final result should be validated. As noted in \textbf{6. Algorithm Implementation} and \textbf{8. Performance and Optimization}, validation can be done by comparing against known digits of \(\pi\), by independent recomputation, or by using a different algorithm. Guard digits should be used to reduce accumulated rounding error, and stopping criteria should be matched to the convergence behavior of the chosen method. The broader significance of \(\pi\) computation is that it provides a common testbed for numerical algorithms, high-precision libraries, performance optimization, and scientific validation. It connects classical mathematics with modern computing, and the techniques developed for \(\pi\) extend naturally to many other constants and functions used in science, engineering, and secure computation. \section{10. Conclusion}\label{conclusion} \subsection{10.1 Summary of Main Findings}\label{summary-of-main-findings} The numerical calculation of \(\pi\) illustrates a broad progression from simple geometric reasoning to highly optimized high-precision computation. As introduced in \textbf{1. Introduction}, \(\pi\) is a fundamental constant whose irrational and transcendental nature makes exact finite representation impossible, so numerical approximation is essential in mathematics, science, and computing. The required precision varies widely: a few digits may suffice for elementary geometry, while scientific computing, benchmarking, and record-scale calculations may require millions or more digits. The historical development reviewed in \textbf{2. Historical Background} shows that early methods were geometric, using inscribed and circumscribed polygons to bound \(\pi\). Archimedes' 96-gon calculation gave the classical bounds \(\frac{223}{71} < \pi < \frac{22}{7}.\) Later, the development of calculus and infinite series shifted the problem from geometry to analysis. The Gregory-Leibniz series, \(\pi = 4\left(1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots\right),\) is mathematically simple but converges slowly. Machin-like formulas improved practical computation by using arctangent identities with small arguments, and modern algorithms such as Brent-Salamin, Chudnovsky, and BBP made extremely high-precision computation feasible. The mathematical foundations presented in \textbf{3. Mathematical Foundations} connect \(\pi\) to circle geometry, trigonometric identities, infinite series, and integral representations. In particular, the identity \(\pi = 4\arctan(1)\) and the arctangent addition formula provide the basis for many efficient series-based methods. Integral representations, such as \(\pi = 4\int_0^1 \frac{dx}{1+x^2},\) also connect \(\pi\) to numerical quadrature and analysis. The classical methods discussed in \textbf{4. Classical Numerical Methods} are historically important and easy to understand, but they are generally inefficient for high precision. Archimedes' polygon method converges as \(O(n^{-2})\) in the number of sides \(n\), while the Gregory-Leibniz series converges only as \(O(N^{-1})\) in the number of terms \(N\). Their main value is pedagogical and conceptual rather than computational. The modern algorithms examined in \textbf{5. Modern Computational Algorithms} are far more efficient. Machin-like formulas converge exponentially and are suitable for low to moderate precision. The Brent-Salamin / AGM method has quadratic convergence, roughly doubling the number of correct digits at each iteration, and is attractive because it is simple and stable. The Chudnovsky algorithm is especially powerful for very high precision, yielding about 14 decimal digits per term, and is the preferred method for record-scale calculations when combined with binary splitting and fast multiplication. The implementation discussion in \textbf{6. Algorithm Implementation} emphasizes that the choice of numeric type is critical. Native floating-point types are adequate only for low precision, while high-precision computation requires arbitrary-precision integers, rationals, or multiprecision floating-point libraries. Guard digits, stable recurrences, careful stopping criteria, and validation against known digits are all essential for reliable results. The convergence and error analysis in \textbf{7. Convergence and Error Analysis} shows that method choice depends strongly on the desired precision and the available arithmetic environment. Classical methods have simple error bounds but poor efficiency. Machin-like formulas are stable and practical for moderate precision. Brent-Salamin is a strong general-purpose method, while Chudnovsky is best for very large computations when implemented with binary splitting, fast multiplication, and careful high-precision division. The performance discussion in \textbf{8. Performance and Optimization} highlights that, for large-scale calculations, the dominant cost is usually high-precision arithmetic rather than the number of algorithmic iterations alone. Multiplication, division, and square-root computation dominate the runtime, and the choice of multiplication algorithm - schoolbook, Karatsuba, or FFT/NTT-based methods - has a major effect on performance. Memory usage is determined primarily by the working precision, and parallelization is especially natural for Chudnovsky-based binary splitting. Finally, \textbf{9. Applications and Extensions} shows that \(\pi\) computation is not merely a historical curiosity. It is a useful test problem for numerical analysis, a benchmark for arbitrary-precision libraries, a tool for validating numerical software, and a gateway to the computation of other mathematical constants. The same principles - rapid convergence, reliable error bounds, guard digits, fast arithmetic, and validation - apply broadly to high-precision constant computation. \subsection{10.2 Comparative Strengths and Limitations}\label{comparative-strengths-and-limitations} The following table summarizes the main strengths and limitations of the methods discussed in the publication. \begin{longtable}[]{@{} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}} >{\raggedright\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.2500}}@{}} \toprule\noalign{} \begin{minipage}[b]{\linewidth}\raggedright Method \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Main Strengths \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Main Limitations \end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright Best Use \end{minipage} \\ \midrule\noalign{} \endhead \bottomrule\noalign{} \endlastfoot Archimedes' polygon method & Geometrically intuitive; gives rigorous bounds; historically important & Slow convergence, \(O(n^{-2})\); impractical for high precision & Educational and historical context \\ Gregory-Leibniz series & Very simple to state and implement; easy error bound & Very slow convergence, \(O(N^{-1})\); many terms required; inefficient for high precision & Low-precision demonstrations; pedagogical use \\ Machin-like formulas & Much faster than Gregory-Leibniz; relatively simple; stable with term recurrences & Less efficient than the best modern algorithms for very high precision & Low to moderate precision \\ Brent-Salamin / AGM & Quadratic convergence; simple; stable; good general-purpose method & Iterations are sequential, limiting algorithmic parallelism; still requires high-precision arithmetic & Moderate to high precision; general-purpose implementation \\ Chudnovsky algorithm & About 14 decimal digits per term; excellent for record-scale computation & More complex to implement; requires binary splitting, fast multiplication, and careful memory management & Very high precision; large-scale and distributed computation \\ \end{longtable} Several important contrasts emerge. First, simplicity and efficiency are often in tension. The Gregory-Leibniz series is one of the easiest methods to implement, but its slow convergence makes it unsuitable for serious high-precision work. Archimedes' polygon method is similarly simple, but increasing the number of polygon sides quickly becomes inefficient. Second, convergence rate is not the only factor. A method may converge quickly in theory but still be difficult to implement efficiently. The Chudnovsky algorithm is extremely efficient in terms of digits per term, but a production implementation must avoid naive factorial evaluation and must use term recurrences, binary splitting, and fast multiplication. Without these optimizations, its theoretical advantage may not be realized in practice. Third, numerical stability depends on implementation. The Brent-Salamin method is generally stable, although the subtraction \(c_n = a_n - b_n\) can lose relative precision if not handled carefully. The Chudnovsky algorithm is stable in exact arithmetic, but practical high-precision computation requires careful control of rounding, guard digits, and division. Fourth, performance is dominated by arithmetic cost. For small precision, the number of iterations may matter more. For large precision, the cost of multiplying, dividing, and taking square roots of large integers becomes the main bottleneck. Thus, the best algorithm is not always the one with the fastest convergence; it is the one whose convergence and arithmetic requirements match the implementation environment. \subsection{10.3 Practical Recommendations}\label{practical-recommendations-1} A practical choice of method depends on the target precision, the available software environment, and the desired balance between simplicity and performance. For low precision, standard floating-point arithmetic is usually sufficient. In this regime, simple formulas such as Machin-like identities or even the Gregory-Leibniz series may be acceptable for educational purposes. However, even for modest precision, Machin-like formulas are generally preferable to Gregory-Leibniz because they converge much faster. For moderate precision, the Brent-Salamin / AGM method is a strong general-purpose choice. It has quadratic convergence, is relatively simple to implement, and is stable when implemented with appropriate guard digits. It is often a good default for libraries that need a reliable method without the full complexity of record-scale Chudnovsky implementations. For very high precision, the Chudnovsky algorithm is usually preferred, especially when combined with binary splitting and fast multiplication. It produces many digits per term and is well suited to large-scale computation. However, it requires more sophisticated implementation techniques, including arbitrary-precision arithmetic, careful memory management, and efficient large-integer multiplication. In all cases, the following implementation practices are important: \begin{itemize} \tightlist \item Use arbitrary-precision arithmetic when the target precision exceeds the reliable range of native floating-point types. \item Include guard digits above the requested output precision to reduce accumulated rounding error. \item Choose stopping criteria that match the convergence behavior of the method. \item Use stable term recurrences rather than naive evaluation of factorials or large powers. \item Validate the final result against known digits of \(\pi\) or by independent recomputation. \item For long computations, use checkpointing to avoid losing progress after failures. \item For large-scale calculations, use fast multiplication algorithms and, where possible, parallelize independent subproblems. \end{itemize} A reasonable practical hierarchy is therefore: \begin{enumerate} \def\labelenumi{\arabic{enumi}.} \tightlist \item \textbf{Low precision:} standard floating-point with simple formulas or Machin-like identities.\\ \item \textbf{Moderate precision:} Brent-Salamin / AGM or well-implemented Machin-like formulas.\\ \item \textbf{Very high precision:} Chudnovsky with binary splitting, fast multiplication, and careful high-precision arithmetic. \end{enumerate} \subsection{10.4 Future Research and Implementation Improvements}\label{future-research-and-implementation-improvements} Although the numerical computation of \(\pi\) is a well-studied problem, several directions remain important for future research and implementation. One major direction is the continued improvement of high-precision arithmetic. Since large-scale \(\pi\) computation is dominated by multiplication, division, and square-root computation, advances in fast multiplication algorithms can have a direct impact on performance. This includes further development of FFT- and NTT-based methods, cache-aware implementations, and algorithms that exploit modern memory hierarchies. Another important area is parallel and distributed computation. The Chudnovsky algorithm, when implemented with binary splitting, naturally decomposes into independent subproblems. Future work can improve load balancing, communication overhead, and checkpointing for distributed systems. GPU and accelerator-based implementations of large-integer arithmetic may also provide significant speedups for very large calculations. Memory efficiency is another practical concern. High-precision computations require large working precision, and binary-splitting implementations can generate substantial intermediate data. More memory-efficient recurrences, streaming evaluation techniques, and better integration with disk-based or distributed storage could make record-scale computations more accessible. Certified and formally verified computation is also an important direction. Because \(\pi\) has known reference digits, it is an excellent test case for verifying arbitrary-precision libraries and numerical algorithms. Future work could combine high-precision computation with formal verification, interval arithmetic, or rigorous error certification to provide stronger guarantees about correctness. Adaptive precision management is another useful improvement. Rather than requiring users to specify a fixed working precision, software could estimate the necessary guard digits automatically based on the algorithm, target precision, and observed error behavior. This would make high-constant computation easier to use and less error-prone. Finally, the techniques developed for \(\pi\) computation extend naturally to other mathematical constants, including \(e\), \(\sqrt{2}\), \(\log 2\), \(\gamma\), zeta values, Catalan's constant, and elliptic integrals. A broader research direction is the development of general-purpose frameworks for high-precision constant computation that combine rapidly convergent representations, reliable error bounds, fast arithmetic, and automatic validation. In summary, the numerical calculation of \(\pi\) remains a valuable benchmark and a rich area for algorithmic development. Classical methods provide historical and conceptual insight, while modern algorithms make extremely high-precision computation practical. Future improvements will likely come not only from new mathematical formulas, but also from better arithmetic, parallelization, memory management, and certified computation. \end{document}