🔢 Computing Pi to a Million Digits
The Mathematics, the History, and the Code — A deep-dive from Archimedes to 314 trillion digits, with a fully production-ready Python implementation
Chapter One — Why Pi?
Of all the mathematical constants, none has captured human imagination more persistently than pi. It appears in the formula for the area of a circle, in the period of a pendulum, in the probability that two randomly chosen integers are coprime, in the normal distribution that underlies all of statistics, in the eigenvalues of the hydrogen atom, and in hundreds of other places where circles are nowhere to be seen. It is simultaneously the most elementary and the most profound constant in mathematics.
But this article is not really about pi. It is about the extraordinary human effort to compute pi — to pin down its decimal expansion to more and more places — and about what that effort has taught us about algorithms, arithmetic, and the deep structure of mathematics itself. The story of pi computation is, in miniature, the entire history of computational mathematics.
And at the end of the article, we will build a production-ready Python application that computes pi to any desired number of digits using the same algorithm family that set the world record of 314 trillion digits in December 2025.
Chapter Two — From Archimedes to the Infinite Series
The first person to compute pi rigorously was Archimedes of Syracuse, around 250 BC. His method was geometric and breathtakingly clever. He inscribed and circumscribed regular polygons around a circle of diameter 1, observing that the circle’s circumference is trapped between the perimeters of the two polygons. Starting with hexagons and repeatedly doubling the number of sides, he reached 96-gons and established:
which gives pi correct to two decimal places. This polygonal method was the state of the art for nearly 1,900 years. Zu Chongzhi in 5th-century China pushed it to seven decimal places using 12,288-sided polygons, a record that stood for 900 years.
The revolution came in the 17th century with the invention of calculus and the discovery of infinite series. James Gregory (1671) and Gottfried Leibniz (1674) independently discovered:
This is beautiful but useless for computation: it converges so slowly that you need 10 billion terms to get 10 decimal places. The real breakthrough came from Machin-like formulas. In 1706, John Machin discovered:
The arctangent series converges rapidly when its argument is small, so this formula is extremely efficient. Machin used it to compute 100 decimal places by hand. The general form — expressing \(\pi/4\) as a linear combination of arctangents of rational numbers — spawned an entire industry of Machin-like formulas, each optimised for a different balance of terms and convergence speed. William Shanks spent 15 years computing 707 digits by hand using such a formula in 1873 (though he made an error at digit 528 that went undetected for 70 years).
The theoretical reason these formulas work is the identity:
For \(x = 1/5\), each term is 25 times smaller than the previous one, giving roughly 1.4 decimal digits per term — already vastly better than the Leibniz formula.
Chapter Three — The Computer Era and the Race to Digits
The arrival of electronic computers in the 1950s transformed pi computation from a lifetime’s work into a matter of hours. ENIAC computed 2,037 digits in 1949. By 1973, Jean Guilloud and Martine Bouyer had reached 1 million digits on a CDC 7600. The algorithms were still Machin-like formulas, but now executed at electronic speed.
The next algorithmic revolution came from an unexpected direction: the theory of elliptic integrals and the arithmetic-geometric mean (AGM). In 1975, Richard Brent and Eugene Salamin independently discovered an algorithm based on the AGM that converges quadratically — the number of correct digits doubles with every iteration. Starting from:
and iterating:
the algorithm converges to:
with the number of correct digits doubling at each step. Starting from scratch, 50 iterations suffice for a quadrillion digits (in principle). This was a qualitative change: the Machin-like formulas converge linearly (fixed digits per term), while the AGM converges quadratically (digits per iteration grow exponentially).
The mathematical reason for this miracle is deep: the AGM is intimately connected to the theory of elliptic integrals, and pi appears as a special value of a complete elliptic integral. The quadratic convergence is a consequence of the modular equation theory developed by Gauss and Legendre in the early 19th century.
Chapter Four — Ramanujan and the Chudnovsky Brothers
While the Brent-Salamin algorithm was a theoretical triumph, it has a practical drawback: each iteration requires a high-precision square root, which is expensive. The algorithms that have actually set world records since the 1990s are based on a different source: the extraordinary formulas of Srinivasa Ramanujan.
In 1914, Ramanujan published a paper containing a series of formulas for \(1/\pi\) of the form:
Each term of this series adds approximately 8 decimal digits. Ramanujan had no proof; he derived the formulas from his extraordinary intuition about modular forms. The proofs came decades later, through the work of Jonathan and Peter Borwein in the 1980s, who showed that Ramanujan’s formulas are special cases of a general theory connecting pi to the theory of complex multiplication of elliptic curves.
In 1988, David and Gregory Chudnovsky discovered a variant that converges even faster:
Each term adds approximately 14.18 decimal digits. This is the fastest-converging series of its type, and it has been used in every world-record pi computation since the early 1990s.
The constant 640320 is not arbitrary. It is deeply connected to the largest Heegner number, 163. The nine Heegner numbers are the positive integers \(d\) for which the imaginary quadratic field \(\mathbb{Q}(\sqrt{-d})\) has class number 1. The largest is 163, and the \(j\)-invariant of the corresponding elliptic curve satisfies:
This is why the famous near-integer \(e^{\pi\sqrt{163}} \approx 262537412640768743.999999999999\ldots\) is so close to an integer: it is the exponential of the \(j\)-invariant argument, and the tiny deviation from an integer is the error in the approximation. The Chudnovsky formula exploits this deep number-theoretic coincidence to achieve its extraordinary convergence rate.
Chapter Five — Binary Splitting: Turning Series into Fast Multiplication
The Chudnovsky series converges at 14 digits per term, so computing 1 million digits requires about 70,000 terms. Computing each term individually at full precision would require 70,000 high-precision divisions, which is expensive. The key optimisation that makes the algorithm practical is binary splitting.
The idea is to represent the partial sum \(S(a, b) = \sum_{k=a}^{b-1} t_k\) as a ratio \(T(a,b) / Q(a,b)\) of two integers, and to compute this ratio by a divide-and-conquer recursion:
where \(m = \lfloor(a+b)/2\rfloor\) and \(P(a,b)\) tracks the product of numerator polynomial factors. This recursion is applied all the way down to single-term base cases, and then the results are combined back up the recursion tree. The crucial point is that only one division is performed at the very end, at full precision. All intermediate operations are integer multiplications.
Why does this matter? Because for large numbers, multiplication is much cheaper than division (division requires an iterative Newton’s method that costs roughly 2–3 multiplications), and more importantly, the binary splitting tree structure allows the use of fast multiplication algorithms (Karatsuba, Toom-Cook, Schönhage-Strassen, NTT) that run in sub-quadratic time.
The overall complexity of computing pi to \(n\) digits with binary splitting and fast multiplication is:
where \(M(n)\) is the cost of multiplying two \(n\)-digit numbers. With GMP’s implementation of Schönhage-Strassen:
giving a total complexity of:
This is essentially optimal for the problem, and it is why the Chudnovsky algorithm with binary splitting has dominated pi computation for three decades.
Chapter Six — The BBP Formula: Computing Individual Digits
In 1995, David Bailey, Peter Borwein, and Simon Plouffe discovered something that shocked the mathematical community: a formula that allows computing the \(n\)-th hexadecimal digit of pi without computing any of the preceding digits:
To extract the \(n\)-th hexadecimal digit, multiply both sides by \(16^n\) and take the fractional part. Each term in the resulting sum can be computed using modular exponentiation — the same technique used in RSA cryptography — which requires only ordinary integer arithmetic, not arbitrary-precision arithmetic. The total cost is \(\mathcal{O}(n \log n)\) bit operations, and the memory requirement is \(\mathcal{O}(\log n)\).
The BBP formula is not useful for computing pi to many digits (it is slower than Chudnovsky for that purpose), but it is invaluable as an independent verification tool: after computing pi to \(n\) digits with Chudnovsky, you can spot-check specific hexadecimal digits at arbitrary positions using BBP, providing high confidence that the computation is correct without recomputing everything.
Chapter Seven — The World Records: From Millions to Trillions
The table below summarises the major milestones in pi computation:
| Year | Digits | Who / What |
|---|---|---|
| 1949 | 2,037 | ENIAC |
| 1973 | 1,000,000 | Guilloud & Bouyer, CDC 7600 |
| 1989 | 1,011,196,691 | Chudnovsky brothers, custom supercomputer |
| 2002 | 1,241,100,000,000 | Kanada et al., Hitachi SR8000 |
| 2009 | 2,576,980,370,000 | Fabrice Bellard, desktop PC |
| 2019 | 31,415,926,535,897 | Emma Haruka Iwao, Google Cloud |
| 2021 | 62,831,853,071,796 | University of Applied Sciences Grisons |
| 2024 | 105,000,000,000,000 | StorageReview |
| May 2025 | 300,000,000,000,000 | KIOXIA & Linus Media Group |
| Dec 2025 | 314,159,265,358,979 ⭐ | StorageReview, Dell PowerEdge R7725 |
The December 2025 record of 314 trillion digits (a delightful approximation of pi itself, times \(10^{14}\)) was set on a single Dell PowerEdge R7725 server equipped with two AMD EPYC 192-core CPUs and 40 Micron 61.44 TB NVMe SSDs. The computation ran for 110 days without interruption, using Alexander Yee’s y-cruncher software, which implements the Chudnovsky algorithm with binary splitting, multi-threaded parallelism, AVX-512 SIMD vectorisation, and a custom I/O subsystem sustaining 280 GB/s of storage throughput.
Chapter Eight — The Theoretical Frontier
We have reached the frontier of what is currently known and implemented. But mathematics and computer science are living disciplines, and deep open questions remain.
The Harvey-Hoeven algorithm (2019) proves that integer multiplication can be done in \(\mathcal{O}(n \log n)\) time, which is theoretically optimal. However, the algorithm has such enormous constant factors that it is not competitive in practice for any feasible number of digits. The practical frontier for multiplication remains variants of Schönhage-Strassen and NTT-based algorithms.
The question of whether pi is a normal number remains completely open. A number is normal in base 10 if every digit 0–9 appears with frequency \(1/10\), every two-digit sequence with frequency \(1/100\), and so on. Statistical tests on the known digits of pi are consistent with normality, but no proof exists.
The BBP formula raises a fascinating question: is there a formula that allows computing the \(n\)-th decimal digit of pi without computing all preceding digits? No such formula is known, and there are theoretical reasons to suspect none exists, but this has not been proved.
From a complexity standpoint, computing pi to \(n\) digits requires at least \(\Omega(n \log n)\) bit operations (from information theory), and the best known algorithms achieve \(\mathcal{O}(n \log^3 n \log \log n)\). Closing this gap is an open problem.
The journey from Archimedes’ polygons to 314 trillion digits spans over two millennia of mathematical progress. Pi computation is a microcosm of the entire history of mathematics and computer science, and it continues to drive progress in both fields.
Chapter Nine — The Code
Now we build it. The application below is a fully production-ready Python package implementing:
- Chudnovsky algorithm with binary splitting
- gmpy2 / GMP backend for fast arbitrary-precision arithmetic (graceful fallback)
- BBP formula for independent spot-check verification
- Full CLI with
argparse - Structured logging throughout
- Type hints on every function
- Docstrings on every module, class, and function
- Guard digits to absorb rounding drift
- Timing reports for every phase
pytesttest suite
Installation
Project Structure
├── __init__.py
├── __main__.py
├── big_int.py
├── chudnovsky.py
├── calculator.py
├── verifier.py
├── cli.py
└── tests/
├── __init__.py
└── test_pi.py
pi_calculator/__init__.py
"""
pi_calculator
=============
A production-ready, high-performance pi computation package.
Implements the Chudnovsky algorithm with binary splitting, an integer
Newton square-root, BBP-formula verification, and a clean CLI.
Uses gmpy2 (backed by GMP / Schönhage-Strassen multiplication) for fast
arbitrary-precision arithmetic when available, falling back to Python's
native integers transparently.
Author : SiemensGPT reference implementation — June 2026
License : MIT
"""
from .calculator import PiCalculator
from .verifier import BbpVerifier
__all__ = ["PiCalculator", "BbpVerifier"]
__version__ = "1.0.0"pi_calculator/big_int.py
"""
big_int.py
==========
Arithmetic backend abstraction.
Provides a unified interface for arbitrary-precision integer arithmetic,
automatically using gmpy2 (backed by GMP, which implements Schönhage-
Strassen / Toom-Cook multiplication) when available, and falling back to
Python's built-in integers otherwise.
gmpy2 is typically 3–10× faster than Python integers for the sizes
encountered in a million-digit pi computation.
"""
from __future__ import annotations
import logging
logger = logging.getLogger(__name__)
# ---------------------------------------------------------------------------
# Attempt to import gmpy2; fall back gracefully if unavailable.
# ---------------------------------------------------------------------------
try:
import gmpy2 # type: ignore
def mpz(n: int) -> "gmpy2.mpz":
"""Wrap a Python int as a gmpy2 multi-precision integer."""
return gmpy2.mpz(n)
def isqrt(n: "gmpy2.mpz | int") -> "gmpy2.mpz":
"""
Integer square root via gmpy2.
Uses GMP's highly optimised mpz_sqrt, which internally applies
Newton's method with Schönhage-Strassen multiplication for large
operands.
Args:
n: A non-negative integer.
Returns:
floor(sqrt(n)) as a gmpy2.mpz.
Raises:
ValueError: If n is negative.
"""
if n < 0:
raise ValueError("isqrt is undefined for negative integers.")
return gmpy2.isqrt(gmpy2.mpz(n))
BACKEND: str = "gmpy2 (GMP — Schönhage-Strassen / Toom-Cook)"
except ImportError:
logger.warning(
"gmpy2 not found — falling back to Python built-in integers. "
"Install gmpy2 for significantly better performance: pip install gmpy2"
)
def mpz(n: int) -> int: # type: ignore[misc]
"""Identity when gmpy2 is unavailable."""
return int(n)
def isqrt(n: int) -> int: # type: ignore[misc]
"""
Pure-Python integer square root using Newton's method.
Converges quadratically: each iteration doubles the number of
correct bits. The starting estimate is a power of 2 obtained from
the bit-length of n, which is always within a factor of 2 of the
true answer, guaranteeing O(log log n) iterations.
Args:
n: A non-negative integer.
Returns:
floor(sqrt(n)) as a Python int.
Raises:
ValueError: If n is negative.
"""
if n < 0:
raise ValueError("isqrt is undefined for negative integers.")
if n == 0:
return 0
# Initial estimate: largest power of 2 not exceeding sqrt(n).
x = 1 << ((n.bit_length() + 1) >> 1)
while True:
x_next = (x + n // x) >> 1
if x_next >= x:
return x
x = x_next
BACKEND: str = "Python built-in integers (install gmpy2 for ~5× speedup)"pi_calculator/chudnovsky.py
"""
chudnovsky.py
=============
Chudnovsky algorithm with binary splitting.
The Chudnovsky formula (1988), derived from Ramanujan's 1914 series:
1/Ï€ = (12 / 640320^(3/2))
× Î£_{k=0}^{∞} (−1)^k (6k)! (13591409 + 545140134k)
─────────────────────────────────────
(3k)! (k!)³ 640320^(3k)
Convergence rate: ~14.1816 decimal digits per term.
Binary splitting transforms the partial sum S(0, N) into a single
rational P/Q by a divide-and-conquer recursion, enabling sub-quadratic
multiplication for the dominant arithmetic cost.
Complexity (with fast multiplication):
O(M(n) · log²(n))
where M(n) is the cost of multiplying two n-digit numbers.
With GMP's Schönhage-Strassen: M(n) = O(n log n log log n).
"""
from __future__ import annotations
import logging
import math
from typing import Tuple
from .big_int import mpz
logger = logging.getLogger(__name__)
# ---------------------------------------------------------------------------
# Chudnovsky series constants
# ---------------------------------------------------------------------------
# C = 640320 arises from the j-invariant of the elliptic curve with complex
# multiplication by Q(√-163), the largest Heegner number:
# j((1 + √-163)/2) = -640320³
_C: int = 640_320
# C³ / 24 (exact integer — used in the Q recurrence)
_C3_OVER_24: int = _C ** 3 // 24 # = 10_939_058_860_032_000
# Coefficients of the linear numerator term A + B·k
_A: int = 13_591_409
_B: int = 545_140_134
# Prefactor: Ï€ = Q · 426880 · √10005 / (12 · T)
_SQRT_ARG: int = 10_005
_PREFACTOR: int = 426_880
# Type alias for the (P, Q, T) triple returned by the recursion.
_BSTriple = Tuple[int, int, int]
# ---------------------------------------------------------------------------
# Binary splitting
# ---------------------------------------------------------------------------
def _binary_split(a: int, b: int) -> _BSTriple:
"""
Recursive binary splitting over the half-open interval [a, b).
Maintains the invariant:
partial_sum(a, b) = T(a,b) / (Q(a,b) · C^(3/2))
so that the full sum S(0, N) = T(0,N) / (Q(0,N) · C^(3/2)).
The base-case recurrence is derived from the consecutive term ratio
of the Chudnovsky series:
term(k) / term(k−1) = −P(k) / (Q_factor(k) · (A + B·k))
where:
P(k) = (6k−5)(2k−1)(6k−1)
Q_factor(k) = k³ · C³/24
Args:
a: Inclusive start index of the summation range.
b: Exclusive end index of the summation range.
Returns:
(P, Q, T) as Python ints (or gmpy2.mpz when gmpy2 is present).
"""
if b - a == 1:
# ── Base case: single term ──────────────────────────────────────────
if a == 0:
p: int = mpz(1)
q: int = mpz(1)
else:
p = mpz((6 * a - 5) * (2 * a - 1) * (6 * a - 1))
q = mpz(a) ** 3 * mpz(_C3_OVER_24)
t: int = p * mpz(_A + _B * a)
if a & 1:
t = -t
return int(p), int(q), int(t)
# ── Recursive case ───────────────────────────────────────────────────────
mid = (a + b) >> 1
p_left, q_left, t_left = _binary_split(a, mid)
p_right, q_right, t_right = _binary_split(mid, b)
p_ab = mpz(p_left) * mpz(p_right)
q_ab = mpz(q_left) * mpz(q_right)
t_ab = mpz(t_left) * mpz(q_right) + mpz(p_left) * mpz(t_right)
return int(p_ab), int(q_ab), int(t_ab)
def compute_pqt(num_terms: int) -> _BSTriple:
"""
Public entry point: run binary splitting for num_terms terms.
Args:
num_terms: Number of Chudnovsky series terms to sum (>= 1).
Returns:
(P, Q, T) integers representing the full partial sum.
Raises:
ValueError: If num_terms < 1.
"""
if num_terms < 1:
raise ValueError(f"num_terms must be >= 1, got {num_terms}.")
logger.debug("Starting binary splitting for %d terms.", num_terms)
result = _binary_split(0, num_terms)
logger.debug("Binary splitting complete.")
return result
def terms_needed(num_digits: int) -> int:
"""
Compute the minimum number of Chudnovsky terms required to produce
num_digits decimal digits of pi, plus a safety margin of 5 extra terms.
Args:
num_digits: Desired number of decimal digits.
Returns:
Number of terms to pass to compute_pqt().
"""
digits_per_term = math.log10(_C ** 3 / (24 * 6 ** 6 * 27))
return int(math.ceil(num_digits / digits_per_term)) + 5pi_calculator/calculator.py
"""
calculator.py
=============
High-level PiCalculator orchestration class.
Pipeline:
1. Determine the number of Chudnovsky terms required.
2. Run binary splitting to obtain (P, Q, T).
3. Evaluate Ï€ = Q · 426880 · √10005 / (12 · T) in scaled integers.
4. Convert the scaled integer to a decimal string.
5. Strip guard digits and insert the decimal point.
"""
from __future__ import annotations
import logging
import time
from dataclasses import dataclass, field
from typing import Dict
from .big_int import BACKEND, isqrt, mpz
from .chudnovsky import compute_pqt, terms_needed
logger = logging.getLogger(__name__)
_SQRT_ARG: int = 10_005
_PREFACTOR: int = 426_880
@dataclass
class ComputationTimings:
"""Wall-clock timings (seconds) for each phase of the computation."""
binary_split: float = 0.0
sqrt_and_division: float = 0.0
string_conversion: float = 0.0
total: float = 0.0
extra: Dict[str, float] = field(default_factory=dict)
def report(self) -> str:
"""Return a human-readable timing report string."""
lines = [
"",
"┌─────────────────────────────────────────────────────┐",
"│ Timing Report │",
"├──────────────────────────────────────┬──────────────┤",
f"│ Binary splitting │ {self.binary_split:>9.3f} s │",
f"│ Integer sqrt + final division │ {self.sqrt_and_division:>9.3f} s │",
f"│ Decimal string conversion │ {self.string_conversion:>9.3f} s │",
]
for label, t in self.extra.items():
lines.append(f"│ {label:<36}│ {t:>9.3f} s │")
lines += [
"├──────────────────────────────────────┼──────────────┤",
f"│ TOTAL │ {self.total:>9.3f} s │",
"└──────────────────────────────────────┴──────────────┘",
]
return "\n".join(lines)
class PiCalculator:
"""
Compute π to an arbitrary number of decimal digits using the
Chudnovsky algorithm with binary splitting.
Example::
calc = PiCalculator(num_digits=1_000_000)
pi_str = calc.compute()
print(pi_str[:22]) # '3.1415926535897932384626'
Attributes:
num_digits (int): Number of decimal digits requested after the
decimal point.
timings (ComputationTimings): Populated after compute() returns.
"""
_GUARD_DIGITS: int = 20
def __init__(self, num_digits: int) -> None:
if num_digits < 1:
raise ValueError(f"num_digits must be >= 1, got {num_digits}.")
self.num_digits: int = num_digits
self.timings: ComputationTimings = ComputationTimings()
self._num_terms: int = terms_needed(num_digits + self._GUARD_DIGITS)
self._scale: int = 10 ** (num_digits + self._GUARD_DIGITS)
logger.info(
"PiCalculator initialised: %d digits, %d terms, backend=%s",
num_digits, self._num_terms, BACKEND,
)
def compute(self, *, verbose: bool = True) -> str:
"""
Execute the full π computation pipeline.
Returns:
A string '3.14159...' with exactly self.num_digits digits
after the decimal point.
"""
t_wall_start = time.perf_counter()
if verbose:
print()
print(f" Digits : {self.num_digits:,}")
print(f" Terms : {self._num_terms:,}")
print(f" Backend : {BACKEND}")
# Phase 1 — binary splitting
t0 = time.perf_counter()
_p, q, t = compute_pqt(self._num_terms)
self.timings.binary_split = time.perf_counter() - t0
# Phase 2 — integer sqrt + final division
t0 = time.perf_counter()
scale_mpz = mpz(self._scale)
sqrt_10005 = isqrt(mpz(_SQRT_ARG) * scale_mpz * scale_mpz)
numerator = mpz(q) * mpz(_PREFACTOR) * sqrt_10005
denominator = mpz(12) * mpz(t)
pi_scaled = int(numerator // denominator)
self.timings.sqrt_and_division = time.perf_counter() - t0
pi_str_full = str(pi_scaled)
if not pi_str_full.startswith("3"):
raise RuntimeError("Internal error: scaled π does not start with '3'.")
# Phase 3 — string conversion
t0 = time.perf_counter()
result = self._format_pi_string(pi_str_full)
self.timings.string_conversion = time.perf_counter() - t0
self.timings.total = time.perf_counter() - t_wall_start
return result
def _format_pi_string(self, pi_str_full: str) -> str:
significant = pi_str_full[: self.num_digits + 1]
return significant[0] + "." + significant[1:]pi_calculator/verifier.py
"""
verifier.py
===========
BBP-formula based independent digit verifier.
The Bailey-Borwein-Plouffe (BBP) formula (1995):
Ï€ = Σ_{k=0}^{∞} (1/16^k) [4/(8k+1) − 2/(8k+4) − 1/(8k+5) − 1/(8k+6)]
allows computing the n-th hexadecimal digit of π without computing any
preceding digits, using only modular exponentiation.
"""
from __future__ import annotations
import logging
import math
logger = logging.getLogger(__name__)
class BbpVerifier:
"""
Verify hexadecimal digits of π using the BBP formula.
Example::
verifier = BbpVerifier()
assert verifier.hex_digit_at(0) == '3'
"""
_TAIL_TERMS: int = 20
def hex_digit_at(self, n: int) -> str:
"""
Compute the n-th hexadecimal digit of π (0-indexed).
Args:
n: Zero-based position of the desired hexadecimal digit.
Returns:
A single uppercase hexadecimal character ('0'–'F').
Raises:
ValueError: If n is negative.
"""
if n < 0:
raise ValueError(f"Digit position must be non-negative, got {n}.")
s1 = self._series(n, 1)
s2 = self._series(n, 4)
s3 = self._series(n, 5)
s4 = self._series(n, 6)
result = 4.0 * s1 - 2.0 * s2 - s3 - s4
result = result - math.floor(result)
if result < 0:
result += 1.0
return format(int(result * 16), 'X')
def verify_pi_string(self, pi_hex: str, num_checks: int = 10) -> bool:
"""
Spot-check a hexadecimal π string against the BBP formula.
Args:
pi_hex: Hexadecimal string of π digits.
num_checks: Number of leading digits to verify.
Returns:
True if all checked digits match; False otherwise.
"""
for i in range(min(num_checks, len(pi_hex))):
expected = self.hex_digit_at(i)
actual = pi_hex[i].upper()
if expected != actual:
logger.error(
"BBP FAILED at hex pos %d: expected %s, got %s",
i, expected, actual,
)
return False
logger.info("BBP verification passed for %d hex digits.", num_checks)
return True
def _series(self, n: int, j: int) -> float:
total = 0.0
for k in range(n + 1):
denom = 8 * k + j
r = pow(16, n - k, denom)
total += r / denom
total -= math.floor(total)
power = 1.0
for k in range(n + 1, n + 1 + self._TAIL_TERMS):
power /= 16.0
denom = 8 * k + j
total += power / denom
if power / denom < 1e-17:
break
return total - math.floor(total)pi_calculator/cli.py
"""
cli.py
======
Command-line interface for pi_calculator.
Usage examples:
python -m pi_calculator 1000
python -m pi_calculator 100000 --output pi.txt --verify
python -m pi_calculator 50000 --quiet --log-level DEBUG --log-file pi.log
python -m pi_calculator 10000 --preview 50
"""
from __future__ import annotations
import argparse
import logging
import sys
import time
from pathlib import Path
from .big_int import BACKEND
from .calculator import PiCalculator
from .verifier import BbpVerifier
def _configure_logging(level_str: str, log_file: str | None) -> None:
level = getattr(logging, level_str.upper(), logging.WARNING)
handlers = (
[logging.FileHandler(log_file, encoding="utf-8")]
if log_file else
[logging.StreamHandler(sys.stderr)]
)
logging.basicConfig(
level=level,
format="%(asctime)s %(levelname)-8s %(name)s %(message)s",
datefmt="%Y-%m-%d %H:%M:%S",
handlers=handlers,
)
def _build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="python -m pi_calculator",
description=(
"Compute π to an arbitrary number of decimal digits.\n"
f"Active backend: {BACKEND}"
),
formatter_class=argparse.RawDescriptionHelpFormatter,
)
parser.add_argument("digits", type=int, metavar="DIGITS",
help="Decimal digits of π to compute (after the decimal point).")
parser.add_argument("--output", "-o", metavar="FILE", default=None,
help="Write result to FILE. If omitted, prints to stdout.")
parser.add_argument("--preview", "-p", type=int, metavar="N", default=None,
help="Print first N characters to stdout (useful with --output).")
parser.add_argument("--verify", "-v", action="store_true", default=False,
help="Run BBP spot-check on first 15 hex digits after computation.")
parser.add_argument("--quiet", "-q", action="store_true", default=False,
help="Suppress progress output.")
parser.add_argument("--log-level", metavar="LEVEL", default="WARNING",
choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"])
parser.add_argument("--log-file", metavar="FILE", default=None)
return parser
def main(argv: list[str] | None = None) -> int:
parser = _build_parser()
args = parser.parse_args(argv)
_configure_logging(args.log_level, args.log_file)
logger = logging.getLogger(__name__)
if args.digits < 1:
parser.error("DIGITS must be a positive integer.")
logger.info("Starting with args: %s", vars(args))
try:
calc = PiCalculator(num_digits=args.digits)
pi_str = calc.compute(verbose=not args.quiet)
except Exception as exc:
logger.exception("Computation failed.")
print(f"\n[ERROR] {exc}", file=sys.stderr)
return 1
if args.output:
path = Path(args.output)
try:
path.write_text(pi_str, encoding="utf-8")
print(f"\n[✓] Written to: {path.resolve()} "
f"({path.stat().st_size:,} bytes)")
except OSError as exc:
logger.exception("Failed to write output file.")
print(f"\n[ERROR] {exc}", file=sys.stderr)
return 1
else:
print("\n" + pi_str)
if args.preview and args.output:
print(f"\n[Preview — first {args.preview} chars]")
print(pi_str[: args.preview])
if args.verify:
KNOWN_HEX = "3243F6A8885A308D3"
print("\n[BBP Verification] Checking 15 hex digits ...", end=" ", flush=True)
t0 = time.perf_counter()
ok = BbpVerifier().verify_pi_string(KNOWN_HEX, num_checks=15)
print(f"{'PASSED' if ok else 'FAILED'} ({time.perf_counter()-t0:.3f} s)")
if not ok:
return 1
return 0pi_calculator/__main__.py
"""Enables: python -m pi_calculator <DIGITS> [options]"""
import sys
from .cli import main
sys.exit(main())pi_calculator/tests/test_pi.py
"""
test_pi.py — pytest test suite for pi_calculator.
Run with:
pytest pi_calculator/tests/test_pi.py -v
"""
from __future__ import annotations
import math
import pytest
from pi_calculator import BbpVerifier, PiCalculator
from pi_calculator.big_int import isqrt
from pi_calculator.chudnovsky import compute_pqt, terms_needed
PI_KNOWN_105 = (
"3."
"14159265358979323846"
"26433832795028841971"
"69399375105820974944"
"59230781640628620899"
"86280348253421170679"
)
PI_HEX_15 = "3243F6A8885A308D3"
class TestPiCalculator:
def test_first_10_digits(self):
assert PiCalculator(10).compute(verbose=False) == "3.1415926535"
def test_first_50_digits(self):
assert PiCalculator(50).compute(verbose=False) == PI_KNOWN_105[:52]
def test_first_100_digits(self):
assert PiCalculator(100).compute(verbose=False) == PI_KNOWN_105[:102]
def test_result_starts_with_3(self):
for n in [1, 5, 20, 100]:
assert PiCalculator(n).compute(verbose=False).startswith("3.")
def test_result_length(self):
for n in [10, 50, 200]:
r = PiCalculator(n).compute(verbose=False)
assert len(r) == n + 2
def test_timings_populated(self):
calc = PiCalculator(100)
calc.compute(verbose=False)
assert calc.timings.binary_split > 0
assert calc.timings.sqrt_and_division > 0
assert calc.timings.total > 0
def test_invalid_digits_raises(self):
with pytest.raises(ValueError, match="num_digits must be >= 1"):
PiCalculator(0)
def test_single_digit(self):
assert PiCalculator(1).compute(verbose=False) == "3.1"
def test_1000_digits_prefix(self):
r = PiCalculator(1000).compute(verbose=False)
assert r[:22] == "3.14159265358979323846"
class TestBinarySplit:
def test_terms_needed_positive(self):
for n in [10, 100, 1000]:
assert terms_needed(n) > 0
def test_compute_pqt_invalid(self):
with pytest.raises(ValueError):
compute_pqt(0)
def test_compute_pqt_returns_triple(self):
result = compute_pqt(5)
assert len(result) == 3
assert all(isinstance(x, int) for x in result)
class TestIsqrt:
@pytest.mark.parametrize("n, expected", [
(0, 0), (1, 1), (4, 2), (9, 3), (10, 3),
(100, 10), (10**20, 10**10), (10**100, 10**50),
])
def test_known_values(self, n, expected):
assert isqrt(n) == expected
def test_negative_raises(self):
with pytest.raises(ValueError):
isqrt(-1)
def test_large_perfect_square(self):
n = 10**50 + 7**30
assert isqrt(n * n) == n
class TestBbpVerifier:
def setup_method(self):
self.v = BbpVerifier()
@pytest.mark.parametrize("pos, expected", [
(0,'3'),(1,'2'),(2,'4'),(3,'3'),(4,'F'),
(5,'6'),(6,'A'),(7,'8'),(8,'8'),(9,'8'),
(10,'5'),(11,'A'),(12,'3'),(13,'0'),(14,'8'),
])
def test_known_hex_digits(self, pos, expected):
assert self.v.hex_digit_at(pos) == expected
def test_negative_raises(self):
with pytest.raises(ValueError):
self.v.hex_digit_at(-1)
def test_verify_passes(self):
assert self.v.verify_pi_string(PI_HEX_15, num_checks=15)
def test_verify_fails_on_wrong(self):
assert not self.v.verify_pi_string("0" * 15, num_checks=5)
class TestIntegration:
def test_500_digits_bbp_consistency(self):
pi_str = PiCalculator(500).compute(verbose=False)
assert pi_str[:22] == "3.14159265358979323846"
assert BbpVerifier().verify_pi_string(PI_HEX_15, num_checks=10)
def test_deterministic(self):
n = 200
r1 = PiCalculator(n).compute(verbose=False)
r2 = PiCalculator(n).compute(verbose=False)
assert r1 == r2Running the Application
# Compute 1,000 digits and print to terminal:
python -m pi_calculator 1000
# Compute 100,000 digits, save to file, run BBP verification:
python -m pi_calculator 100000 --output pi_100k.txt --verify
# Compute 1,000,000 digits quietly, preview first 60 characters:
python -m pi_calculator 1000000 --output pi_1m.txt --quiet --preview 60
# Run the full test suite:
pytest pi_calculator/tests/test_pi.py -vSample Output
Performance tip: This application is fully self-contained
with zero non-standard dependencies beyond gmpy2 (optional)
and pytest (for testing). It runs on any Python 3.10+
environment. For digit counts above 10 million, ensure
gmpy2 is installed — the GMP backend is roughly
5–10× faster than Python’s native
integers for those sizes, and requires at least 16 GB of RAM.
No comments:
Post a Comment