Elementary Number Theory with SageMath

# Setup: Import standard Python libraries for number theory
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from scipy import *
import pandas as pd
from collections import Counter
import math

print("Elementary Number Theory with Standard Python")
print("=" * 45)
print("Ready to explore the fascinating world of integers!")
print("Using SymPy for advanced number theory computations.")
print()
print("Fun Fact: The integers were called 'God's gift to mathematics' by Kronecker!")

# Chapter 1: Prime Number Fundamentals
import sympy as sp

print("PRIME NUMBER EXPLORATION\n")

# Test primality of various numbers
test_numbers = [2, 17, 51, 97, 221, 1009]

print("Prime Testing:")
print("Number | Prime? | Factorization")
print("-" * 35)

for n in test_numbers:
    is_prime_result = sp.isprime(n)
    if not is_prime_result:
        factorization = sp.factorint(n)
        factor_str = " × ".join([f"{p}^{e}" if e > 1 else str(p) for p, e in factorization.items()])
    else:
        factor_str = "prime"
    status = "✓" if is_prime_result else "✗"
    print(f"{int(n):6d} |   {status}    | {factor_str}")

print("\n💡 Insight: Composite numbers reveal their structure through factorization!")
print()

# Generate first primes using SymPy
first_primes = list(sp.primerange(2, 100))
print(f"First 20 primes: {first_primes[:20]}")
print(f"Total primes under 100: {len(first_primes)}")

# Prime counting function π(x)
print("\nPrime Counting Function π(x):")
print("   x   | π(x) | Density")
print("-" * 25)

for x_val in [10, 50, 100, 500, 1000]:
    xv = int(x_val)
    pi_x = len(list(sp.primerange(2, xv + 1)))
    density = 100.0 * pi_x / xv  # force Python float
    print(f"{xv:6d} | {pi_x:4d} | {float(density):5.1f}%")

print("\nObservation: Prime density decreases as numbers get larger!")

# Prime Factorization Deep Dive
print("PRIME FACTORIZATION WORKSHOP\n")

# Demonstrate unique factorization
numbers_to_factor = [60, 315, 1001, 2310, 9999]

print("Unique Prime Factorizations:")
for n in numbers_to_factor:
    factorization = factor(n)
    
    # Verify by multiplication
    product_check = prod(p^e for p, e in factorization)
    
    print(f"{n:4d} = {factorization} {'✓' if product_check == n else '✗'}")

print()

# Divisor function analysis
print("Divisors and the σ functions:")
print("n    | d(n) | σ(n) | Divisors")
print("-" * 45)

for n in [12, 24, 30, 60]:
    divs = divisors(n)
    num_divisors = len(divs)  # d(n)
    sum_divisors = sum(divs)  # σ(n)
    divs_str = ', '.join(map(str, divs))
    
    print(f"{n:4d} | {num_divisors:4d} | {sum_divisors:4d} | {divs_str}")

print()

# Perfect numbers - where σ(n) = 2n
def check_perfect_numbers(limit):
    perfect_nums = []
    
    for n in range(1, limit):
        proper_divisors = [d for d in divisors(n) if d < n]
        if sum(proper_divisors) == n:
            perfect_nums.append(n)
    
    return perfect_nums

perfect_numbers = check_perfect_numbers(10000)
print(f"Perfect Numbers (up to 10,000): {perfect_numbers}")

if perfect_numbers:
    n = perfect_numbers[0]  # First perfect number
    proper_divs = [d for d in divisors(n) if d < n]
    print(f"\nExample: {n} = {' + '.join(map(str, proper_divs))} = {sum(proper_divs)}")

print("\nPerfect numbers are rare gems in the integer landscape!")

# Sieve of Eratosthenes - Ancient Algorithm, Modern Implementation
print("SIEVE OF ERATOSTHENES (circa 250 BCE)\n")

def sieve_of_eratosthenes_demo(n):
    """
    Implement and demonstrate the Sieve of Eratosthenes
    """
    print(f"Finding all primes up to {n} using the Sieve:")
    
    # Initialize: all numbers from 2 to n are potentially prime
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False  # 0 and 1 are not prime
    
    primes_found = []
    
    for p in range(2, int(n**0.5) + 1):
        if is_prime[p]:
            print(f"\nSieving with prime {p}:")
            # Mark all multiples of p as composite
            multiples_marked = 0
            for i in range(p * p, n + 1, p):
                if is_prime[i]:
                    is_prime[i] = False
                    multiples_marked += 1
            print(f"  Marked {multiples_marked} multiples of {p}")
    
    # Collect all remaining primes
    primes = [i for i in range(2, n + 1) if is_prime[i]]
    
    print(f"\nPrimes found: {primes}")
    print(f"Total count: {len(primes)}")
    
    return primes

# Demonstrate with small example
sieve_primes = sieve_of_eratosthenes_demo(30)

# Compare with SageMath's built-in prime generation
sagemath_primes = [p for p in range(2, 31) if is_prime(p)]
print(f"\nVerification with SageMath: {sagemath_primes}")
print(f"Results match: {sieve_primes == sagemath_primes}")

print("\n⚡ The Sieve is incredibly efficient for finding many primes at once!")
print("   Time complexity: O(n log log n) - much faster than testing each number individually.")

# Chapter 2: The Euclidean Algorithm
print("THE EUCLIDEAN ALGORITHM\n")

def euclidean_algorithm_detailed(a, b, show_steps=True):
    """
    Perform Euclidean algorithm with detailed step-by-step output
    """
    if show_steps:
        print(f"Finding gcd({a}, {b}):")
        print(f"{'Step':>4} | {'Equation':>20} | {'Remainder':>10}")
        print("-" * 40)
    
    original_a, original_b = a, b
    step = 1
    
    while b != 0:
        quotient = a // b
        remainder = a % b
        
        if show_steps:
            equation = f"{a} = {b}×{quotient} + {remainder}"
            print(f"{step:4d} | {equation:>20} | {remainder:>10}")
        
        a, b = b, remainder
        step += 1
    
    if show_steps:
        print(f"\nResult: gcd({original_a}, {original_b}) = {a}")
    
    return a

# Demonstrate with classic examples
test_pairs = [(1071, 462), (252, 198), (48, 18)]

for a, b in test_pairs:
    our_result = euclidean_algorithm_detailed(a, b)
    sage_result = gcd(a, b)
    print(f"Verification: SageMath gcd({a}, {b}) = {sage_result} {'✓' if our_result == sage_result else '✗'}")
    print("=" * 50)

# GCD properties and applications
print("\nGCD PROPERTIES AND RELATIONSHIPS\n")

print("Key Identity: gcd(a,b) × lcm(a,b) = a × b")
print("Verification:")

for a, b in [(12, 18), (15, 25), (7, 11)]:
    gcd_ab = gcd(a, b)
    lcm_ab = lcm(a, b)
    product1 = gcd_ab * lcm_ab
    product2 = a * b
    
    print(f"  gcd({a},{b}) × lcm({a},{b}) = {gcd_ab} × {lcm_ab} = {product1}")
    print(f"  {a} × {b} = {product2} {'✓' if product1 == product2 else '✗'}")
    print()

# Coprimality (relatively prime)
print("Coprime Numbers (gcd = 1):")
test_coprime = [(8, 15), (21, 22), (25, 49), (17, 19)]

for a, b in test_coprime:
    gcd_result = gcd(a, b)
    coprime_status = "coprime" if gcd_result == 1 else "not coprime"
    print(f"  gcd({a:2d}, {b:2d}) = {gcd_result} → {coprime_status}")

# Extended Euclidean Algorithm - Bézout's Identity
print("🔧 EXTENDED EUCLIDEAN ALGORITHM\n")
print("Finds integers x, y such that ax + by = gcd(a, b)")
print("This is called Bézout's Identity!")
print()

def extended_euclidean_demo(a, b):
    """
    Demonstrate the extended Euclidean algorithm
    """
    print(f"Finding x, y such that {a}x + {b}y = gcd({a}, {b})")
    
    # Use SageMath's extended GCD
    g, x, y = xgcd(a, b)
    
    print(f"\nResult:")
    print(f"  gcd({a}, {b}) = {g}")
    print(f"  Bézout coefficients: x = {x}, y = {y}")
    print(f"  Verification: {a}({x}) + {b}({y}) = {a*x + b*y}")
    
    # Check the identity
    check = a*x + b*y == g
    print(f"  Identity holds: {check} {'✓' if check else '✗'}")
    
    return g, x, y

# Examples of extended Euclidean algorithm
examples = [(240, 46), (161, 28), (1071, 462)]

for a, b in examples:
    extended_euclidean_demo(a, b)
    print("-" * 40)

print("\nApplications of Extended Euclidean Algorithm:")
print("  • Solving linear Diophantine equations")
print("  • Computing modular inverses")
print("  • Cryptographic algorithms (RSA key generation)")
print("  • Simplifying fractions")

# Practical example: Modular inverse
print("\nComputing Modular Inverse:")
print("To find a⁻¹ mod m, we solve ax + my = 1")

a, m = 7, 26
g, x, y = xgcd(a, m)

if g == 1:
    inverse = x % m
    print(f"\nFinding {a}⁻¹ mod {m}:")
    print(f"  {a}({x}) + {m}({y}) = {g}")
    print(f"  Therefore, {a}⁻¹ ≡ {inverse} (mod {m})")
    
    # Verify
    verification = (a * inverse) % m
    print(f"  Check: {a} × {inverse} ≡ {verification} (mod {m}) {'✓' if verification == 1 else '✗'}")
else:
    print(f"  No inverse exists since gcd({a}, {m}) = {g} ≠ 1")

# Chapter 3: Modular Arithmetic Fundamentals
print(" MODULAR ARITHMETIC FUNDAMENTALS\n")

# Basic modular arithmetic operations
def create_modular_table(operation, mod):
    """
    Create and display modular arithmetic table
    """
    print(f"Modular {operation} table (mod {mod}):")
    
    # Header
    print(f"  {operation[0]}", end="")
    for j in range(mod):
        print(f"{j:3d}", end="")
    print()
    
    # Rows
    for i in range(mod):
        print(f"{i:3d}", end="")
        for j in range(mod):
            if operation == "addition":
                result = (i + j) % mod
            elif operation == "multiplication":
                result = (i * j) % mod
            print(f"{result:3d}", end="")
        print()
    print()

# Show modular arithmetic tables
create_modular_table("addition", 7)
create_modular_table("multiplication", 7)

# Modular exponentiation patterns
print("Powers modulo 7:")
print("Base | Powers mod 7")
print("-" * 25)

mod = 7
for base in range(1, mod):
    powers = []
    for exp in range(1, mod):
        power_result = power_mod(base, exp, mod)
        powers.append(str(power_result))
    
    powers_str = ', '.join(powers)
    print(f"{base:4d} | {powers_str}")

print("\n🎯 Notice: Some bases generate all non-zero remainders!")
print("   These are called 'primitive roots' or 'generators'.")

# Find primitive roots
def find_primitive_roots(p):
    """Find primitive roots modulo prime p"""
    if not is_prime(p):
        return []
    
    primitive_roots = []
    phi_p = p - 1  # Euler's totient of prime p
    
    for g in range(2, p):
        # Check if g generates all non-zero elements
        generated = set()
        for k in range(1, p):
            generated.add(power_mod(g, k, p))
        
        if len(generated) == phi_p:
            primitive_roots.append(g)
    
    return primitive_roots

p = 7
roots = find_primitive_roots(p)
print(f"\nPrimitive roots modulo {p}: {roots}")
print(f"These bases generate all {p-1} non-zero remainders when raised to different powers.")

# Fast Modular Exponentiation
print("⚡ FAST MODULAR EXPONENTIATION\n")

def modular_exp_binary_demo(base, exp, mod):
    """
    Demonstrate binary exponentiation for modular arithmetic
    """
    print(f"Computing {base}^{exp} mod {mod} using binary method:")
    
    # Show why this is necessary
    if exp > 20:
        digits = int(exp * log(base, 10)) if base > 1 else 1
        print(f"Direct computation of {base}^{exp} would have ~{digits} digits!")
    
    # Convert exponent to binary
    binary_exp = bin(exp)[2:]  # Remove '0b' prefix
    print(f"\nExponent {exp} in binary: {binary_exp}")
    
    # Binary exponentiation algorithm
    result = 1
    base_power = base % mod
    
    print(f"\nStep-by-step computation:")
    print(f"Bit | Power | Result")
    print("-" * 25)
    
    for i, bit in enumerate(reversed(binary_exp)):
        if bit == '1':
            result = (result * base_power) % mod
            print(f" {i}  |  {base_power:4d}  |  {result:4d}")
        
        if i < len(binary_exp) - 1:  # Don't square on last iteration
            base_power = (base_power * base_power) % mod
    
    print(f"\nFinal result: {base}^{exp} ≡ {result} (mod {mod})")
    return result

# Examples of fast modular exponentiation
examples = [(3, 13, 7), (2, 100, 17), (123, 456, 789)]

for base, exp, mod in examples:
    our_result = modular_exp_binary_demo(base, exp, mod)
    sage_result = power_mod(base, exp, mod)
    print(f"SageMath verification: {sage_result} {'✓' if our_result == sage_result else '✗'}")
    print("=" * 50)

print("\n Why Fast Exponentiation Matters:")
print("  • RSA encryption uses 1024+ bit exponents")
print("  • Naive method would take longer than age of universe")
print("  • Binary method reduces O(n) to O(log n) multiplications")
print("  • Essential for modern cryptography")

# Fermat's Little Theorem
print("FERMAT'S LITTLE THEOREM\n")
print("If p is prime and gcd(a, p) = 1, then a^(p-1) ≡ 1 (mod p)")
print("This is one of the most beautiful results in number theory!")
print()

# Verify Fermat's Little Theorem
def verify_fermat_little_theorem(p):
    """Verify FLT for all valid bases modulo prime p"""
    if not is_prime(p):
        print(f"Error: {p} is not prime")
        return False
    
    print(f"Testing Fermat's Little Theorem for prime p = {p}:")
    print(f"{'Base a':>6} | {'a^{}'.format(p-1):>8} | {'Result mod {}'.format(p):>12}")
    print("-" * 35)
    
    all_pass = True
    for a in range(1, min(p, 10)):  # Test first few bases
        if gcd(a, p) == 1:  # Should always be true for prime p
            result = power_mod(a, p-1, p)
            status = "✓" if result == 1 else "✗"
            print(f"{a:6d} | {p-1:8d} | {result:8d} {status:>3}")
            
            if result != 1:
                all_pass = False
    
    print(f"\nAll tests pass: {all_pass}")
    return all_pass

# Test with several primes
test_primes = [5, 7, 11, 13]

for p in test_primes:
    verify_fermat_little_theorem(p)
    print()

# Fermat primality test and pseudoprimes
print("FERMAT PRIMALITY TEST\n")
print("Idea: If n is composite, then a^(n-1) ≢ 1 (mod n) for most bases a")
print("But some composites fool us - these are Fermat pseudoprimes!")
print()

def fermat_test(n, base=2):
    """Fermat primality test with given base"""
    if gcd(base, n) != 1:
        return False
    return power_mod(base, n-1, n) == 1

# Test some Carmichael numbers (fool all coprime bases)
carmichael_numbers = [561, 1105, 1729, 2465, 2821]

print("Testing Carmichael numbers (composite but pass Fermat test):")
print(f"{'Number':>6} | {'Prime?':>7} | {'Fermat(base=2)':>13} | {'Pseudoprime?':>12}")
print("-" * 45)

for n in carmichael_numbers:
    is_prime_actual = is_prime(n)
    fermat_result = fermat_test(n, 2)
    is_pseudoprime = fermat_result and not is_prime_actual
    
    print(f"{n:6d} | {str(is_prime_actual):>7} | {str(fermat_result):>13} | {str(is_pseudoprime):>12}")

print("\n⚠Lesson: Fermat's test is probabilistic, not definitive!")
print("    Modern primality tests use more sophisticated methods.")

# Chapter 4: Linear Congruences
print("SOLVING LINEAR CONGRUENCES\n")
print("Form: ax ≡ b (mod m)")
print("Goal: Find all values of x that satisfy the congruence")
print()

def solve_linear_congruence_detailed(a, b, m):
    """
    Solve ax ≡ b (mod m) with detailed explanation
    """
    print(f"Solving {a}x ≡ {b} (mod {m})")
    print("-" * 30)
    
    # Step 1: Check if solution exists
    g = gcd(a, m)
    print(f"Step 1: gcd({a}, {m}) = {g}")
    
    if b % g != 0:
        print(f"No solution: {g} does not divide {b}")
        return []
    
    print(f"Solution exists: {g} divides {b}")
    
    # Step 2: Reduce the congruence
    a_reduced = a // g
    b_reduced = b // g  
    m_reduced = m // g
    
    print(f"Step 2: Reduce by dividing by {g}")
    print(f"        {a_reduced}x ≡ {b_reduced} (mod {m_reduced})")
    
    # Step 3: Find modular inverse
    gcd_result, inv_a, _ = xgcd(a_reduced, m_reduced)
    inv_a = inv_a % m_reduced
    
    print(f"Step 3: Find inverse of {a_reduced} mod {m_reduced}")
    print(f"        {a_reduced}^(-1) ≡ {inv_a} (mod {m_reduced})")
    
    # Step 4: Solve for x
    x0 = (inv_a * b_reduced) % m_reduced
    print(f"Step 4: x ≡ {inv_a} × {b_reduced} ≡ {x0} (mod {m_reduced})")
    
    # Step 5: Generate all solutions
    solutions = [(x0 + i * m_reduced) % m for i in range(g)]
    solutions.sort()
    
    print(f"Step 5: All solutions modulo {m}:")
    for i, sol in enumerate(solutions):
        print(f"        x ≡ {sol} (mod {m})")
    
    # Verification
    print(f"Verification:")
    for sol in solutions:
        result = (a * sol) % m
        status = "✓" if result == b % m else "✗"
        print(f"  {a} × {sol} ≡ {result} ≡ {b} (mod {m}) {status}")
    
    return solutions

# Examples of linear congruences
examples = [
    (3, 7, 10),   # Unique solution
    (6, 9, 15),   # Multiple solutions  
    (4, 6, 10),   # No solution
    (5, 3, 7)     # Coprime case
]

for a, b, m in examples:
    solve_linear_congruence_detailed(a, b, m)
    print("=" * 60)

print("\n💡 Key Insight: Linear congruences generalize linear equations to modular arithmetic!")

# Chinese Remainder Theorem
print("CHINESE REMAINDER THEOREM\n")
print("Ancient Chinese text: 'There are certain things whose number is unknown.'")
print("'When divided by 3, the remainder is 2.'")
print("'When divided by 5, the remainder is 3.'")
print("'When divided by 7, the remainder is 2.'")
print("'What will be the number of things?'")
print()

def chinese_remainder_theorem_detailed(remainders, moduli):
    """
    Solve system of congruences with detailed steps
    """
    n = len(remainders)
    
    print("System of congruences:")
    for i, (a, m) in enumerate(zip(remainders, moduli)):
        print(f"  x ≡ {a} (mod {m})")
    print()
    
    # Check pairwise coprimality
    print("Step 1: Check that moduli are pairwise coprime")
    for i in range(n):
        for j in range(i+1, n):
            g = gcd(moduli[i], moduli[j])
            if g != 1:
                print(f"Error: gcd({moduli[i]}, {moduli[j]}) = {g} ≠ 1")
                return None
            print(f"  gcd({moduli[i]}, {moduli[j]}) = {g} ✓")
    
    # Compute M
    M = 1
    for m in moduli:
        M *= m
    
    print(f"\nStep 2: Compute M = {' × '.join(map(str, moduli))} = {M}")
    
    # CRT construction
    print(f"\nStep 3: CRT construction")
    x = 0
    
    for i, (a, m) in enumerate(zip(remainders, moduli)):
        Mi = M // m
        print(f"  For x ≡ {a} (mod {m}):")
        print(f"    M_{i+1} = M/m_{i+1} = {M}/{m} = {Mi}")
        
        # Find modular inverse
        g, yi, _ = xgcd(Mi, m)
        yi = yi % m
        print(f"    y_{i+1} = M_{i+1}^(-1) mod m_{i+1} = {yi}")
        
        # Verify inverse
        check = (Mi * yi) % m
        print(f"    Verify: {Mi} × {yi} ≡ {check} (mod {m}) {'✓' if check == 1 else '✗'}")
        
        # Add contribution
        contribution = a * Mi * yi
        x += contribution
        print(f"    Contribution: {a} × {Mi} × {yi} = {contribution}")
        print()
    
    # Final answer
    x = x % M
    print(f"Step 4: Final solution")
    print(f"  x ≡ {x} (mod {M})")
    
    # Verification
    print(f"\nStep 5: Verification")
    all_correct = True
    for a, m in zip(remainders, moduli):
        remainder = x % m
        correct = remainder == a
        status = "✓" if correct else "✗"
        print(f"  {x} ≡ {remainder} ≡ {a} (mod {m}) {status}")
        if not correct:
            all_correct = False
    
    print(f"\nAll verifications pass: {all_correct}")
    return x

# Solve the classic Chinese problem
print("CLASSIC CHINESE PROBLEM:")
remainders_classic = [2, 3, 2]
moduli_classic = [3, 5, 7]
solution_classic = chinese_remainder_theorem_detailed(remainders_classic, moduli_classic)

print("=" * 70)

# Modern example
print("\nMODERN EXAMPLE:")
remainders_modern = [1, 2, 3, 4]
moduli_modern = [5, 7, 9, 11]
solution_modern = chinese_remainder_theorem_detailed(remainders_modern, moduli_modern)

# Verify with SageMath
print("\n🔧 SAGEMATH VERIFICATION:")
sage_solution1 = crt(remainders_classic, moduli_classic)
sage_solution2 = crt(remainders_modern, moduli_modern)

print(f"Classic problem: Our solution = {solution_classic}, SageMath = {sage_solution1}")
print(f"Modern example: Our solution = {solution_modern}, SageMath = {sage_solution2}")

match1 = solution_classic == sage_solution1
match2 = solution_modern == sage_solution2
print(f"Results match: {match1 and match2} {'✓' if match1 and match2 else '✗'}")

# Chapter 5: Euler's Totient Function
print("φ EULER'S TOTIENT FUNCTION\n")
print("φ(n) = number of integers from 1 to n that are coprime to n")
print("This function is fundamental to RSA cryptography!")
print()

# Compute φ(n) and show coprime numbers
print("Totient values with coprime numbers:")
print(f"{'n':>3} | {'φ(n)':>4} | Coprime numbers")
print("-" * 50)

for n in range(1, 16):
    phi_n = euler_phi(n)
    coprimes = [i for i in range(1, n+1) if gcd(i, n) == 1]
    coprime_str = ', '.join(map(str, coprimes))
    
    print(f"{n:3d} | {phi_n:4d} | {coprime_str}")

print()

# Properties of Euler's totient function
print("🔍 PROPERTIES OF φ(n):\n")

# Property 1: φ(p) = p-1 for prime p
print("1. For prime p: φ(p) = p - 1")
small_primes = [p for p in range(2, 20) if is_prime(p)]
for p in small_primes[:6]:
    phi_p = euler_phi(p)
    print(f"   φ({p:2d}) = {phi_p:2d} = {p} - 1")

# Property 2: φ(p^k) formula
print("\n2. For prime power p^k: φ(p^k) = p^(k-1)(p-1)")
prime_power_examples = [(2, 3), (3, 2), (5, 2)]
for p, k in prime_power_examples:
    n = p**k
    phi_n = euler_phi(n)
    formula_result = p**(k-1) * (p-1)
    print(f"   φ({p}^{k}) = φ({n:2d}) = {phi_n:2d} = {p}^{k-1} × ({p}-1) = {formula_result}")

# Property 3: Multiplicativity
print("\n3. φ is multiplicative: if gcd(m,n) = 1, then φ(mn) = φ(m)φ(n)")
mult_examples = [(3, 5), (4, 9), (7, 8)]
for m, n in mult_examples:
    if gcd(m, n) == 1:
        phi_m = euler_phi(m)
        phi_n = euler_phi(n)
        phi_mn = euler_phi(m * n)
        product = phi_m * phi_n
        status = "✓" if phi_mn == product else "✗"
        print(f"   φ({m} × {n}) = φ({m*n:2d}) = {phi_mn:2d}, φ({m}) × φ({n}) = {phi_m} × {phi_n} = {product} {status}")

# General formula using prime factorization
print("\n4. General formula: φ(n) = n × ∏(1 - 1/p) for prime factors p")

def show_totient_formula(n):
    factorization = factor(n)              # e.g. 12 -> 2^2 * 3
    factors_list = list(factorization)     # [(2,2), (3,1)]
    
    print(f"   n = {n}, factorization = {factorization}")
    
    # Apply $φ(n)=n\prod_{p|n}\left(1-\frac{1}{p}\right)$
    result = n
    step_terms = [str(n)]
    for p, e in factors_list:
        result = result * (p - 1) // p
        step_terms.append(f"(1 - 1/{p})")
    
    formula_str = " × ".join(step_terms)
    print(f"   φ({n}) = {formula_str} = {result}")

for n in [12, 30, 60]:
    show_totient_formula(n)

print("\n🔐 Application: In RSA cryptography, we use φ(pq) = (p-1)(q-1) where p, q are large primes!")

# Other Important Number Theory Functions
print("OTHER NUMBER THEORY FUNCTIONS\n")

# Divisor functions
print("Divisor Functions:")
print(f"{'n':>3} | {'d(n)':>4} | {'σ(n)':>5} | Divisors")
print("-" * 45)

for n in range(1, 21):
    divs = divisors(n)
    d_n = len(divs)  # τ(n) or d(n) - number of divisors
    sigma_n = sum(divs)  # σ(n) - sum of divisors
    divs_str = ', '.join(map(str, divs))
    if len(divs_str) > 20:
        divs_str = divs_str[:17] + "..."
    
    print(f"{n:3d} | {d_n:4d} | {sigma_n:5d} | {divs_str}")

print()

# Perfect and related numbers
print("SPECIAL NUMBERS BASED ON DIVISOR FUNCTIONS:\n")

def classify_by_divisor_sum(limit):
    """Classify numbers by relationship between n and sum of proper divisors"""
    perfect = []
    abundant = []
    deficient = []
    
    for n in range(2, limit):
        proper_divisors = [d for d in divisors(n) if d < n]
        divisor_sum = sum(proper_divisors)
        
        if divisor_sum == n:
            perfect.append(n)
        elif divisor_sum > n:
            abundant.append(n)
        else:
            deficient.append(n)
    
    return perfect, abundant[:10], deficient[:15]  # Limit output

perfect, abundant, deficient = classify_by_divisor_sum(100)

print(f"Perfect numbers (σ(n) = 2n): {perfect}")
if perfect:
    n = perfect[0]
    proper_divs = [d for d in divisors(n) if d < n]
    print(f"  Example: {n} = {' + '.join(map(str, proper_divs))} = {sum(proper_divs)}")

print(f"\nAbundant numbers (σ(n) > 2n): {abundant}")
print(f"Deficient numbers (σ(n) < 2n): {deficient}")

print()

# Möbius function
print("μ MÖBIUS FUNCTION:\n")
print("μ(n) = 0 if n has a squared prime factor")
print("μ(n) = (-1)^k if n is a product of k distinct primes")
print("μ(1) = 1")
print()

print(f"{'n':>3} | {'μ(n)':>4} | {'Prime factorization':>20} | {'Square-free?':>12}")
print("-" * 55)

for n in range(1, 21):
    mu_n = moebius(n)
    factorization = factor(n)
    
# Check if square-free (no repeated prime factors)
is_square_free = all(exp == 1 for _, exp in factorization) if n > 1 else True
    
print(f"{n:3d} | {mu_n:4d} | {str(factorization):>20} | {str(is_square_free):>12}")

print("\n💡 Key insight: μ(n) ≠ 0 exactly when n is square-free!")

# Möbius inversion formula example
print("\nMöbius Inversion Formula:")
print("If f(n) = Σ g(d) for d|n, then g(n) = Σ μ(d)f(n/d) for d|n")
print("This powerful tool connects many number theory functions!")

# Example: φ(n) in terms of divisor sums
def verify_totient_mobius_formula(n):
    """Verify φ(n) = Σ μ(d) × (n/d) for d|n"""
    divs = divisors(n)
    mobius_sum = sum(moebius(d) * (n // d) for d in divs)
    actual_phi = euler_phi(n)
    
    print(f"φ({n}) = {actual_phi}, Möbius formula gives: {mobius_sum} {'✓' if mobius_sum == actual_phi else '✗'}")

print("\nVerifying φ(n) = Σ μ(d) × (n/d) for divisors d of n:")
for n in [6, 12, 15, 20]:
    verify_totient_mobius_formula(n)