Number Theory

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def extended_gcd(a, b):
    if a == 0:
        return b, 0, 1
    g, x, y = extended_gcd(b % a, a)
    return g, y - (b // a) * x, x

def mod_inverse(a, m):
    g, x, _ = extended_gcd(a % m, m)
    if g != 1:
        raise ValueError("Modular inverse does not exist")
    return x % m

def lcm(a, b):
    return a * b // gcd(a, b)

print(f"gcd(48, 18) = {gcd(48, 18)}")
print(f"lcm(4, 6) = {lcm(4, 6)}")
print(f"mod_inverse(3, 7) = {mod_inverse(3, 7)}")

def mod_pow(base, exp, mod):
    result = 1
    base %= mod
    while exp > 0:
        if exp % 2 == 1:
            result = (result * base) % mod
        exp //= 2
        base = (base * base) % mod
    return result

def is_prime_miller_rabin(n, k=10):
    if n < 2:
        return False
    if n == 2 or n == 3:
        return True
    if n % 2 == 0:
        return False

    d, r = n - 1, 0
    while d % 2 == 0:
        d //= 2
        r += 1

    for _ in range(k):
        a = 2 + (hash(str(_)) % (n - 4))
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True

print(f"7^222 mod 11 = {mod_pow(7, 222, 11)}")
print(f"Is 97 prime: {is_prime_miller_rabin(97)}")

def sieve_of_eratosthenes(limit):
    is_prime = [True] * (limit + 1)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(limit**0.5) + 1):
        if is_prime[i]:
            for j in range(i*i, limit + 1, i):
                is_prime[j] = False
    return [i for i in range(2, limit + 1) if is_prime[i]]

primes = sieve_of_eratosthenes(100)
print(f"Primes below 100: {primes}")
print(f"Count: {len(primes)}")

def euler_totient(n):
    result = n
    p = 2
    temp = n
    while p * p <= temp:
        if temp % p == 0:
            while temp % p == 0:
                temp //= p
            result -= result // p
        p += 1
    if temp > 1:
        result -= result // temp
    return result

for n in range(1, 21):
    print(f"phi({n}) = {euler_totient(n)}")

def chinese_remainder_theorem(remainders, moduli):
    M = 1
    for m in moduli:
        M *= m
    x = 0
    for i in range(len(moduli)):
        Mi = M // moduli[i]
        yi = mod_inverse(Mi, moduli[i])
        x += remainders[i] * Mi * yi
    return x % M

# Solve: x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7)
remainders = [2, 3, 2]
moduli = [3, 5, 7]
print(f"Solution: x = {chinese_remainder_theorem(remainders, moduli)}")

def trial_division(n):
    if n < 2:
        return False
    if n < 4:
        return True
    if n % 2 == 0 or n % 3 == 0:
        return False
    i = 5
    while i * i <= n:
        if n % i == 0 or n % (i + 2) == 0:
            return False
        i += 6
    return True

def fermat_test(n, k=10):
    if n < 2:
        return False
    if n == 2 or n == 3:
        return True
    for _ in range(k):
        a = 2 + (hash(str(_)) % (n - 3))
        if pow(a, n - 1, n) != 1:
            return False
    return True

# Compare methods
test_numbers = [97, 561, 1105, 1729, 2147483647]
for n in test_numbers:
    td = trial_division(n)
    fr = fermat_test(n)
    print(f"{n}: trial={td}, fermat={fr}")

def legendre_symbol(a, p):
    if p < 2 or p % 2 == 0:
        raise ValueError("p must be an odd prime")
    a = a % p
    if a == 0:
        return 0
    ls = pow(a, (p - 1) // 2, p)
    return -1 if ls == p - 1 else ls

def quadratic_residues(p):
    return sorted(set(pow(x, 2, p) for x in range(p)))

for p in [7, 11, 13]:
    qr = quadratic_residues(p)
    print(f"p={p}: QR={qr}, QNR={[x for x in range(1, p) if x not in qr]}")