from lib_prime import primes
from lib_misc import modinv


def s(p: int) -> int:
    a = p - 1
    fp = p - 1
    r = fp
    # Calculate (!(p - 1) + !(p - 2) + ... + !(p - 5)) % p
    for _ in range(4):
        fp = fp * modinv(a, p) % p
        a -= 1
        r += fp
        r %= p
    return r


def euler_381():
    assert s(5) == 4
    assert s(7) == 4

    # Example given by problem statement
    t = 0
    for p in primes(100):
        if p < 5:
            continue
        t += s(p)
    assert t == 480

    # Actual solution (#slow)
    t = 0
    for p in primes(10**8):
        if p < 5:
            continue
        t += s(p)
    return t


if __name__ == "__main__":
    solution = euler_381()
    print(f"e381.py: {solution}")
    assert solution == 139602943319822