Modular Inverse Bug — Fermat's Fails on Composite Modulus

# io.thecodeforge: Modular Arithmetic Properties — Complete Reference
# Every property verified with concrete examples.

MOD = 10**9 + 7

# ──────────────────────────────────────────────────────────────────────
# ADDITION: (a + b) mod m = ((a mod m) + (b mod m)) mod m
# ──────────────────────────────────────────────────────────────────────

a, b, m = 10**18, 10**18, 10**9 + 7

result_naive = (a + b) % m
result_distributed = ((a % m) + (b % m)) % m

print('=== Addition ===')
print(f'({a} + {b}) mod {m} = {result_naive}')
print(f'(({a} mod {m}) + ({b} mod {m})) mod {m} = {result_distributed}')
print(f'Match: {result_naive == result_distributed}')
print()

# ──────────────────────────────────────────────────────────────────────
# SUBTRACTION: (a - b) mod m = ((a mod m) - (b mod m) + m) mod m
# The '+m' prevents negative results.
# ──────────────────────────────────────────────────────────────────────

a, b, m = 5, 10, 7

result = (a - b) % m
result_safe = ((a % m) - (b % m) + m) % m

print('=== Subtraction ===')
print(f'({a} - {b}) mod {m} = {result} (Python handles negatives correctly)')
print(f'(({a} mod {m}) - ({b} mod {m}) + {m}) mod {m} = {result_safe}')
print(f'Match: {result == result_safe}')
print()

# ──────────────────────────────────────────────────────────────────────
# MULTIPLICATION: (a × b) mod m = ((a mod m) × (b mod m)) mod m
# This is the property that makes RSA possible.
# ──────────────────────────────────────────────────────────────────────

a, b, m = 10**18, 10**18, 10**9 + 7

result_naive = (a * b) % m
result_distributed = ((a % m) * (b % m)) % m

print('=== Multiplication ===')
print(f'({a} × {b}) mod {m} = {result_naive}')
print(f'(({a} mod {m}) × ({b} mod {m})) mod {m} = {result_distributed}')
print(f'Match: {result_naive == result_distributed}')
print()

# ──────────────────────────────────────────────────────────────────────
# DIVISION: (a / b) mod m = (a × b⁻¹) mod m
# Division does NOT distribute over mod. Use modular inverse.
# ──────────────────────────────────────────────────────────────────────

a, b, m = 10, 3, 10**9 + 7

# Can't just do (a / b) mod m — division doesn't work in mod
# Instead: find b⁻¹ such that b × b⁻¹ ≡ 1 (mod m)
b_inv = pow(b, m - 2, m)  # Fermat's little theorem (m must be prime)
result = a * b_inv % m

# Verify: result × b mod m should equal a
verify = result * b % m

print('=== Division ===')
print(f'{a} / {b} mod {m} = {result}')
print(f'Verify: {result} × {b} mod {m} = {verify} (should be {a})')
print(f'Match: {verify == a}')
print()

# ──────────────────────────────────────────────────────────────────────
# NEGATIVE MODULO: Python vs C++ behavior
# This is a common interview trap.
# ──────────────────────────────────────────────────────────────────────

print('=== Negative Modulo ===')
print(f'Python: (-3) % 7 = {(-3) % 7} (always non-negative)')
print(f'Python: (-10) % 3 = {(-10) % 3}')
print()
print('C/C++: (-3) % 7 = -3 (keeps the sign)')
print('C/C++: (-10) % 3 = -1')
print()
print('Fix for C/C++: ((a % m) + m) % m')
print(f'Python equivalent: ((-3) % 7 + 7) % 7 = {((-3) % 7 + 7) % 7}')
print()

# ──────────────────────────────────────────────────────────────────────
# PROPERTIES SUMMARY TABLE
# ──────────────────────────────────────────────────────────────────────

print('=== Properties Summary ===')
properties = [
    ('(a + b) mod m', 'Distributes', 'Yes'),
    ('(a - b) mod m', 'Distributes (add m to fix negatives)', 'Yes'),
    ('(a × b) mod m', 'Distributes', 'Yes'),
    ('(a^b) mod m', 'Use fast modexp', 'Yes'),
    ('(a / b) mod m', 'Does NOT distribute — use inverse', 'No'),
    ('a mod (-m)', 'Undefined convention — avoid negative moduli', 'N/A'),
]

for expr, note, distributes in properties:
    print(f'{expr:20s} | Distributes: {distributes:3s} | {note}')

# io.thecodeforge: Fast Modular Exponentiation
# O(log exp) using repeated squaring. The most important algorithm in this guide.

import time

def mod_exp_manual(base: int, exp: int, mod: int) -> int:
    """
    Compute base^exp mod mod using repeated squaring.
    O(log exp) multiplications.
    
    How it works:
    - Write exp in binary
    - For each bit (right to left):
      - If bit is 1: result = result * base mod m
      - Square: base = base * base mod m
    - Return result
    """
    result = 1
    base %= mod  # reduce base first (handles base > mod)
    while exp > 0:
        if exp & 1:  # if current bit is 1
            result = result * base % mod
        base = base * base % mod  # square the base
        exp >>= 1  # shift to next bit
    return result

# ──────────────────────────────────────────────────────────────────────
# VERIFICATION: manual vs Python built-in
# ──────────────────────────────────────────────────────────────────────

test_cases = [
    (2, 10, 1000),          # 2^10 mod 1000 = 1024 mod 1000 = 24
    (3, 200, 13),           # 3^200 mod 13
    (7, 1000000, 10**9+7),  # large exponent
    (2, 10**18, 10**9+7),   # massive exponent
]

print('=== Manual vs Built-in ===')
for base, exp, mod in test_cases:
    manual = mod_exp_manual(base, exp, mod)
    builtin = pow(base, exp, mod)
    match = manual == builtin
    exp_str = f'{exp}' if exp < 10000 else f'10^{len(str(exp))-1}'
    print(f'{base}^{exp_str} mod {mod} = {manual} | Match: {match}')
print()

# ──────────────────────────────────────────────────────────────────────
# PERFORMANCE: built-in pow() is ~100x faster (C implementation)
# ──────────────────────────────────────────────────────────────────────

base, exp, mod = 7, 10**18, 10**9 + 7
iterations = 10000

start = time.perf_counter()
for _ in range(iterations):
    mod_exp_manual(base, exp, mod)
manual_time = (time.perf_counter() - start) / iterations

start = time.perf_counter()
for _ in range(iterations):
    pow(base, exp, mod)
builtin_time = (time.perf_counter() - start) / iterations

print('=== Performance ===')
print(f'Manual mod_exp: {manual_time*1e6:.2f} µs per call')
print(f'Built-in pow(): {builtin_time*1e6:.2f} µs per call')
print(f'Speedup: {manual_time/builtin_time:.1f}x')
print('Always use pow(base, exp, mod) in production and competitions.')
print()

# ──────────────────────────────────────────────────────────────────────
# STEP-BY-STEP TRACE: understand what happens inside
# ──────────────────────────────────────────────────────────────────────

print('=== Step-by-Step Trace: 3^13 mod 100 ===')
base, exp, mod = 3, 13, 100
result = 1
base = base % mod
step = 0

while exp > 0:
    step += 1
    bit = exp & 1
    print(f'Step {step}: exp={exp:>4} ({bin(exp):>8}) bit={bit} | result={result:>3} base={base:>3}', end='')
    if bit:
        result = result * base % mod
        print(f' → result={result} (multiply)')
    else:
        print()
    base = base * base % mod
    exp >>= 1

print(f'Final answer: {result}')

// io.thecodeforge: Fast Modular Exponentiation in C++
// Uses unsigned long long and __int128 for safe multiplication

#include <cstdint>
#include <iostream>
using namespace std;

typedef unsigned long long ull;
typedef __int128 int128;

// Modular exponentiation: (base^exp) % mod
ull mod_pow(ull base, ull exp, ull mod) {
    ull result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) {
            result = (int128)result * base % mod;
        }
        base = (int128)base * base % mod;
        exp >>= 1;
    }
    return result;
}

// If __int128 is not available (e.g., MSVC), use more portable approach:
// ull mul_mod(ull a, ull b, ull m) {
//     ull res = 0;
//     a %= m;
//     while (b) {
//         if (b & 1) res = (res + a) % m;
//         a = (a << 1) % m;
//         b >>= 1;
//     }
//     return res;
// }

int main() {
    // Test cases
    cout << "3^13 mod 100 = " << mod_pow(3, 13, 100) << endl;
    cout << "2^10 mod 1000 = " << mod_pow(2, 10, 1000) << endl;
    cout << "7^(10^6) mod (10^9+7) = " << mod_pow(7, 1000000, 1000000007ULL) << endl;
    cout << "3^1000000005 mod (10^9+7) = " << mod_pow(3, 1000000005ULL, 1000000007ULL) << endl; // Fermat inverse
    return 0;
}

# io.thecodeforge: Modular Inverse — Three Methods Compared

from math import gcd

# ──────────────────────────────────────────────────────────────────────
# METHOD 1: FERMAT'S LITTLE THEOREM
# a^(-1) ≡ a^(p-2) mod p  (only when p is prime)
# O(log p) — one modular exponentiation
# ──────────────────────────────────────────────────────────────────────

def mod_inv_fermat(a: int, p: int) -> int:
    """Modular inverse using Fermat's little theorem. p must be prime."""
    return pow(a, p - 2, p)

# ──────────────────────────────────────────────────────────────────────
# METHOD 2: EXTENDED EUCLIDEAN ALGORITHM
# Finds x such that a*x + m*y = gcd(a, m)
# If gcd(a, m) = 1, then a*x ≡ 1 (mod m), so x = a^(-1) mod m
# Works for ANY m (not just primes). O(log m).
# ──────────────────────────────────────────────────────────────────────

def extended_gcd(a: int, b: int) -> tuple:
    """Returns (gcd, x, y) such that a*x + b*y = gcd(a, b)."""
    if b == 0:
        return a, 1, 0
    g, x, y = extended_gcd(b, a % b)
    return g, y, x - (a // b) * y

def mod_inv_extended_gcd(a: int, m: int) -> int:
    """Modular inverse using Extended GCD. Works for any m where gcd(a,m)=1."""
    g, x, y = extended_gcd(a, m)
    if g != 1:
        raise ValueError(f"Modular inverse does not exist: gcd({a}, {m}) = {g}")
    return x % m

# ──────────────────────────────────────────────────────────────────────
# METHOD 3: EULER'S THEOREM
# a^phi(m) ≡ 1 mod m, so a^(-1) ≡ a^(phi(m)-1) mod m
# Requires computing phi(m), which factors m. Educational only.
# ──────────────────────────────────────────────────────────────────────

def euler_totient_simple(n: int) -> int:
    """Compute Euler's totient phi(n) by trial division."""
    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

def mod_inv_euler(a: int, m: int) -> int:
    """Modular inverse using Euler's theorem. Works for any m where gcd(a,m)=1."""
    phi = euler_totient_simple(m)
    return pow(a, phi - 1, m)

# ──────────────────────────────────────────────────────────────────────
# COMPARISON: all three methods give the same result
# ──────────────────────────────────────────────────────────────────────

print('=== Comparison (m = 101, prime) ===')
m = 101  # prime
for a in [3, 7, 42, 99]:
    fermat = mod_inv_fermat(a, m)
    ext_gcd = mod_inv_extended_gcd(a, m)
    euler = mod_inv_euler(a, m)
    print(f'a={a:>3}: Fermat={fermat:>4} | ExtGCD={ext_gcd:>4} | Euler={euler:>4} | Match: {fermat == ext_gcd == euler}')
print()

print('=== Extended GCD works for non-prime moduli ===')
m = 100  # NOT prime
for a in [3, 7, 11, 99]:
    try:
        inv = mod_inv_extended_gcd(a, m)
        verify = a * inv % m
        print(f'a={a:>3}, m={m}: inverse={inv:>4} | verify: {a}×{inv} mod {m} = {verify}')
    except ValueError as e:
        print(f'a={a:>3}, m={m}: {e}')

print()
print('Note: a=99 has no inverse mod 100 because gcd(99, 100) = 1... wait, it does!')
print('But a=2 has no inverse mod 100 because gcd(2, 100) = 2 ≠ 1')
inv_2 = None
try:
    inv_2 = mod_inv_extended_gcd(2, 100)
except ValueError as e:
    print(f'a=2, m=100: {e}')
print()

# ──────────────────────────────────────────────────────────────────────
# PYTHON 3.8+ BUILT-IN: pow(a, -1, m)
# ──────────────────────────────────────────────────────────────────────

print('=== Python 3.8+ Built-in ===')
a, m = 3, 10**9 + 7
builtin_inv = pow(a, -1, m)
fermat_inv = pow(a, m - 2, m)
print(f'pow({a}, -1, {m}) = {builtin_inv}')
print(f'pow({a}, {m}-2, {m}) = {fermat_inv}')
print(f'Match: {builtin_inv == fermat_inv}')
print(f'Verify: {a} × {builtin_inv} mod {m} = {a * builtin_inv % m}')

# io.thecodeforge: Comparison of Fermat vs Extended GCD

MOD_PRIME = 998244353  # common CP prime
MOD_COMPOSITE = 1000000  # not prime

def modinv_fermat(a, mod):
    return pow(a, mod - 2, mod)

def modinv_extgcd(a, mod):
    # Extended GCD: solves a*x + mod*y = 1
    def egcd(a, b):
        if b == 0:
            return a, 1, 0
        g, x, y = egcd(b, a % b)
        return g, y, x - (a // b) * y
    g, x, y = egcd(a, mod)
    if g != 1:
        return None
    return x % mod

print("=== On prime modulus ===")
a = 3
print(f"Fermat: {modinv_fermat(a, MOD_PRIME)}")
print(f"ExtGCD: {modinv_extgcd(a, MOD_PRIME)}")

print("\n=== On composite modulus ===")
print(f"Fermat: {modinv_fermat(a, MOD_COMPOSITE)}")
print(f"ExtGCD: {modinv_extgcd(a, MOD_COMPOSITE)}")
# Verify
inv = modinv_extgcd(a, MOD_COMPOSITE)
print(f"Verify: {a} * {inv} % {MOD_COMPOSITE} = {a * inv % MOD_COMPOSITE} (should be 1)")
# Fermat gives wrong result
fermat_inv = modinv_fermat(a, MOD_COMPOSITE)
print(f"Fermat's inverse computed: {fermat_inv}")
print(f"Verify Fermat result: {a} * {fermat_inv} % {MOD_COMPOSITE} = {a * fermat_inv % MOD_COMPOSITE} (should be 1)")

# io.thecodeforge: Modular Division — The Pattern That Trips Everyone Up

MOD = 10**9 + 7

def mod_inv(a, mod=MOD):
    return pow(a, mod - 2, mod)

# ──────────────────────────────────────────────────────────────────────
# WRONG WAY: integer division after mod — gives wrong answer
# ──────────────────────────────────────────────────────────────────────

a, b = 10, 3

wrong = (a % MOD) // b  # integer division — WRONG in modular arithmetic
right = (a % MOD) * mod_inv(b) % MOD  # modular inverse — CORRECT

print('=== Modular Division ===')
print(f'{a} / {b} mod {MOD} = ')
print(f'  Integer division (WRONG): {wrong}')
print(f'  Modular inverse (CORRECT): {right}')
print(f'  Verify: {right} * {b} % {MOD} = {right * b % MOD} (should be {a})')
print()

# Large number example
a_large = 123456789012345678
print(f'{a_large} / {b} mod {MOD} = {a_large * mod_inv(b) % MOD}')
print(f'Verify: {a_large % MOD} = {a_large * mod_inv(b) % MOD * b % MOD}')
print()

print('=== Floor Division (Not Defined) ===')
print('In modular arithmetic, division yields the unique result in [0, m-1] such that')
print('a = (a/b) * b mod m. There is no floor or truncation.')
print('If you need floor division, you must work in integers and only reduce at the end.')

# io.thecodeforge: Fermat's Little Theorem — Applications and Proof Verification

import random

MOD = 13  # prime

# Verify FLT for random a
def check_flt(a, p):
    return pow(a, p-1, p) == 1

print('=== Fermat Verification (mod 13) ===')
for a in range(1, MOD):
    print(f'a={a:2d}: a^{MOD-1} mod {MOD} = {pow(a, MOD-1, MOD)}')
print()

# Exponent reduction: a^b mod p = a^(b mod (p-1)) mod p
print('=== Exponent Reduction ===')
a = 5
b_large = 10**18
b_reduced = b_large % (MOD - 1)
full = pow(a, b_large, MOD)
reduced = pow(a, b_reduced, MOD)
print(f'5^(10^18) mod {MOD} = {full}')
print(f'5^(10^18 mod {MOD-1}) mod {MOD} = {reduced}')
print(f'Match: {full == reduced}')
print()

# Caution: FLT gives wrong result on composite modulus
print('=== FLT on Composite Modulus ===')
a = 2
n = 15  # not prime
print(f'2^{n-1} mod {n} = {pow(a, n-1, n)} (should be 1 if n prime, but is {pow(a, n-1, n)})')
print('This is why Fermat's test can have false positives for Carmichael numbers.')

# io.thecodeforge: Chinese Remainder Theorem — Complete Implementation

from math import prod

def extended_gcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x, y = extended_gcd(b, a % b)
    return g, x, y  # returns (g, x, y) with a*x + b*y = g

def crt(remainders, moduli):
    """
    Solve the system of congruences:
    x ≡ remainders[i] (mod moduli[i])
    Assumes moduli are pairwise coprime.
    Returns (x, M) where x is the unique solution modulo M.
    """
    M = prod(moduli)
    result = 0
    for a, m in zip(remainders, moduli):
        Mi = M // m
        g, x, y = extended_gcd(Mi, m)  # Mi*x + m*y = 1
        inv = x % m
        result += a * Mi * inv
    return result % M, M

# Example
remainders = [2, 3, 2]
moduli = [3, 5, 7]
x, M = crt(remainders, moduli)
print(f'Solution: x ≡ {x} (mod {M})')
for a, m in zip(remainders, moduli):
    print(f'  Check: {x} mod {m} = {x % m} (expected {a})')
print()

# Large example: x mod 11 = 10, mod 13 = 12, mod 17 = 16
remainders2 = [10, 12, 16]
moduli2 = [11, 13, 17]
x2, M2 = crt(remainders2, moduli2)
print(f'Solution: x ≡ {x2} (mod {M2})')
for a, m in zip(remainders2, moduli2):
    print(f'  Check: {x2} mod {m} = {x2 % m} (expected {a})')

# io.thecodeforge: Fast nCr Mod p — Precomputation and Queries

MOD = 10**9 + 7
MAX_N = 10**6  # adjust as needed

# Precompute factorials and inverse factorials
fact = [1] * (MAX_N + 1)
inv_fact = [1] * (MAX_N + 1)

for i in range(1, MAX_N + 1):
    fact[i] = fact[i-1] * i % MOD

inv_fact[MAX_N] = pow(fact[MAX_N], MOD - 2, MOD)
for i in range(MAX_N, 0, -1):
    inv_fact[i-1] = inv_fact[i] * i % MOD

def nCr(n: int, r: int) -> int:
    if r < 0 or r > n:
        return 0
    return fact[n] * inv_fact[r] % MOD * inv_fact[n-r] % MOD

# Test queries
queries = [(5,2), (10,3), (100,50), (0,0), (10,10)]
for n,r in queries:
    print(f'C({n},{r}) mod {MOD} = {nCr(n,r)}')
print()

# Verify small values manually
from math import comb
print('=== Verification with Python int comb ===')
for n,r in [(5,2), (10,3), (0,0)]:
    expected = comb(n,r) % MOD
    result = nCr(n,r)
    print(f'C({n},{r}) = {result}, expected {expected}, match: {result == expected}')

// io.thecodeforge — dsa tutorial

public class ModInt {
    private static final long MOD = 1_000_000_007L;
    private final long value;

    private ModInt(long value) {
        this.value = ((value % MOD) + MOD) % MOD;
    }

    public static ModInt of(long raw) {
        return new ModInt(raw);
    }

    public ModInt add(ModInt other) {
        return new ModInt(this.value + other.value);
    }

    public ModInt multiply(ModInt other) {
        // Java long multiplication won't overflow for MOD up to 1e9+7
        return new ModInt((this.value * other.value) % MOD);
    }

    public long toLong() {
        return value;
    }

    public static void main(String[] args) {
        ModInt a = ModInt.of(500000004);
        ModInt b = ModInt.of(3);
        // 500000004 * 3 = 1500000012 mod 1e9+7 = 500000005
        System.out.println(a.multiply(b).toLong());
    }
}

// io.thecodeforge — dsa tutorial

public class RollingHashSearch {
    private static final long BASE = 131L;
    private static final long MOD = 1_000_000_007L;

    public static int findPattern(String text, String pattern) {
        int n = text.length(), m = pattern.length();
        if (m > n) return -1;

        long[] pow = new long[n + 1];
        long[] hash = new long[n + 1];
        pow[0] = 1;
        for (int i = 1; i <= n; i++) {
            pow[i] = (pow[i-1] * BASE) % MOD;
            hash[i] = (hash[i-1] * BASE + (text.charAt(i-1) - 'a' + 1)) % MOD;
        }

        long patternHash = 0;
        for (int i = 0; i < m; i++)
            patternHash = (patternHash * BASE + (pattern.charAt(i) - 'a' + 1)) % MOD;

        for (int i = 0; i <= n - m; i++) {
            long subHash = (hash[i + m] - hash[i] * pow[m] % MOD + MOD) % MOD;
            if (subHash == patternHash) return i;
        }
        return -1;
    }

    public static void main(String[] args) {
        System.out.println(findPattern("competitiveprogramming", "prog"));
    }
}

// io.thecodeforge — dsa tutorial

public class SafeOperations {
    static final long MOD = 1_000_000_007L;

    static long safeAdd(long a, long b) {
        return ((a % MOD) + (b % MOD)) % MOD;
    }

    static long safeMul(long a, long b) {
        return ((a % MOD) * (b % MOD)) % MOD;
    }

    // Fermat inverse: b^(MOD-2) mod MOD (MOD must be prime)
    static long modInverse(long b) {
        return fastPow(b, MOD - 2);
    }

    static long fastPow(long base, long exp) {
        long result = 1;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1) == 1) result = (result * base) % MOD;
            base = (base * base) % MOD;
            exp >>= 1;
        }
        return result;
    }

    public static void main(String[] args) {
        long a = 500000004, b = 3;
        System.out.println(safeMul(a, b));           // 500000005
        System.out.println(modInverse(2));           // 500000004
        System.out.println(safeMul(10, modInverse(3))); // 333333336 (10/3 mod 1e9+7)
    }
}