Solve problem 381

python/e381.py

@@ -0,0 +1,42 @@
+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

python/lib_misc.py

@@ -230,3 +230,18 @@ def cache(f):
         return r
     return func_cached
 
+
+def ext_gcd(a: int, b: int) -> tuple[int, int, int]:
+    if a == 0:
+        return (b, 0, 1)
+    else:
+        gcd, x, y = ext_gcd(b % a, a)
+        return (gcd, y - (b // a) * x, x)
+
+
+def modinv(a: int, b: int) -> int:
+    gcd, x, _ = ext_gcd(a, b)
+    if gcd != 1:
+        raise Exception('Modular inverse does not exist')
+    else:
+        return x % b