diff --git a/python/e134.py b/python/e134.py index 0f5452f..df6d08e 100644 --- a/python/e134.py +++ b/python/e134.py @@ -13,7 +13,7 @@ def s(p1, p2): def s2(p1, p2): - # Invert, always invert ;) + # Invert, always invert... but still slow, haha p1l = ceil(log(p1, 10)) base = 10**p1l c = p2 @@ -23,12 +23,45 @@ def s2(p1, p2): c += p2 +def ext_gcd(a, b): + 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, b): + gcd, x, _ = ext_gcd(a, b) + if gcd != 1: + raise Exception('Modular inverse does not exist') + else: + return x % b + + +def s3(p1, p2): + # Okay with some math we are fast. + d = 10**ceil(log(p1, 10)) + dinv = modinv(d, p2) + k = ((p2 - p1) * dinv) % p2 + r = k * d + p1 + return r + + def euler_134(): + ps = primes(10000) + for i in range(2, len(ps) - 1): + p1, p2 = ps[i], ps[i + 1] + sc = s(p1, p2) + sc2 = s2(p1, p2) + sc3 = s3(p1, p2) + assert sc == sc2 and sc2 == sc3 + r = 0 ps = primes(10**6 + 10) # p1 < 10**6 but p2 is the first p > 10**6 for i in range(2, len(ps) - 1): p1, p2 = ps[i], ps[i + 1] - sc = s(p1, p2) + sc = s3(p1, p2) r += sc return r @@ -37,4 +70,3 @@ if __name__ == "__main__": solution = euler_134() print("e134.py: " + str(solution)) assert solution == 18613426663617118 -