Elementary Number Theory with SageMath in CoCalc

# Chapter 6: Linear Diophantine Equations
print("LINEAR DIOPHANTINE EQUATIONS\n")
print("Form: ax + by = c")
print("Goal: Find integer solutions (x, y)")
print("Method: Use Extended Euclidean Algorithm")
print()

def solve_linear_diophantine_complete(a, b, c):
    """
    Completely solve linear Diophantine equation ax + by = c
    """
    print(f"Solving {a}x + {b}y = {c}")
    print("-" * 40)
    
    # Step 1: Check solvability
    g = gcd(a, b)
    print(f"Step 1: gcd({a}, {b}) = {g}")
    
    if c % g != 0:
        print(f"No integer solutions: {g} does not divide {c}")
        return None
    
    print(f"Solutions exist: {g} divides {c}")
    
    # Step 2: Find particular solution using extended Euclidean algorithm
    gcd_result, x0, y0 = xgcd(a, b)
    
    # Scale to match c
    scale = c // g
    x_particular = x0 * scale
    y_particular = y0 * scale
    
    print(f"Step 2: Extended GCD gives {a}×{x0} + {b}×{y0} = {gcd_result}")
    print(f"        Scaling by {scale}: {a}×{x_particular} + {b}×{y_particular} = {c}")
    
    # Verify particular solution
    check = a * x_particular + b * y_particular
    print(f"        Verification: {a}×{x_particular} + {b}×{y_particular} = {check} {'✓' if check == c else '✗'}")
    
    # Step 3: General solution
    b_coeff = b // g
    a_coeff = a // g
    
    print(f"Step 3: General solution (for any integer t):")
    print(f"        x = {x_particular} + {b_coeff}t")
    print(f"        y = {y_particular} - {a_coeff}t")
    
    # Step 4: Show specific solutions
    print(f"Step 4: Specific solutions:")
    print(f"{'t':>3} | {'x':>6} | {'y':>6} | {'Check':>8}")
    print("-" * 30)
    
    for t in range(-3, 4):
        x_t = x_particular + b_coeff * t
        y_t = y_particular - a_coeff * t
        check_t = a * x_t + b * y_t
        print(f"{t:3d} | {x_t:6d} | {y_t:6d} | {check_t:8d}")
    
    return (x_particular, y_particular, b_coeff, -a_coeff)

# Examples with different characteristics
examples = [
    (3, 5, 1),      # Coprime coefficients
    (6, 9, 15),     # Common factor that divides c
    (12, 18, 30),   # Larger common factor
    (4, 6, 5),      # No solution case
]

for a, b, c in examples:
    result = solve_linear_diophantine_complete(a, b, c)
    print("="*60)

print("\nApplications of Linear Diophantine Equations:")
print("  • Making change with different coin denominations")
print("  • Scheduling problems with time constraints")
print("  • Resource allocation in optimization")
print("  • Solving systems of congruences")

# Practical example: Coin change problem
print("\nPRACTICAL EXAMPLE: Coin Change")
print("How many ways can we make $1.00 using quarters (25¢) and dimes (10¢)?")
print("Equation: 25x + 10y = 100 or simplified: 5x + 2y = 20")

result = solve_linear_diophantine_complete(5, 2, 20)
if result:
    x_part, y_part, b_coeff, a_coeff = result
    print("\nNon-negative solutions (valid coin counts):")
    valid_solutions = []
    
    for t in range(-10, 11):
        x = x_part + b_coeff * t
        y = y_part + a_coeff * t
        if x >= 0 and y >= 0:
            valid_solutions.append((x, y))
            print(f"  {x} quarters + {y} dimes = ${float((25*x + 10*y)/100):.2f}")
    
    print(f"\nTotal ways: {len(valid_solutions)}")

# Pythagorean Triples
print("PYTHAGOREAN TRIPLES\n")
print("Finding integer solutions to x² + y² = z²")
print("These represent right triangles with integer side lengths!")
print()

def generate_primitive_pythagorean_triples(limit):
    """
    Generate primitive Pythagorean triples using the parametric method
    Formula: a = m² - n², b = 2mn, c = m² + n²
    where m > n > 0, gcd(m,n) = 1, and m,n not both odd
    """
    print("Parametric Generation Method:")
    print("For m > n > 0 with gcd(m,n) = 1 and m,n not both odd:")
    print("  a = m² - n²")
    print("  b = 2mn")
    print("  c = m² + n²")
    print()
    
    triples = []
    
    print(f"{'m':>2} {'n':>2} | {'a':>3} {'b':>3} {'c':>3} | {'a²+b²':>6} {'c²':>6} | {'Check':>5}")
    print("-" * 50)
    
    m = 2
    while m * m + 1 <= limit:  # Ensure c ≤ limit
        for n in range(1, m):
            # Check conditions for primitive triple
            if gcd(m, n) == 1 and (m % 2 != n % 2):  # Not both odd
                a = m*m - n*n
                b = 2*m*n
                c = m*m + n*n
                
                if c <= limit:
                    # Ensure a ≤ b by convention
                    if a > b:
                        a, b = b, a
                    
                    # Verify it's a Pythagorean triple
                    check = a*a + b*b == c*c
                    check_symbol = "✓" if check else "✗"
                    
                    print(f"{m:2d} {n:2d} | {a:3d} {b:3d} {c:3d} | {a*a + b*b:6d} {c*c:6d} | {check_symbol:>5}")
                    
                    if check:
                        triples.append((a, b, c, m, n))
        m += 1
    
    return triples

# Generate primitive triples
limit = 50
primitive_triples = generate_primitive_pythagorean_triples(limit)

print(f"\nFound {len(primitive_triples)} primitive Pythagorean triples with c ≤ {limit}")
print()

# Famous examples
print("FAMOUS PYTHAGOREAN TRIPLES:\n")

famous_triples = [(3, 4, 5), (5, 12, 13), (8, 15, 17), (7, 24, 25), (20, 21, 29)]

print(f"{'Triple':>12} | {'a²':>4} {'b²':>4} {'c²':>4} | {'Sum':>6} | {'Area':>5}")
print("-" * 45)

for a, b, c in famous_triples:
    area = (a * b) // 2  # Right triangle area
    sum_squares = a*a + b*b
    c_squared = c*c
    status = "✓" if sum_squares == c_squared else "✗"
    
    print(f"({a:2d},{b:2d},{c:2d}) | {a*a:4d} {b*b:4d} {c*c:4d} | {sum_squares:6d} | {area:5d}")

print()

# All Pythagorean triples from primitive ones
print("ALL PYTHAGOREAN TRIPLES:\n")
print("Every Pythagorean triple is either primitive or a multiple of a primitive triple.")

# Show multiples of the first primitive triple
if primitive_triples:
    a, b, c, m, n = primitive_triples[0]
    print(f"\nMultiples of primitive triple ({a}, {b}, {c}):")
    
    for k in range(1, 5):
        ka, kb, kc = k*a, k*b, k*c
        print(f"  k={k}: ({ka:2d}, {kb:2d}, {kc:2d})")

print("\n🏛️  Historical Note: Pythagorean triples were known to Babylonians around 1800 BCE!")
print("   Tablet Plimpton 322 lists many examples, predating Pythagoras by over 1000 years.")

# Connection to rational points on unit circle
print("\nGEOMETRIC INTERPRETATION:")
print("Pythagorean triples (a,b,c) correspond to rational points on the unit circle:")
print("  (x,y) = (a/c, b/c) satisfies x² + y² = 1")

if primitive_triples:
    a, b, c, m, n = primitive_triples[0]
    x_coord = Rational(a, c)
    y_coord = Rational(b, c)
    print(f"\nExample: ({a}, {b}, {c}) → ({x_coord}, {y_coord}) on unit circle")
    print(f"Check: ({x_coord})² + ({y_coord})² = {x_coord**2 + y_coord**2} = 1 ✓")

# Practice Problems
print("🎯 PRACTICE PROBLEMS\n")

# Problem 1: Prime factorization and divisors
print("Problem 1: Prime Factorization Analysis")
print("Find the prime factorization of 2520 and count its divisors.")
print()

n = 2520
factorization = factor(n)
print(f"Solution:")
print(f"  2520 = {factorization}")

# Count divisors using formula: d(n) = ∏(αᵢ + 1) where n = ∏pᵢ^αᵢ
divisor_count = 1
for prime, exp in factorization:   # CHANGED
    divisor_count *= (exp + 1)

actual_divisors = divisors(n)
print(f"  Number of divisors: d(2520) = {divisor_count}")
print(f"  Verification: {len(actual_divisors)} divisors found ✓")
print(f"  First 10 divisors: {actual_divisors[:10]}")
print()

# Problem 2: Extended Euclidean Algorithm
print("Problem 2: Bézout's Identity")
print("Find integers x, y such that 252x + 198y = gcd(252, 198)")
print()

a, b = 252, 198
g, x, y = xgcd(a, b)
print(f"Solution:")
print(f"  gcd(252, 198) = {g}")
print(f"  252 × {x} + 198 × {y} = {g}")
print(f"  Verification: 252 × {x} + 198 × {y} = {252*x + 198*y} ✓")
print()

# Problem 3: Chinese Remainder Theorem
print("Problem 3: System of Congruences")
print("Solve the system:")
print("  x ≡ 2 (mod 5)")
print("  x ≡ 3 (mod 7)")
print("  x ≡ 4 (mod 11)")
print()

remainders = [2, 3, 4]
moduli = [5, 7, 11]
solution = crt(remainders, moduli)
M = 5 * 7 * 11

print(f"Solution:")
print(f"  x ≡ {solution} (mod {M})")
print(f"  Verification:")
for r, m in zip(remainders, moduli):
    check = solution % m
    status = "✓" if check == r else "✗"
    print(f"    {solution} ≡ {check} ≡ {r} (mod {m}) {status}")
print()

# Problem 4: Euler's totient function
print("Problem 4: Totient Function")
print("Compute φ(720) using the prime factorization method.")
print()

n = 720
factorization = factor(n)
phi_n = euler_phi(n)

print(f"Solution:")
print(f"  720 = {factorization}")

# Apply φ(n) = n × ∏(1 - 1/p)
result = n
calculation_steps = [f"φ(720) = 720"]

for prime, exp in factorization:   # CHANGED
    result = result * (prime - 1) // prime
    calculation_steps.append(f" × (1 - 1/{prime})")

formula = "".join(calculation_steps)
print(f"  {formula} = {phi_n}")
print(f"  This means {phi_n} integers from 1 to 720 are coprime to 720.")
print()

# Problem 5: Linear Diophantine equation
print("Problem 5: Diophantine Equation")
print("Find all integer solutions to 15x + 25y = 5")
print()

a, b, c = 15, 25, 5
g = gcd(a, b)
print(f"Solution:")
print(f"  gcd(15, 25) = {g}")
print(f"  Since {g} divides {c}, solutions exist.")

# Find particular solution
gcd_result, x0, y0 = xgcd(a, b)
scale = c // g
x_part = x0 * scale
y_part = y0 * scale

print(f"  Particular solution: x = {x_part}, y = {y_part}")
print(f"  General solution: x = {x_part} + {b//g}t, y = {y_part} - {a//g}t")
print(f"  Check: 15×{x_part} + 25×{y_part} = {15*x_part + 25*y_part} = 5 ✓")

print("\nChallenge Problems:")
print("1. Find the next primitive Pythagorean triple after (20, 21, 29)")
print("2. Prove that φ(n) is even for all n > 2")
print("3. Show that if gcd(a,n) = 1, then a^φ(n) ≡ 1 (mod n) (Euler's theorem)")
print("4. Find the smallest positive integer that leaves remainder 1 when divided by")
print("   2, 3, 4, 5, and 6")
print("5. Characterize all integers n such that φ(n) = n/2")