"""
Euclid's GCD and Extended Euclidean Algorithm — MAT5EJ302(1) Module I
Includes modular inverse computation.
"""


def gcd(a, b):
    """
    Basic Euclidean algorithm.
    gcd(a, b) = gcd(b, a mod b)
    Time: O(log min(a, b)) steps.
    """
    while b != 0:
        a, b = b, a % b
    return a


def gcd_verbose(a, b):
    """Same as gcd() but prints each step."""
    print(f"  gcd({a}, {b})")
    step = 1
    while b != 0:
        q, r = divmod(a, b)
        print(f"  Step {step}: gcd({a}, {b}) = gcd({b}, {a} mod {b}) = gcd({b}, {r})")
        a, b = b, r
        step += 1
    print(f"  Result: {a}")
    return a


def extended_gcd(a, b):
    """
    Extended Euclidean algorithm.
    Returns (d, x, y) such that a*x + b*y = d = gcd(a, b).
    """
    if b == 0:
        return a, 1, 0
    d, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return d, x, y


def modular_inverse(a, N):
    """
    Returns x such that a*x ≡ 1 (mod N).
    Requires gcd(a, N) = 1.
    """
    d, x, _ = extended_gcd(a, N)
    if d != 1:
        raise ValueError(f"gcd({a}, {N}) = {d} ≠ 1 — inverse does not exist")
    return x % N


def mod_exp(base, exp, mod):
    """Fast modular exponentiation: base^exp mod mod. Time: O(log exp)."""
    result = 1
    base %= mod
    while exp > 0:
        if exp % 2 == 1:
            result = result * base % mod
        exp //= 2
        base = base * base % mod
    return result


if __name__ == "__main__":
    print("=== Basic GCD ===")
    pairs = [(1071, 462), (840, 252), (17, 13), (100, 75), (48, 18)]
    for a, b in pairs:
        print(f"  gcd({a}, {b}) = {gcd(a, b)}")

    print()
    print("=== GCD Step-by-Step ===")
    gcd_verbose(1071, 462)

    print()
    print("=== Extended GCD ===")
    test_pairs = [(35, 15), (17, 13), (100, 37)]
    for a, b in test_pairs:
        d, x, y = extended_gcd(a, b)
        print(f"  {a}*({x}) + {b}*({y}) = {d}  [verify: {a*x + b*y}]")

    print()
    print("=== Modular Inverse ===")
    test_cases = [(3, 7), (5, 11), (7, 26)]
    for a, N in test_cases:
        inv = modular_inverse(a, N)
        print(f"  Inverse of {a} mod {N} = {inv}  [check: {a}*{inv} mod {N} = {(a*inv) % N}]")

    print()
    print("=== Modular Exponentiation ===")
    print(f"  3^100 mod 17 = {mod_exp(3, 100, 17)}")
    print(f"  7^10  mod 11 = {mod_exp(7, 10, 11)}  (should be 1 — Fermat's Little Theorem)")