Chapter 25 — Irrational Numbers

# --- Visualization: Rational approximations converging to √2 and π ---
import numpy as np
import matplotlib.pyplot as plt
from fractions import Fraction
plt.style.use('seaborn-v0_8-whitegrid')

def continued_fraction_convergents(coeffs):
    """
    Compute convergents of a continued fraction [a0; a1, a2, ...].

    Args:
        coeffs: list of int continued fraction coefficients

    Returns:
        list of Fraction: convergents p_k/q_k
    """
    convergents = []
    p_prev, p_curr = 1, coeffs[0]
    q_prev, q_curr = 0, 1
    convergents.append(Fraction(p_curr, q_curr))

    for a in coeffs[1:]:
        p_prev, p_curr = p_curr, a * p_curr + p_prev
        q_prev, q_curr = q_curr, a * q_curr + q_prev
        convergents.append(Fraction(p_curr, q_curr))

    return convergents

# √2 = [1; 2, 2, 2, 2, ...]
sqrt2_coeffs = [1] + [2] * 15
sqrt2_convergents = continued_fraction_convergents(sqrt2_coeffs)

# π ≈ [3; 7, 15, 1, 292, 1, 1, 1, 2, ...]
pi_coeffs = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14]
pi_convergents = continued_fraction_convergents(pi_coeffs)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

for ax, convergents, true_val, name, color in [
    (axes[0], sqrt2_convergents, np.sqrt(2), '√2', 'steelblue'),
    (axes[1], pi_convergents,   np.pi,      'π',  'crimson')
]:
    errors = [abs(float(c) - true_val) for c in convergents]
    denoms = [c.denominator for c in convergents]

    ax.semilogy(range(len(errors)), errors, 'o-', color=color, linewidth=2)
    for i, (c, e) in enumerate(zip(convergents, errors)):
        if i < 8:
            ax.annotate(f'{c}', (i, e), textcoords='offset points',
                        xytext=(5, 5), fontsize=7, color='gray')

    ax.axhline(1e-15, color='gray', linestyle=':', linewidth=1,
               label='float64 precision limit')
    ax.set_xlabel('Convergent index k')
    ax.set_ylabel(f'|convergent - {name}|')
    ax.set_title(f'Continued Fraction Convergents to {name}')
    ax.legend(fontsize=8)

plt.tight_layout()
plt.show()

print("Best approximations to π:")
for c in pi_convergents[:6]:
    print(f"  {c} = {float(c):.10f}  (error: {abs(float(c)-np.pi):.2e})")

# --- Implementation: Continued fraction expansion and irrationality test ---
import math
from fractions import Fraction

def to_continued_fraction(x, max_terms=20, tol=1e-10):
    """
    Compute the continued fraction expansion of float x.

    Args:
        x: float
        max_terms: int, maximum number of coefficients
        tol: float, stop when fractional part is smaller than this

    Returns:
        list of int: continued fraction coefficients [a0, a1, a2, ...]
    """
    coeffs = []
    for _ in range(max_terms):
        a = int(x)         # floor
        coeffs.append(a)
        frac = x - a       # fractional part
        if frac < tol:
            break
        x = 1.0 / frac     # reciprocal of fractional part
    return coeffs


def is_likely_rational(x, max_denom=10000):
    """
    Heuristic: check if x is 'close enough' to a rational with small denominator.

    Args:
        x: float
        max_denom: int

    Returns:
        Fraction or None
    """
    f = Fraction(x).limit_denominator(max_denom)
    if abs(float(f) - x) < 1e-10:
        return f
    return None


print("Continued fraction expansions:")
test_vals = [
    (math.sqrt(2), '√2'),
    (math.pi, 'π'),
    (math.e, 'e'),
    ((1 + math.sqrt(5)) / 2, 'φ (golden ratio)'),
    (22/7, '22/7'),
    (355/113, '355/113'),
]
for val, name in test_vals:
    cf = to_continued_fraction(val, max_terms=10)
    rat = is_likely_rational(val)
    print(f"  {name:20s}: [{'; '.join(map(str, cf[:6]))}...]  rational? {rat}")

# --- Experiment 1: The golden ratio has the 'worst' rational approximations ---
# Hypothesis: φ = [1;1,1,1,...] — all 1s — converges slowest to irrational
# Try changing: TARGET to other irrationals and compare convergence rate

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

phi = (1 + np.sqrt(5)) / 2  # golden ratio
TARGET = phi  # <-- try np.sqrt(2), np.pi, np.e

cf_coeffs = to_continued_fraction(TARGET, max_terms=20)
convs = continued_fraction_convergents(cf_coeffs)
errors = [abs(float(c) - TARGET) for c in convs]

fig, ax = plt.subplots(figsize=(9, 4))
ax.semilogy(errors, 'o-', color='purple')
ax.set_xlabel('Convergent index')
ax.set_ylabel('Approximation error (log scale)')
ax.set_title(f'Convergence of continued fraction to {TARGET:.6f}')
plt.tight_layout()
plt.show()

print(f"CF coefficients: {cf_coeffs}")
print("(All 1s = golden ratio — slowest possible convergence)")

# --- Experiment 2: Verify √n is irrational iff n is not a perfect square ---
# Hypothesis: for non-perfect-square n, continued fraction of √n is periodic
# Try changing: N_MAX

N_MAX = 20  # <-- modify this

import math
print(f"{'n':>4}  {'Is perfect square?':>20}  {'CF coefficients':>40}")
print('-' * 70)
for n in range(2, N_MAX + 1):
    root = math.sqrt(n)
    is_perf = int(root) ** 2 == n
    cf = to_continued_fraction(root, max_terms=12)
    print(f"{n:>4}  {'Yes' if is_perf else 'No':>20}  {str(cf[:10]):>40}")

# --- Mini Project: π Approximation Benchmark ---
# Problem: Multiple historical formulas for π converge at different rates.
#          Which is best per computation unit?
# Task: Implement and benchmark 3 algorithms, plot their convergence.

import math
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

def leibniz_pi(n_terms):
    """Leibniz formula: π/4 = 1 - 1/3 + 1/5 - 1/7 + ..."""
    # TODO: implement. Returns approximation of π after n_terms.
    pass

def nilakantha_pi(n_terms):
    """Nilakantha series: π = 3 + 4/(2·3·4) - 4/(4·5·6) + ..."""
    # TODO: implement.
    pass

def wallis_pi(n_terms):
    """Wallis product: π/2 = (2/1)(2/3)(4/3)(4/5)(6/5)(6/7)..."""
    # TODO: implement.
    pass

N_TERMS = 1000
ns = range(1, N_TERMS + 1)

methods = [
    ('Leibniz', leibniz_pi, 'steelblue'),
    ('Nilakantha', nilakantha_pi, 'crimson'),
    ('Wallis', wallis_pi, 'green'),
]

fig, ax = plt.subplots(figsize=(10, 5))
for name, fn, color in methods:
    if fn(1) is not None:  # only plot implemented methods
        errors = [abs(fn(n) - math.pi) for n in ns]
        ax.semilogy(list(ns), errors, label=name, color=color, linewidth=1.5)

ax.set_xlabel('Number of terms')
ax.set_ylabel('|approximation - π|')
ax.set_title('π Approximation: Convergence Comparison')
ax.legend()
plt.tight_layout()
plt.show()