Chapter 37 — Large Numbers in Computing

# --- Visualization: Integer overflow in fixed-width types ---
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')

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

# Left: int8 overflow (range -128 to 127)
x = np.arange(-200, 200, 1, dtype=np.int32)  # true values
x_int8 = x.astype(np.int8)                    # will overflow outside [-128, 127]

axes[0].plot(x, x, 'gray', linewidth=1, linestyle='--', label='True value (int32)')
axes[0].plot(x, x_int8, 'crimson', linewidth=1.5, label='int8 (overflows!)')
axes[0].set_xlabel('Input', fontsize=11)
axes[0].set_ylabel('Stored value', fontsize=11)
axes[0].set_title('int8 Overflow: wraps at ±128', fontsize=12)
axes[0].legend()
axes[0].axvline(127, color='blue', linestyle=':', alpha=0.5, label='MAX_INT8')
axes[0].axvline(-128, color='green', linestyle=':', alpha=0.5, label='MIN_INT8')

# Right: uint64 factorial overflow
factorials_true = []
factorials_uint64 = []
f = 1
MAX_UINT64 = 2**64 - 1

for n in range(1, 25):
    f *= n
    factorials_true.append(f)
    try:
        factorials_uint64.append(int(np.uint64(f % (MAX_UINT64 + 1))))
    except (OverflowError, ValueError):
        factorials_uint64.append(float('nan'))

ns = list(range(1, 25))
axes[1].semilogy(ns, factorials_true, 'steelblue', linewidth=2, label='True n! (Python int)')
overflow_idx = next((i for i, v in enumerate(factorials_true) if v > MAX_UINT64), None)
if overflow_idx:
    axes[1].axvline(overflow_idx + 1, color='crimson', linestyle='--',
                    label=f'uint64 overflows at n={overflow_idx+1}')
axes[1].axhline(MAX_UINT64, color='orange', linestyle=':', linewidth=1, label='uint64 max')
axes[1].set_xlabel('n', fontsize=11)
axes[1].set_ylabel('n! (log scale)', fontsize=11)
axes[1].set_title('Factorial growth vs uint64 capacity', fontsize=12)
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.show()

# --- Implementation: Safe large number operations ---
import math
import sys

def log_factorial(n):
    """
    Compute log(n!) in float64 without overflow.
    Uses sum of logs instead of computing n! directly.

    Args:
        n: non-negative integer

    Returns:
        log(n!) as float64
    """
    return sum(math.log(k) for k in range(1, n + 1))


def stirling_log_factorial(n):
    """
    Stirling's approximation: log(n!) ≈ 0.5*log(2πn) + n*log(n) - n
    Accurate to ~0.1% for n > 10, improves as n grows.
    """
    if n == 0:
        return 0.0
    return 0.5 * math.log(2 * math.pi * n) + n * math.log(n) - n


def digit_count(n):
    """Number of decimal digits of positive integer n."""
    return math.floor(math.log10(n)) + 1 if n > 0 else 1


def python_int_size_bytes(n):
    """Memory used by Python integer n."""
    return sys.getsizeof(n)


def safe_binomial(n, k):
    """
    Compute C(n, k) = n! / (k! (n-k)!) via log-space to avoid intermediate overflow,
    then return as Python int (exact).
    """
    if k < 0 or k > n:
        return 0
    log_result = log_factorial(n) - log_factorial(k) - log_factorial(n - k)
    return round(math.exp(log_result))


# Demonstrations
print("=== Factorial digit counts ===")
for n in [10, 50, 100, 1000, 10000]:
    log_f = log_factorial(n)
    stir  = stirling_log_factorial(n)
    digits = int(log_f / math.log(10)) + 1
    stir_err = abs(stir - log_f) / abs(log_f) * 100
    print(f"  {n}!: {digits} digits, Stirling error: {stir_err:.4f}%")

print("\n=== Memory cost of big integers ===")
f = 1
for n in [1, 10, 100, 1000, 10000]:
    while True:
        f *= n if f == 1 else 1  # we'll just compute directly
        break

import math as _math
for n in [1, 100, 1000, 10_000]:
    big = _math.factorial(n)
    print(f"  {n}!: {sys.getsizeof(big)} bytes, {big.bit_length()} bits")

print("\n=== Binomial coefficients (large) ===")
for n, k in [(52, 5), (100, 50), (1000, 100)]:
    result = safe_binomial(n, k)
    print(f"  C({n},{k}) ≈ {result:.6e}")

# --- Experiment 1: Cost of big integer arithmetic ---
# Hypothesis: multiplication time grows faster than linearly with digit count.
# Try changing: MAX_N to measure at larger scales.

import time
import math
import matplotlib.pyplot as plt

MAX_N = 5000   # <-- modify: try 1000, 10000
step  = MAX_N // 20
ns    = list(range(step, MAX_N + 1, step))

times = []
for n in ns:
    # Time the factorial computation (builds up through multiplication)
    t0 = time.perf_counter()
    _ = math.factorial(n)
    t1 = time.perf_counter()
    times.append(t1 - t0)

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(ns, times, 'o-', color='steelblue', markersize=4)
ax.set_xlabel('n (computing n!)', fontsize=11)
ax.set_ylabel('Time (seconds)', fontsize=11)
ax.set_title(f'Time to compute n! for n up to {MAX_N}', fontsize=12)
plt.tight_layout()
plt.show()

# Fit a power law: time ~ n^alpha
import numpy as np
log_ns = np.log(ns)
log_ts = np.log(times)
alpha, _ = np.polyfit(log_ns, log_ts, 1)
print(f"Fitted exponent: time ~ n^{alpha:.2f}")
print("Naive multiplication: O(n^2). Python uses more advanced algorithms for large n.")

# --- Experiment 2: Float overflow and inf arithmetic ---
# Hypothesis: float overflow saturates to inf; subsequent arithmetic follows IEEE rules.
# Try changing: BASE and observe when overflow occurs.

import numpy as np

BASE = 2.0   # <-- modify

print("Float exponent ladder (base 2):")
for exp in [100, 500, 1000, 1022, 1023, 1024, 1100]:
    try:
        val = float(BASE) ** exp
    except OverflowError:
        val = float('inf')
    print(f"  {BASE}^{exp} = {val:.6e}  (inf={np.isinf(val)})")

print("\nInf arithmetic rules (IEEE 754):")
inf = float('inf')
nan = float('nan')
cases = [
    (inf + 1, 'inf + 1'),
    (inf - inf, 'inf - inf'),
    (inf * 0, 'inf * 0'),
    (inf / inf, 'inf / inf'),
    (1 / inf, '1 / inf'),
    (inf * -1, 'inf * -1'),
]
for result, expr in cases:
    print(f"  {expr} = {result}")

# --- Mini Project: Big Integer Arithmetic ---
# Represent non-negative integers as lists of decimal digits, least significant first.
# e.g., 12345 → [5, 4, 3, 2, 1]

def int_to_digits(n):
    """Convert non-negative int to list of digits, LSB first."""
    if n == 0: return [0]
    digits = []
    while n:
        digits.append(n % 10)
        n //= 10
    return digits


def digits_to_int(digits):
    """Convert digit list (LSB first) to int."""
    return sum(d * (10 ** i) for i, d in enumerate(digits))


def bigint_add(a_digits, b_digits):
    """
    Add two big integers represented as digit lists (LSB first).
    Uses standard carry-based addition.

    Args:
        a_digits, b_digits: lists of single digits (0-9), LSB first

    Returns:
        result as digit list (LSB first)
    """
    result = []
    carry = 0
    n = max(len(a_digits), len(b_digits))
    for i in range(n):
        a = a_digits[i] if i < len(a_digits) else 0
        b = b_digits[i] if i < len(b_digits) else 0
        total = a + b + carry
        result.append(total % 10)
        carry = total // 10
    if carry:
        result.append(carry)
    return result


def bigint_mul(a_digits, b_digits):
    """
    Multiply two big integers using schoolbook O(n^2) algorithm.

    Args:
        a_digits, b_digits: digit lists (LSB first)

    Returns:
        product as digit list (LSB first)
    """
    n, m = len(a_digits), len(b_digits)
    result = [0] * (n + m)
    for i, a in enumerate(a_digits):
        for j, b in enumerate(b_digits):
            result[i + j] += a * b
    # Normalize carries
    for k in range(len(result) - 1):
        result[k + 1] += result[k] // 10
        result[k] %= 10
    # Remove leading zeros
    while len(result) > 1 and result[-1] == 0:
        result.pop()
    return result


# Validate
test_pairs = [(12345, 67890), (999, 1), (100, 100), (12345678, 87654321)]
print("Validation:")
for a, b in test_pairs:
    da, db = int_to_digits(a), int_to_digits(b)
    sum_result = digits_to_int(bigint_add(da, db))
    mul_result = digits_to_int(bigint_mul(da, db))
    print(f"  {a} + {b} = {sum_result} (correct: {a+b==sum_result})")
    print(f"  {a} * {b} = {mul_result} (correct: {a*b==mul_result})")

# Compute 100! using our big integer library
import math
result = int_to_digits(1)
for k in range(1, 101):
    result = bigint_mul(result, int_to_digits(k))
our_factorial = digits_to_int(result)
true_factorial = math.factorial(100)
print(f"\n100! matches Python: {our_factorial == true_factorial}")
print(f"100! has {len(result)} decimal digits")