Make 134 much faster with math(TM).

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
-