Chapter 26 — Real Numbers

# --- Visualization: Float64 as a discrete approximation of ℝ ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

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

# Plot 1: spacing between consecutive float64 values (machine epsilon)
ax = axes[0]
magnitudes = np.logspace(-5, 15, 500)
epsilons = [float(np.finfo(float).eps) * x for x in magnitudes]
ax.loglog(magnitudes, epsilons, color='steelblue', linewidth=2)
ax.set_xlabel('x')
ax.set_ylabel('Spacing to next float64')
ax.set_title('Float64 Spacing: Relative Error is Constant')

# Plot 2: coverage of float64 on [0, 2]
ax = axes[1]
# Show float64 values in [1, 1 + 1e-14]
start = 1.0
step = np.finfo(float).eps  # machine epsilon at 1.0
n_vals = 20
vals = [start + i * step for i in range(n_vals)]
ax.scatter(vals, [0]*len(vals), s=40, color='crimson')
ax.set_xlabel('Value')
ax.set_yticks([])
ax.set_title(f'Float64 values near 1.0: spacing = {step:.2e}')

plt.tight_layout()
plt.show()

# Show key float64 properties
print(f"Machine epsilon at 1.0: {np.finfo(float).eps:.2e}")
print(f"Max float64:            {np.finfo(float).max:.2e}")
print(f"Min positive float64:   {np.finfo(float).tiny:.2e}")
print(f"Approx # float64 vals:  2^64 ≈ {2**64:.2e}")

# --- Implementation: Cauchy sequence convergence checker and interval arithmetic ---
import numpy as np

def cauchy_test(seq_fn, n_terms=100, tol=1e-10):
    """
    Test if a sequence appears to be Cauchy (converging in ℝ).

    Args:
        seq_fn: function int -> float, computes nth term
        n_terms: int
        tol: float, tolerance for convergence

    Returns:
        (bool, float): (is_cauchy, apparent_limit)
    """
    terms = [seq_fn(n) for n in range(1, n_terms + 1)]
    # Check if later terms are close together
    last_diff = abs(terms[-1] - terms[-2])
    is_cauchy = last_diff < tol
    return is_cauchy, terms[-1]

# Test with sequences converging to irrationals
sequences = [
    ('√2 via Newton', lambda n: sum(1/(2**k) * ((-1)**k if k > 0 else 1)
                                    for k in range(n)) if False else
                      lambda n, x=1.0: [x := (x + 2/x)/2 for _ in range(n)][-1]),
    ('Leibniz π/4',   lambda n: 4 * sum((-1)**k / (2*k+1) for k in range(n))),
    ('e via series',  lambda n: sum(1/np.math.factorial(k) for k in range(n)) if False
                      else lambda n: sum(1/np.prod(range(1,k+1)) if k > 0 else 1
                                        for k in range(n+1))),
]

# Just demonstrate Newton's method for √2
print("Newton's method for √2:")
x = 1.0
for i in range(10):
    x = (x + 2/x) / 2
    print(f"  step {i+1}: {x:.15f}  error: {abs(x - np.sqrt(2)):.2e}")

# --- Experiment: Completeness failure in ℚ ---
# Hypothesis: a Cauchy sequence in ℚ can converge to an irrational
# Try changing: N_TERMS

import numpy as np
N_TERMS = 30  # <-- modify this

# Newton's method for √2, starting from 1
# Each step: x_{n+1} = (x_n + 2/x_n) / 2
# All x_n are rational, but limit is irrational

from fractions import Fraction
x = Fraction(1)
true_val = float(np.sqrt(2))
print(f"{'n':>4}  {'x_n':>30}  {'float(x_n)':>18}  {'error':>12}")
for n in range(1, min(N_TERMS, 10) + 1):
    x = (x + Fraction(2, 1) / x) / 2
    err = abs(float(x) - true_val)
    print(f"{n:>4}  {str(x)[:30]:>30}  {float(x):>18.15f}  {err:>12.3e}")
print("All x_n are rational. Limit is √2 — irrational. ℚ is not complete.")

# --- Mini Project: Bisection Root Finder ---
# Problem: Find real roots of polynomials that have no closed-form expression.
# Task: Implement bisection, Newton's method, and compare convergence rates.

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

def bisection(f, a, b, tol=1e-12, max_iter=100):
    """
    Find root of f in [a, b] using bisection.
    Requires f(a)*f(b) < 0.
    Returns (root, list_of_midpoints).
    """
    # TODO: implement
    pass

def newton(f, df, x0, tol=1e-12, max_iter=50):
    """
    Find root of f using Newton's method starting from x0.
    df is the derivative of f.
    Returns (root, list_of_iterates).
    """
    # TODO: implement
    pass

# Test: find √2 as root of f(x) = x^2 - 2
f  = lambda x: x**2 - 2
df = lambda x: 2*x

bis_root, bis_hist = bisection(f, 1.0, 2.0)
newt_root, newt_hist = newton(f, df, 1.0)

if bis_hist and newt_hist:
    true_root = np.sqrt(2)
    fig, ax = plt.subplots(figsize=(9, 4))
    if bis_hist:
        ax.semilogy([abs(x - true_root) for x in bis_hist],
                    'o-', label='Bisection (linear)', color='steelblue')
    if newt_hist:
        ax.semilogy([abs(x - true_root) for x in newt_hist],
                    's-', label="Newton's (quadratic)", color='crimson')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Error |x - √2|')
    ax.set_title('Bisection vs Newton: Convergence to √2')
    ax.legend()
    plt.tight_layout()
    plt.show()