73 lines
1.5 KiB
Python
73 lines
1.5 KiB
Python
from lib_prime import primes
|
|
from math import log, ceil
|
|
|
|
|
|
def s(p1, p2):
|
|
p1l = ceil(log(p1, 10))
|
|
base = 10**p1l
|
|
for lhs in range(base, 10**12, base):
|
|
r = lhs + p1
|
|
if r % p2 == 0:
|
|
return r
|
|
assert False
|
|
|
|
|
|
def s2(p1, p2):
|
|
# Invert, always invert... but still slow, haha
|
|
p1l = ceil(log(p1, 10))
|
|
base = 10**p1l
|
|
c = p2
|
|
while True:
|
|
if c % base == p1:
|
|
return c
|
|
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 = s3(p1, p2)
|
|
r += sc
|
|
return r
|
|
|
|
|
|
if __name__ == "__main__":
|
|
solution = euler_134()
|
|
print("e134.py: " + str(solution))
|
|
assert solution == 18613426663617118
|