Chapter 43 — Logarithms Intuition

# --- Visualization: Logarithm as inverse of exponential ---
import numpy as np
import matplotlib.pyplot as plt

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

x_exp = np.linspace(-3, 3, 300)
x_log = np.linspace(0.01, 20, 300)

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

# Panel 1: e^x and ln(x) are reflections across y=x
ax = axes[0]
ax.plot(x_exp, np.exp(x_exp), color='red',   label='y = e^x (exponential)', linewidth=2)
ax.plot(x_log, np.log(x_log), color='blue',  label='y = ln(x) (logarithm)', linewidth=2)
ax.plot(x_log, x_log,         color='gray',  label='y = x (mirror)',         linestyle='--')
ax.set_xlim(-3, 8)
ax.set_ylim(-3, 8)
ax.set_aspect('equal')
ax.set_title('Exponential and Logarithm\nare Reflections Across y=x')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(fontsize=9)

# Panel 2: Three log bases on the same plot
ax = axes[1]
x_pos = np.linspace(0.01, 100, 500)
ax.plot(x_pos, np.log2(x_pos),  color='steelblue', label='log₂(x) — bits')
ax.plot(x_pos, np.log10(x_pos), color='green',     label='log₁₀(x) — decades')
ax.plot(x_pos, np.log(x_pos),   color='red',       label='ln(x) — nats')
ax.axhline(0, color='black', linewidth=0.5)
ax.set_title('Three Standard Logarithm Bases')
ax.set_xlabel('x')
ax.set_ylabel('logᵦ(x)')
ax.legend()

plt.tight_layout()
plt.show()

# --- Visualization: Logarithmic compression — the zoom-out effect ---
import numpy as np
import matplotlib.pyplot as plt

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

# Real-world quantities spanning many orders of magnitude
LABELS = [
    'Bacterium (1μm)', 'Red blood cell (8μm)', 'Human hair (70μm)',
    'Grain of sand (1mm)', 'Pencil width (7mm)', 'Human (1.7m)',
    'Football field (100m)', 'Mt Everest (8800m)', 'Earth radius (6.4Mm)',
    'Earth–Moon (384Mm)'
]
VALUES_M = [1e-6, 8e-6, 70e-6, 1e-3, 7e-3, 1.7, 100, 8800, 6.4e6, 3.84e8]

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

ax = axes[0]
ax.barh(LABELS, VALUES_M, color='steelblue', alpha=0.7)
ax.set_title('Sizes on Linear Scale (nearly useless — everything looks like zero)')
ax.set_xlabel('Size (meters)')

ax = axes[1]
log_values = np.log10(VALUES_M)
ax.barh(LABELS, log_values, color='orange', alpha=0.7)
ax.set_title('Sizes on log₁₀ Scale (all differences visible)')
ax.set_xlabel('log₁₀(size in meters)')

# Add tick labels in real units
ticks = range(-6, 9, 2)
ax.set_xticks(ticks)
ax.set_xticklabels([f'10^{t}m' for t in ticks], fontsize=8)

plt.tight_layout()
plt.show()

# --- Verify laws of logarithms numerically ---
import numpy as np

x, y, p = 12.0, 7.0, 3.0
b = 2.0  # base

def logb(val, base):
    return np.log(val) / np.log(base)

print("Laws of Logarithms — Numerical Verification  (b=2, x=12, y=7, p=3)")
print("-" * 65)

# Law 1
lhs = logb(x * y, b)
rhs = logb(x, b) + logb(y, b)
print(f"1. log(xy)  = log(x)+log(y):  {lhs:.6f} == {rhs:.6f}  {'✓' if np.isclose(lhs,rhs) else '✗'}")

# Law 2
lhs = logb(x / y, b)
rhs = logb(x, b) - logb(y, b)
print(f"2. log(x/y) = log(x)-log(y): {lhs:.6f} == {rhs:.6f}  {'✓' if np.isclose(lhs,rhs) else '✗'}")

# Law 3
lhs = logb(x**p, b)
rhs = p * logb(x, b)
print(f"3. log(x^p) = p·log(x):      {lhs:.6f} == {rhs:.6f}  {'✓' if np.isclose(lhs,rhs) else '✗'}")

# Change of base
log2_direct = np.log2(x)
log2_change = np.log(x) / np.log(2)
print(f"\nChange of base: log₂({x}) direct={log2_direct:.6f}  via ln={log2_change:.6f}  {'✓' if np.isclose(log2_direct,log2_change) else '✗'}")

# --- Implementation: Computing logarithms from scratch via bisection ---
#
# We cannot use np.log. We find log_b(x) by solving b^y = x for y.
# Since b^y is strictly increasing in y (for b>1), we can binary-search for y.

def log_bisection(x, base, tol=1e-10, max_iter=200):
    """
    Compute log_base(x) by binary search on y such that base^y = x.

    This works because f(y) = base^y is strictly monotone, so we
    can bisect for the root of f(y) - x = 0.

    Args:
        x:        positive real number (argument)
        base:     base > 1
        tol:      tolerance for convergence
        max_iter: maximum bisection steps

    Returns:
        Approximate log_base(x)
    """
    if x <= 0:
        raise ValueError("Logarithm undefined for x <= 0")
    if base <= 0 or base == 1:
        raise ValueError("Base must be positive and not 1")

    # Find initial bracket [lo, hi] such that base^lo < x < base^hi
    lo, hi = -1000.0, 1000.0
    # Shrink the bracket (optional but faster)
    while True:
        try:
            if base ** lo <= x:
                break
        except OverflowError:
            # underflow means base**lo is effectively 0, which is < x
            break
        lo *= 2

    while True:
        try:
            if base ** hi >= x:
                break
        except OverflowError:
            # overflow means base**hi is effectively +inf, which is > x
            break
        hi *= 2

    # Binary search
    for _ in range(max_iter):
        mid = (lo + hi) / 2
        val = base ** mid
        if abs(val - x) < tol * x:  # relative tolerance
            return mid
        if val < x:
            lo = mid
        else:
            hi = mid

    return (lo + hi) / 2


# Validate against numpy
import numpy as np

test_cases = [(8, 2), (1000, 10), (1, 2), (0.125, 2), (np.e**3, np.e)]
print(f"{'x':>10} {'base':>6} {'bisection':>14} {'numpy':>14} {'error':>12}")
print("-" * 60)
for x, b in test_cases:
    our   = log_bisection(x, b)
    ref   = np.log(x) / np.log(b)
    err   = abs(our - ref)
    print(f"{x:>10.4f} {b:>6.4f} {our:>14.8f} {ref:>14.8f} {err:>12.2e}")

# --- Experiment 1: O(log n) vs O(n) — binary search step counts ---
# Hypothesis: Binary search takes log₂(n) steps while linear search takes n steps.
# The difference becomes staggering as n grows.
# Try changing: MAX_N

import numpy as np
import matplotlib.pyplot as plt

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

MAX_N = 1_000_000   # <-- modify this

n_values = np.logspace(0, np.log10(MAX_N), 400)
linear_steps = n_values
log_steps    = np.log2(n_values)

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

ax = axes[0]
ax.plot(n_values, linear_steps, color='red',       label='Linear search: O(n)')
ax.plot(n_values, log_steps,    color='steelblue', label='Binary search: O(log₂n)', linewidth=2)
ax.set_title('Search Steps: Linear vs Binary')
ax.set_xlabel('Array size n')
ax.set_ylabel('Steps (worst case)')
ax.legend()

ax = axes[1]
ax.semilogx(n_values, linear_steps, color='red',       label='O(n)')
ax.semilogx(n_values, log_steps,    color='steelblue', label='O(log₂n)', linewidth=2)
ax.set_title('Same Plot on Semilog Scale')
ax.set_xlabel('Array size n (log scale)')
ax.set_ylabel('Steps (worst case)')
ax.legend()

plt.tight_layout()
plt.show()

for n in [100, 1_000, 1_000_000, 1_000_000_000]:
    print(f"n={n:>12,}: linear={n:>12,} steps, binary search={int(np.log2(n)):>2} steps")

# --- Experiment 2: Where does log₂(n) equal n? ---
# Hypothesis: log₂(n) < n for all n > 1. Find the crossover numerically.
# Also: the boundary log₂(n!) vs n² via Stirling's approximation.
# Try changing: BASE

import numpy as np

BASE = 2.0  # <-- modify this (try 3.0, 10.0, np.e)

print(f"Comparing log_{BASE}(n) vs n:")
print(f"{'n':>10} {'log_b(n)':>12} {'n':>12} {'log<n':>8}")
print("-" * 46)
for n in [0.1, 0.5, 1, 2, 4, 10, 100, 1000]:
    logval = np.log(n) / np.log(BASE)
    print(f"{n:>10.3f} {logval:>12.4f} {n:>12.4f} {str(logval < n):>8}")

# For base > 1, log_b(n) < n for all n > 1. But log can exceed n for n < 1.
print(f"\nConclusion: log_{BASE}(n) < n for n > 1 (log grows much slower)")

# --- Experiment 3: Real-world decibels and pH — log scales in practice ---
# Hypothesis: Adding 10 dB doesn't add 10 to the power — it multiplies by 10.
# Try changing: POWER_WATTS to explore the dB scale

import numpy as np

P_REFERENCE = 1e-12  # reference power in watts (threshold of human hearing)

print("Sound Pressure Levels (dB = 10 × log₁₀(P / P_ref))")
print("-" * 50)
sounds = {
    'Threshold of hearing': 1e-12,
    'Whisper':              1e-10,
    'Normal conversation':  1e-6,
    'Vacuum cleaner':       1e-4,
    'Rock concert':         1e-1,
    'Jet engine (nearby)':  1e2,
}
for label, power in sounds.items():
    dB = 10 * np.log10(power / P_REFERENCE)
    print(f"{label:30s}: {power:.1e} W  →  {dB:6.1f} dB")

print("\n+10 dB means power × 10. +20 dB means power × 100. +3 dB ≈ power × 2.")
print("This compression is only legible because of the logarithm.")

# --- Mini Project: Log-linearization as a diagnostic tool ---
#
# Problem:
#   Given a dataset of (time, value) pairs, determine whether the growth is
#   linear, polynomial (quadratic), or exponential — using logarithms.
#
#   Key insight:
#     - Exponential:  log(y) is linear in x         → plot log(y) vs x
#     - Power law:    log(y) is linear in log(x)    → plot log(y) vs log(x)
#     - Linear:       y is linear in x (no transform needed)
#
# Dataset: Three synthetic series with noise.

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-v0_8-whitegrid')
np.random.seed(7)

t = np.linspace(1, 20, 50)
noise = lambda scale: np.random.normal(0, scale, len(t))

# Three datasets
y_linear = 3 * t + 5           + noise(3)
y_power  = 2 * t**2.3          + noise(t**2.3 * 0.05)
y_exp    = 5 * np.exp(0.2 * t) + noise(5 * np.exp(0.2 * t) * 0.05)

# TODO 1: For each dataset, plot:
#   (a) y vs t (raw)
#   (b) log(y) vs t (semi-log: linearizes exponential)
#   (c) log(y) vs log(t) (log-log: linearizes power law)
# Identify which transform makes each dataset look most linear.

fig, axes = plt.subplots(3, 3, figsize=(14, 11))
datasets = [('Linear: 3t+5',  y_linear, 'green'),
            ('Power: 2t^2.3', y_power,  'blue'),
            ('Exponential: 5e^0.2t', y_exp, 'red')]

for row, (title, y, color) in enumerate(datasets):
    # Raw
    axes[row, 0].scatter(t, y, s=15, color=color, alpha=0.7)
    axes[row, 0].set_title(f'{title}\nRaw: y vs t')
    axes[row, 0].set_xlabel('t')
    axes[row, 0].set_ylabel('y')

    # Semi-log: log(y) vs t
    axes[row, 1].scatter(t, np.log(np.maximum(y, 1e-9)), s=15, color=color, alpha=0.7)
    axes[row, 1].set_title(f'Semi-log: log(y) vs t')
    axes[row, 1].set_xlabel('t')
    axes[row, 1].set_ylabel('log(y)')

    # Log-log: log(y) vs log(t)
    axes[row, 2].scatter(np.log(t), np.log(np.maximum(y, 1e-9)), s=15, color=color, alpha=0.7)
    axes[row, 2].set_title(f'Log-log: log(y) vs log(t)')
    axes[row, 2].set_xlabel('log(t)')
    axes[row, 2].set_ylabel('log(y)')

plt.suptitle('Identifying Growth Type via Log Transformations', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

# TODO 2: For the exponential dataset, fit a line to log(y) vs t
# and recover the growth rate k and initial value P0.
log_y_exp = np.log(y_exp)
coeffs = np.polyfit(t, log_y_exp, 1)
k_est  = coeffs[0]
P0_est = np.exp(coeffs[1])
print(f"Exponential fit from log-linearization:")
print(f"  Estimated k:  {k_est:.4f}  (true: 0.2000)")
print(f"  Estimated P₀: {P0_est:.4f}  (true: 5.0000)")