Arbitrary Precision Computation
Standard computer registers represent floating-point numbers using the IEEE 754 standard, allocating a fixed 64 bits (double precision) of memory. This limits representations to roughly 15-17 significant decimal digits. If you calculate $ \pi $ or any mathematical constant using standard variables, error accumulation quickly destroys accuracy at higher targets.
To calculate constants to millions of decimal places, we use Arbitrary Precision Arithmetic. This technique stores numbers as dynamic arrays of digital blocks (often 32-bit or 64-bit integers representing large bases like $ 2^{32} $ or $ 2^{64} $). Using advanced backend engines (like rug/GMP or dashu), we perform basic math operations (addition, multiplication, division) on these large arrays. This page explains the algorithms used to construct high-precision constants from these array-based backends.
Pi (π) — Chudnovsky Algorithm
The calculation of Pi in apcn uses the famous Chudnovsky Algorithm, published by the Chudnovsky brothers in 1988. It is a rapid Ramanujan-like hypergeometric series that yields approximately 14 decimal digits per term, making it the industry standard for world-record computations.
To implement this formula at extreme precisions efficiently, the hypergeometric series is rewritten into a fraction of integer terms and computed using Binary Splitting. This reduces the problem of dividing large floating-point numbers at every step to a single final high-precision division, parallelizing the integer multiplications across CPU cores.
Euler's Number (e) — Taylor Series
Euler's number $ e $ is calculated using its natural infinite Taylor series definition:
This series converges extremely fast because factorials grow super-exponentially. For a target of $ N $ digits of precision, the number of terms needed is estimated by solving for $ k $ in:
Using binary splitting, the sum of fractions is grouped recursively into a single numerator and denominator before performing one final division:
- Numerator P: Computed recursively via divide-and-conquer.
- Denominator Q: Represents the factorial products.
- Final Division: $ P / Q $ produces the high-precision floating-point constant.
Euler-Mascheroni Constant (γ) — Brent-McMillan B3
The Euler-Mascheroni constant $ \gamma $ is the limiting difference between the harmonic series and the natural logarithm. It is one of the most challenging constants to calculate to high precision.
apcn uses the highly optimized Brent-McMillan B3 Algorithm, which utilizes modified Bessel functions:
Where $ I_0 $ and $ K_0 $ are modified Bessel functions, and $ S_0 $ is a related series. The algorithm works by:
- Bessel Series: Computing the Bessel functions $ I_0 $ and $ S_0 $ using binary splitting for a chosen parameter $ N $.
- Asymptotic Expansion: Evaluating the asymptotic expansion of $ I_0(2N) K_0(2N) $ to avoid calculating the extremely slow $ K_0 $ series directly.
- Logarithm Correction: Subtracting $ \ln(N) $ using a high-precision logarithm calculator.
Natural Logarithms (ln) — AGM & Taylor
APCN computes logarithms using two complementary algorithms based on the magnitude of the input parameter:
- Logarithmic Taylor Series: For values near 1, we use the classic Mercator series combined with binary splitting for fast convergence:
$$\ln(1+x) = \sum_{n=1}^{\infty} \frac{(-1)^{n-1} x^n}{n} \quad \text{for } |x| < 1$$
- Arithmetic-Geometric Mean (AGM): For larger values, Taylor series converge too slowly. Instead, we use the AGM method, which displays quadratic convergence:
$$\ln(x) \approx \frac{\pi}{2 \cdot \text{AGM}\left(1, \frac{4}{2^m}\right)} - m \ln(2)$$where $ m $ is a large scaling factor. This allows calculating the natural logarithm of any number with high efficiency by reducing it to iterative arithmetic averages.
Golden Ratio (φ) — Direct Calculation
The Golden Ratio $ \phi $ has a simple algebraic definition:
Instead of computing an infinite series directly for $ \phi $, we calculate the square root of 5 ($ \sqrt{5} $) to the target precision using the Newton-Raphson iteration, add 1 (a cheap integer operation), and shift the division bits by 1 (equivalent to dividing by 2). This makes the performance of $ \phi $ closely match the computation speed of square roots.
Square Roots (√) — Newton-Raphson
Square roots are computed using the Newton-Raphson Method, which is equivalent to solving the function $ f(y) = y^2 - x = 0 $. The iterative step is defined as:
This method exhibits quadratic convergence, meaning the number of correct decimal digits doubles with every single iteration. If you start with a 15-digit double estimate, the precision grows as $ 15 \to 30 \to 60 \to 120 \to 240 \to \cdots $ meaning only 16 iterations are required to calculate 1,000,000 digits!
Binary Splitting & Parallelization
Binary Splitting is a divide-and-conquer optimization technique for evaluating infinite series with rational terms.
Consider a series where each term can be expressed as a fraction. If computed sequentially, each step requires a floating-point division and multiplication of large numbers, which is extremely expensive. Binary splitting recursively splits the terms of the series into two halves:
Instead of dividing floating-point numbers at each step, we merge the rational fractions together using integer additions and multiplications. We only perform a single high-precision division at the very end of the calculation.
Why Parallelization Works: Since the recursion tree divides the work into independent sub-trees, we can run the calculations of the left and right sub-trees in parallel across multiple CPU threads using Rayon. This speeds up computation on multi-core processors.