Chapter 41 — Exponents and Powers

# --- Visualization: Exponent landscape — how b^n behaves across bases and exponents ---
import numpy as np
import matplotlib.pyplot as plt

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

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# --- Panel 1: b^n for fixed base b, varying exponent n ---
ax = axes[0]
exponents = np.arange(0, 11)
for base in [2, 3, 10]:
    values = base ** exponents
    ax.plot(exponents, values, marker='o', label=f'base = {base}')
ax.set_title('b^n: Fixed Base, Varying Exponent')
ax.set_xlabel('Exponent n')
ax.set_ylabel('b^n')
ax.legend()

# --- Panel 2: b^n for fixed exponent n, varying base b ---
ax = axes[1]
bases = np.linspace(0.1, 4, 200)
for exp in [1, 2, 3, 4]:
    ax.plot(bases, bases ** exp, label=f'n = {exp}')
ax.set_title('b^n: Fixed Exponent, Varying Base')
ax.set_xlabel('Base b')
ax.set_ylabel('b^n')
ax.legend()
ax.set_ylim(0, 20)

# --- Panel 3: Negative and fractional exponents ---
ax = axes[2]
exps = np.linspace(-3, 3, 300)
for base in [2, 0.5]:
    ax.plot(exps, base ** exps, label=f'base = {base}')
ax.axhline(1, color='gray', linestyle='--', linewidth=0.8, label='y = 1 (any b^0)')
ax.set_title('b^n: Fractional and Negative Exponents')
ax.set_xlabel('Exponent n (real-valued)')
ax.set_ylabel('b^n')
ax.legend()

plt.suptitle('Exponentiation: Landscape Across Bases and Exponents', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

# --- Visualization: The geometric interpretation (squares and cubes) ---
import matplotlib.patches as patches

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

BASE = 3  # b = 3

# b^1: line
ax = axes[0]
ax.plot([0, BASE], [0.5, 0.5], 'b-', linewidth=6)
ax.set_xlim(-0.5, BASE + 0.5)
ax.set_ylim(0, 1)
ax.set_title(f'b^1 = {BASE}^1 = {BASE}\n(length)')
ax.set_xlabel('units')
ax.set_yticks([])

# b^2: square
ax = axes[1]
square = patches.Rectangle((0, 0), BASE, BASE, linewidth=1.5,
                             edgecolor='blue', facecolor='lightblue', alpha=0.5)
ax.add_patch(square)
ax.set_xlim(-0.5, BASE + 0.5)
ax.set_ylim(-0.5, BASE + 0.5)
ax.set_aspect('equal')
ax.set_title(f'b^2 = {BASE}^2 = {BASE**2}\n(area)')
ax.set_xlabel('units')
ax.set_ylabel('units')

# b^n: value plot showing geometric growth
ax = axes[2]
ns = np.arange(0, 7)
vals = BASE ** ns
ax.bar(ns, vals, color='steelblue', alpha=0.7)
for i, v in zip(ns, vals):
    ax.text(i, v + 0.5, str(v), ha='center', fontsize=9)
ax.set_title(f'b^n for b={BASE}, n=0..6')
ax.set_xlabel('Exponent n')
ax.set_ylabel('Value')

plt.tight_layout()
plt.show()

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

b, c = 3.0, 2.0
m, n = 4.0, 3.0

print("Laws of Exponents — Numerical Verification")
print("-" * 50)

# Law 1: b^m × b^n = b^(m+n)
lhs = b**m * b**n
rhs = b**(m + n)
print(f"1. b^m × b^n = b^(m+n):  {lhs:.4f} == {rhs:.4f}  {'✓' if np.isclose(lhs, rhs) else '✗'}")

# Law 2: b^m / b^n = b^(m-n)
lhs = b**m / b**n
rhs = b**(m - n)
print(f"2. b^m / b^n = b^(m-n):  {lhs:.4f} == {rhs:.4f}  {'✓' if np.isclose(lhs, rhs) else '✗'}")

# Law 3: (b^m)^n = b^(m×n)
lhs = (b**m)**n
rhs = b**(m * n)
print(f"3. (b^m)^n  = b^(m×n):   {lhs:.4f} == {rhs:.4f}  {'✓' if np.isclose(lhs, rhs) else '✗'}")

# Law 4: (b×c)^n = b^n × c^n
lhs = (b * c)**n
rhs = b**n * c**n
print(f"4. (b×c)^n = b^n × c^n: {lhs:.4f} == {rhs:.4f}  {'✓' if np.isclose(lhs, rhs) else '✗'}")

# Rational exponent: 8^(2/3)
val = 8 ** (2/3)
expected = 4.0
print(f"\nRational: 8^(2/3) = {val:.6f}  (expected {expected})  {'✓' if np.isclose(val, expected) else '✗'}")

# Real exponent: 2^π via exp(π × ln(2))
val_direct = 2 ** np.pi
val_via_exp = np.exp(np.pi * np.log(2))
print(f"Real: 2^π = {val_direct:.8f}  via exp(π×ln2) = {val_via_exp:.8f}  {'✓' if np.isclose(val_direct, val_via_exp) else '✗'}")

# --- Implementation: Naive integer exponentiation ---

def power_naive(base, exp):
    """
    Compute base^exp using repeated multiplication.
    Requires exactly (exp - 1) multiplications.

    Args:
        base: numeric base
        exp: non-negative integer exponent

    Returns:
        base raised to exp
    """
    if exp == 0:
        return 1
    result = 1
    for _ in range(exp):
        result *= base  # multiply exp times
    return result


# --- Implementation: Fast exponentiation (exponentiation by squaring) ---

def power_fast(base, exp):
    """
    Compute base^exp using exponentiation by squaring.
    Requires O(log exp) multiplications instead of O(exp).

    Key insight: 
        if exp is even:  base^exp = (base^2)^(exp/2)
        if exp is odd:   base^exp = base × base^(exp-1)

    Args:
        base: numeric base
        exp: non-negative integer exponent

    Returns:
        base raised to exp
    """
    if exp == 0:
        return 1
    if exp % 2 == 0:
        # Even: square the base, halve the exponent
        return power_fast(base * base, exp // 2)
    else:
        # Odd: pull out one factor, recurse on even
        return base * power_fast(base, exp - 1)


# Validate both implementations
test_cases = [(2, 10), (3, 7), (5, 0), (7, 1), (2, 32)]
print("base  exp  naive        fast         python")
print("-" * 55)
for b, e in test_cases:
    naive = power_naive(b, e)
    fast  = power_fast(b, e)
    ref   = b ** e
    match = '✓' if naive == fast == ref else '✗'
    print(f"{b:4d}  {e:3d}  {naive:12d}  {fast:12d}  {ref:12d}  {match}")

# --- Count multiplications: naive vs fast ---

def power_fast_counting(base, exp, count=None):
    """Same as power_fast but counts the number of multiplications."""
    if count is None:
        count = [0]
    if exp == 0:
        return 1, count[0]
    if exp % 2 == 0:
        count[0] += 1  # the squaring: base * base
        result, _ = power_fast_counting(base * base, exp // 2, count)
        return result, count[0]
    else:
        count[0] += 1  # the outer base * ...
        result, _ = power_fast_counting(base, exp - 1, count)
        return result, count[0]

import numpy as np
import matplotlib.pyplot as plt

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

exponents = np.arange(1, 65)
naive_mults = exponents - 1  # always exp-1 multiplications
fast_mults  = [power_fast_counting(2, int(e))[1] for e in exponents]
log2_bound  = np.log2(exponents) * 2  # theoretical O(log n) bound

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(exponents, naive_mults, label='Naive: O(n) multiplications', color='tomato')
ax.plot(exponents, fast_mults,  label='Fast: actual multiplications', color='steelblue', marker='.')
ax.plot(exponents, log2_bound,  label='2·log₂(n) bound', color='green', linestyle='--')
ax.set_title('Naive vs Fast Exponentiation: Multiplication Count')
ax.set_xlabel('Exponent n')
ax.set_ylabel('Number of multiplications')
ax.legend()
plt.tight_layout()
plt.show()

print(f"\nFor n=64: naive needs {64-1} multiplications, fast needs {power_fast_counting(2, 64)[1]}")

# --- Experiment 1: Where do the laws of exponents break down? ---
# Hypothesis: The laws assume b > 0 and real exponents.
# Negative bases with fractional exponents cause problems.
# Try changing: BASE to a negative number, EXP to a fraction

import numpy as np

BASE = -8   # <-- modify this (try -2, -1, 0)
EXP  = 1/3  # <-- modify this (try 2, -1, 0.5)

# Python's float arithmetic
try:
    result_float = BASE ** EXP
    print(f"Python float: ({BASE})^({EXP}) = {result_float}")
except Exception as e:
    print(f"Python float error: {e}")

# NumPy's complex-aware version
result_complex = complex(BASE) ** EXP
print(f"Complex result: ({BASE})^({EXP}) = {result_complex}")

# The real cube root of -8 should be -2 (since (-2)^3 = -8)
# Python cannot compute it as a float because x^(1/3) via exp(ln(x)/3)
# requires ln of a negative number — which is complex.
print(f"\nVerification: (-2)^3 = {(-2)**3}")
print(f"So the real cube root of -8 is -2, but Python gives: {BASE**EXP}")
print("This is why the laws require b > 0 for real-valued exponentiation.")

# --- Experiment 2: Power towers (iterated exponentiation) ---
# Hypothesis: a^(b^c) ≠ (a^b)^c in general — exponentiation is NOT associative.
# Try changing: A, B, C to explore when they happen to be equal

A = 2  # <-- modify this
B = 3  # <-- modify this
C = 2  # <-- modify this

# Right-associative: a^(b^c) — how power towers work by convention
right = A ** (B ** C)

# Left-associative: (a^b)^c — equivalent to a^(b×c) by Law 3
left = (A ** B) ** C

print(f"a={A}, b={B}, c={C}")
print(f"Right-associative: a^(b^c) = {A}^({B}^{C}) = {A}^{B**C} = {right}")
print(f"Left-associative: (a^b)^c  = ({A}^{B})^{C} = {A**B}^{C}  = {left}")
print(f"Equal? {right == left}")
print(f"\nNote: Law 3 says (a^b)^c = a^(b×c) = {A}^{B*C} = {A**(B*C)}")
print(f"But a^(b^c) = a^{B**C} = {right}  — very different for large exponents!")

# --- Experiment 3: Modular exponentiation (RSA preview) ---
# Hypothesis: (base^exp) mod m can be computed efficiently WITHOUT
# computing the full (astronomically large) base^exp first.
# This is the core operation in RSA encryption.

BASE = 7
EXP  = 1000    # <-- modify this (try 10_000, 100_000)
MOD  = 13      # <-- modify this (try other primes: 17, 29, 97)

# Python's built-in 3-argument pow is optimized for this
result = pow(BASE, EXP, MOD)  # computes (BASE^EXP) mod MOD efficiently

# Naive: compute the full number first (catastrophically slow for large EXP)
# result_naive = (BASE ** EXP) % MOD  # <-- try this for EXP=1000, then comment out

print(f"({BASE}^{EXP}) mod {MOD} = {result}")
print(f"Using Python's 3-arg pow: O(log {EXP}) ≈ {int(EXP.bit_length())} multiplications")
print(f"The full number {BASE}^{EXP} has ~{int(EXP * len(str(BASE)))} digits — never computed.")

# --- Mini Project: Smooth Numbers and the Cost of Representing Powers ---
#
# Problem:
#   A number is called "B-smooth" if all its prime factors are ≤ B.
#   Powers of small primes are the smoothest numbers — and they appear
#   constantly in algorithm analysis (2^n in binary search, 2^32 in
#   32-bit integers, etc.).
#
#   Your task: Build a tool that, given a number N, finds its representation
#   as a product of prime powers, and computes how many bits it needs.
#
# Dataset: Generated inline (integer powers of small primes)
# Task: Factor a number into prime powers, count its bit length, and
#       visualize the relationship between prime-power structure and bit length.

import numpy as np
import matplotlib.pyplot as plt

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

# --- Starter code: complete the TODOs ---

def prime_factorization(n):
    """
    Return the prime factorization of n as a dict {prime: exponent}.
    Example: prime_factorization(360) → {2: 3, 3: 2, 5: 1}
    """
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

def bit_length(n):
    """Number of bits needed to represent n in binary."""
    return n.bit_length()

# TODO 1: Generate 50 numbers that are powers of 2 and 3
#         (i.e., numbers of the form 2^a * 3^b for a,b in 0..10)
#         Hint: nested loop, collect unique values, sort
numbers = sorted(set(
    2**a * 3**b
    for a in range(12)
    for b in range(8)
))

# TODO 2: For each number, compute its bit length and factorization
bit_lengths = [bit_length(n) for n in numbers]
factorizations = [prime_factorization(n) for n in numbers]

# TODO 3: Plot bit_length vs log2(number)
# What do you expect this relationship to look like?
log2_values = np.log2(numbers)

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

ax = axes[0]
ax.scatter(log2_values, bit_lengths, alpha=0.6, color='steelblue')
ax.set_title('Bit Length vs log₂(n) for 2ᵃ × 3ᵇ numbers')
ax.set_xlabel('log₂(n)')
ax.set_ylabel('Bit length')

# TODO 4: Show which numbers have the most "prime power" structure
#         Color the scatter by the total sum of exponents (a + b)
total_exponents = [sum(f.values()) for f in factorizations]
sc = axes[1].scatter(np.array(numbers), bit_lengths, 
                     c=total_exponents, cmap='viridis', alpha=0.7)
plt.colorbar(sc, ax=axes[1], label='Total exponent sum (a+b)')
axes[1].set_xscale('log')
axes[1].set_title('Bit Length vs Number (log scale), colored by exponent sum')
axes[1].set_xlabel('n (log scale)')
axes[1].set_ylabel('Bit length')

plt.tight_layout()
plt.show()

# Inspect a few
print("Sample factorizations:")
for n in [1, 8, 12, 72, 576, 1152]:
    f = prime_factorization(n)
    expr = ' × '.join(f'{p}^{e}' for p, e in sorted(f.items()))
    print(f"  {n:6d} = {expr:20s}  ({bit_length(n)} bits)")