Solve problem 387

python/e387.py

@@ -0,0 +1,84 @@
+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)
+

python/lib_prime.py

@@ -68,7 +68,7 @@ def is_prime(n):
 
 def is_prime_rabin_miller(number):
     """ Rabin-Miller Primality Test """
-    witnesses = [2, 3, 5, 7, 11, 13, 17, 23, 29, 31, 37, 41, 43, 47, 53]
+    witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
 
     if number < 2:
         return False