from lib_prime import is_prime_rabin_miller
from functools import cache

def euler_387():
    harshads = []
    for n in range(10, 100):
        ds = n // 10 + n % 10
        if n % ds == 0:
            harshads.append((n, ds))
    limit = 10**14

    r = 0
    for _ in range(13):
        new_harshads = []
        for h, ds in harshads:
            is_strong = is_prime_rabin_miller(h // ds)
            for d in range(10):
                harshad_candidate = h * 10 + d
                if is_strong and d in (1, 3, 7, 9) and harshad_candidate < limit:
                    if is_prime_rabin_miller(harshad_candidate):
                        r += harshad_candidate
                nds = ds + d
                if harshad_candidate % nds != 0:
                    continue
                new_harshads.append((harshad_candidate, nds))
        harshads = new_harshads
    return r

@cache
def digit_sum(n: int) -> int:
    return sum(map(int, str(n)))

@cache
def is_harshad(n: int) -> bool:
    return n % digit_sum(n) == 0

@cache
def trunc_right(n: int) -> int:
    return int(str(n)[:-1])

@cache
def is_strong_harshad(n: int) -> bool:
    ds = digit_sum(n)
    if n % ds != 0:
        return False
    return is_prime_rabin_miller(n // ds)

@cache
def is_right_trunk_harshad(n: int) -> bool:
    if n == 0:
        return False
    if n < 10:
        return True
    if not is_harshad(n):
        return False
    return is_right_trunk_harshad(trunc_right(n))

@cache
def is_strong_right_harshad_prime(n: int) -> bool:
    if not is_prime_rabin_miller(n):
        return False
    nt = trunc_right(n)
    if not is_right_trunk_harshad(nt):
        return False
    if not is_strong_harshad(nt):
        return False
    return True

def asserts():
    assert is_harshad(201)
    assert not is_harshad(122)
    assert trunc_right(201) == 20
    assert trunc_right(10) == 1
    assert is_strong_harshad(201)
    assert is_strong_right_harshad_prime(2011)
    assert is_right_trunk_harshad(201)
    assert not is_right_trunk_harshad(1220)

if __name__ == "__main__":
    solution = euler_387()
    asserts()
    print("e387.py: " + str(solution))
    assert(solution == 696067597313468)