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