From 63a76b01f8ec07e41f1b49a0124541d9277ea502 Mon Sep 17 00:00:00 2001
From: Tibor Bizjak <tibor.bizjak@gmail.com>
Date: Sat, 25 Mar 2023 16:35:17 +0100
Subject: [PATCH] Adding solutions from highest to lowest. Done for N>200.
 Updated to python3.

---
 203_squarefree_binomial.py                    | 38 +++++++++++
 204_generalised_hamming_numbers.py            | 22 ++++++
 205_dice_game.py                              | 17 +++++
 206_concealed_square.py                       | 23 +++++++
 211_divisor_square_sum.py                     | 42 ++++++++++++
 217_balanced_numbers.py                       | 40 +++++++++++
 ..._factorisation_of_binomial_coefficients.py | 42 ++++++++++++
 243_resilience.py                             | 48 +++++++++++++
 315_digital_clock_roots.py                    | 58 ++++++++++++++++
 329_prime_frog.py                             | 42 ++++++++++++
 346_strong_repunits.py                        | 16 +++++
 347_largest_integer.py                        | 29 ++++++++
 357_prime_generating_integers.py              | 54 +++++++++++++++
 381_prime-k_factorial.py                      | 10 +++
 387_harshad_numbers.py                        | 44 ++++++++++++
 429_sum_of_squares.py                         | 21 ++++++
 493_under_the_rainbow.py                      | 20 ++++++
 504_square_on_the_inside.py                   | 48 +++++++++++++
 518_prime_triples.py                          | 38 +++++++++++
 549_divisibility_of_factorials.py             | 30 ++++++++
 data/digs.txt                                 | 30 ++++++++
 lib.py                                        | 68 +++++++++++++++++++
 22 files changed, 780 insertions(+)
 create mode 100644 203_squarefree_binomial.py
 create mode 100644 204_generalised_hamming_numbers.py
 create mode 100644 205_dice_game.py
 create mode 100644 206_concealed_square.py
 create mode 100644 211_divisor_square_sum.py
 create mode 100644 217_balanced_numbers.py
 create mode 100644 231_prime_factorisation_of_binomial_coefficients.py
 create mode 100644 243_resilience.py
 create mode 100644 315_digital_clock_roots.py
 create mode 100644 329_prime_frog.py
 create mode 100644 346_strong_repunits.py
 create mode 100644 347_largest_integer.py
 create mode 100644 357_prime_generating_integers.py
 create mode 100644 381_prime-k_factorial.py
 create mode 100644 387_harshad_numbers.py
 create mode 100644 429_sum_of_squares.py
 create mode 100644 493_under_the_rainbow.py
 create mode 100644 504_square_on_the_inside.py
 create mode 100644 518_prime_triples.py
 create mode 100644 549_divisibility_of_factorials.py
 create mode 100644 data/digs.txt
 create mode 100644 lib.py

diff --git a/203_squarefree_binomial.py b/203_squarefree_binomial.py
new file mode 100644
index 0000000..32bc5e9
--- /dev/null
+++ b/203_squarefree_binomial.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+
+from operator import mul
+from functools import reduce
+
+def factorize(bound):
+    sieve = [[] for i in range(bound)]
+    for n in range(2, bound):
+        if sieve[n] == []:
+            sieve[n] = [n]
+        for i in range(2, n+1):
+            if i*n >= bound:
+                break
+            sieve[i*n] = sieve[i] + sieve[n]
+    return sieve
+
+def bin(row, col, fnums):
+    num = [f for a in fnums[row-col+1:row+1] for f in a]
+    den = [f for a in fnums[1:col+1] for f in a]
+    for f in set(num):
+        if num.count(f) - den.count(f) > 1:
+            return 0
+    return reduce(mul, num) // reduce(mul, den)
+
+
+rn = 51
+fnums = [[1] + i for i in factorize(rn)]
+t = 0
+dup = []
+for row in range(2, rn):
+    for col in range(1, row//2+1):
+        n = bin(row, col, fnums)
+        if n in dup:
+            continue
+        t += n
+        dup.append(n)
+
+print(t+1)
diff --git a/204_generalised_hamming_numbers.py b/204_generalised_hamming_numbers.py
new file mode 100644
index 0000000..f1d0dc4
--- /dev/null
+++ b/204_generalised_hamming_numbers.py
@@ -0,0 +1,22 @@
+#!/usr/bin/env python3
+
+def primes(n):
+    sieve = [False]*n
+    for p in range(2, n):
+        if sieve[p]:
+            continue
+        for i in range(p**2, n, p):
+            sieve[i] = True
+    return [p for p, t in enumerate(sieve) if not t][2:]
+
+def rec(n, nums, i, lim):
+    c = 1
+    for j, f in enumerate(nums[i:], start=i):
+        new = f*n
+        if new > lim:
+            break
+        c += rec(new, nums, j, lim)
+    return c
+
+nums = primes(100)
+print(rec(1, nums, 0, 10**9))
diff --git a/205_dice_game.py b/205_dice_game.py
new file mode 100644
index 0000000..a63e198
--- /dev/null
+++ b/205_dice_game.py
@@ -0,0 +1,17 @@
+#!/usr/bin/env python3
+
+from itertools import product
+
+peter = [0]*(4*9 + 1)
+colin = [0]*(6*6 + 1)
+
+for  p in product([1,2,3,4], repeat=9):
+    peter[sum(p)] += 1
+for p in product([1,2,3,4,5,6], repeat=6):
+    colin[sum(p)] += 1
+
+c = 0.0
+for i, p in enumerate(peter):
+    c += p * sum(colin[:i])
+
+print("{0:.7f}".format(c/(6**6 * 4**9)))
diff --git a/206_concealed_square.py b/206_concealed_square.py
new file mode 100644
index 0000000..d7681de
--- /dev/null
+++ b/206_concealed_square.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+
+# Brute force
+from math import sqrt
+
+def main():
+    low = 10*(int(sqrt(10203040506070809))//10)
+    upp = int(sqrt(19293949596979899))
+    for i in range(low+3, upp+1, 10):
+        if is_num(i**2):
+            return i*10
+    for i in range(low+7, upp+1, 10):
+        if is_num(i**2):
+            return i*10
+    return False
+
+def is_num(n):
+    s = str(n)
+    for i in range(9):
+        if int(s[2*i]) != i+1: return False
+    return True
+
+print(main())
diff --git a/211_divisor_square_sum.py b/211_divisor_square_sum.py
new file mode 100644
index 0000000..ec884e3
--- /dev/null
+++ b/211_divisor_square_sum.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+from math import ceil
+
+def is_square(n):
+    return int(n**0.5)**2 == n
+
+def getpw(n, p):
+    pw = 0
+    while not n % p:
+        pw += 1
+        n //= p
+    return pw
+
+N = 64 * 10**6
+
+sieve = [1] * N
+
+for p in range(2, N, 2):
+    sieve[p] = 2
+
+for p in range(3, int(ceil(N**0.5 - 1)), 2):
+    if sieve[p] > 1:
+        continue
+    sieve[p] = p
+    for i in range(p**2, N, p):
+        sieve[i] = p
+
+print("Sieved")
+s = 1
+
+for n in range(2, N):
+    if sieve[n] == 1:
+        sieve[n] = n
+    p = sieve[n]
+    pw = p ** getpw(n, p)
+    sieve[n] = sieve[n // p] + \
+               pw**2 * sieve[n // pw]
+    if is_square(sieve[n]):
+        s += n
+
+print(s)
diff --git a/217_balanced_numbers.py b/217_balanced_numbers.py
new file mode 100644
index 0000000..07c326c
--- /dev/null
+++ b/217_balanced_numbers.py
@@ -0,0 +1,40 @@
+#!/usr/bin/env python3
+
+def next(sums, counts, mod=0):
+    ln = len(sums) + 9
+    news = [0] * ln
+    newc = [0] * ln
+    for i, (s, c) in enumerate(zip(sums, counts)):
+        for n in range(10):
+            news[i + n] += 10 * s + c * n
+            newc[i + n] += c
+    return news, newc
+
+
+left_sums  = range(10)
+left_count = [0] + [1] * 9
+
+right_sums  = range(10)
+right_count = [1] * 10
+
+
+length = 47
+mod = 3**15
+s = 45
+
+for nofd in range(2, length+1):
+    partial = sum(l * 10**(nofd//2 + nofd%2) * rc + r * lc
+                  for l, r, lc, rc in
+                  zip(left_sums, right_sums, left_count, right_count))
+    if not nofd % 2:
+        s += partial
+        s %= mod
+        continue
+    total = sum(lc * rc
+                for lc, rc in zip(left_count, right_count))
+    s += 10 * partial + 45 * 10**(nofd//2) * total
+    s %= mod
+    left_sums, left_count = next(left_sums, left_count)
+    right_sums, right_count = next(right_sums, right_count)
+
+print(s)
diff --git a/231_prime_factorisation_of_binomial_coefficients.py b/231_prime_factorisation_of_binomial_coefficients.py
new file mode 100644
index 0000000..6822de3
--- /dev/null
+++ b/231_prime_factorisation_of_binomial_coefficients.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+from math import ceil
+
+def reducetwo(n):
+    c = 0
+    while not n%2:
+        n //= 2
+        c += 1
+    return n, c
+
+def primes(n):
+    sieve = [0] * (n//2)
+    for p in range(3, n, 2):
+        if sieve[p//2] > 0:
+            continue
+        for k in range(1, n//p+1, 2):
+            sieve[(p*k)//2] += p + sieve[k//2] - sieve[(p*k)//2]
+    return sieve
+
+
+n = 2*10**7
+k = 15*10**6
+sieve = primes(n)
+print("Primed.")
+s = 0
+
+for m in range(n, n-k, -1):
+    if not m % 2:
+        m, c = reducetwo(m)
+        s += 2*c
+    s += sieve[m//2]
+
+print("Calculated.")
+
+for m in range(2, k+1):
+    if not m % 2:
+        m, c = reducetwo(m)
+        s -= 2*c
+    s -= sieve[m//2]
+
+print(s)
diff --git a/243_resilience.py b/243_resilience.py
new file mode 100644
index 0000000..e1a64d7
--- /dev/null
+++ b/243_resilience.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+
+from itertools import combinations
+from operator import mul
+from functools import reduce
+
+def factorize(n, primes):
+    return [p for p in primes if not n%p]
+
+def totient(factors, n):
+    r = 0
+    for i in range(1, len(factors)+1):
+        for f in combinations(factors, i):
+            r += (-1)**i * (n // reduce(mul, f))
+    return n + r
+
+def get_primes(n):
+    sieve = [False] * n
+    r = []
+    for p in range(2, n):
+        if sieve[p]:
+            continue
+        r.append(p)
+        for i in range(p**2, n, p):
+            sieve[i] = True
+    return r
+
+a, b = 15499, 94744
+primes = get_primes(100)
+n = 1
+for i, p in enumerate(primes):
+    n *= p
+    t = totient(primes[:i+1], n)
+    if t*b < (n-1)*a:
+        break
+else:
+    print("Not found")
+    quit()
+
+n //= primes[i]
+factors = primes[:i]
+
+for i in range(2, 100):
+    t = totient(set(factors + factorize(i, primes)), n*i)
+    if t*b < (n*i - 1) * a:
+        print(n*i)
+        quit()
+print("Not found")
diff --git a/315_digital_clock_roots.py b/315_digital_clock_roots.py
new file mode 100644
index 0000000..d3c92ff
--- /dev/null
+++ b/315_digital_clock_roots.py
@@ -0,0 +1,58 @@
+#!/usr/bin/env python3
+
+from lib import primegen
+from itertools import product
+
+segments = {"th" : (0,1), "mh" : (1,1),
+            "bh" : (2,1), "tl" : (1,0),
+            "bl" : (2,0), "tr" : (1,2),
+            "br" : (2,2)
+           }
+
+def rootdiff(n):
+    chain = list(rootchain(n))
+    return sum(transdiff(n1, n2)
+               for n1, n2 in zip(chain, chain[1:]))
+
+def rootchain(n):
+    n = str(n)
+    while len(n) > 1:
+        yield n
+        n = str(sum(map(int, n)))
+    yield n
+
+def transdiff(n1, n2):
+    return sum(cache[(d1, d2)]
+               for d1, d2 in zip(n1[::-1], n2[::-1]))
+
+def init_cache(digs):
+    return {(d1, d2) : 2 * len(digs[d1] & digs[d2])
+            for d1, d2 in product(digs, repeat=2)
+           }
+
+def parsenum(text):
+    return set(s for s, (rw, cl) in segments.items()
+                 if str(text[rw] + 3*' ')[cl] != ' '
+              )
+
+def parse(fname):
+    with open(fname) as f:
+        r = f.read().split('\n')
+    for n in range(10):
+        yield str(n), parsenum(r[3*n : 3*n+3])
+
+def main(A, B):
+    s = 0
+    gen = primegen(B)
+    for p in gen:
+        if p > A: break
+    s = rootdiff(p)
+    for p in gen:
+        s += rootdiff(p)
+    return s
+
+digs = {d : segs for d, segs in parse("data/digs.txt")}
+digs[' '] = set()
+cache = init_cache(digs)
+
+print(main(10**7, 2*10**7))
diff --git a/329_prime_frog.py b/329_prime_frog.py
new file mode 100644
index 0000000..46e9912
--- /dev/null
+++ b/329_prime_frog.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+from lib import primegen
+from fractions import Fraction
+
+def is_prime(n):
+    return n in primes
+
+def Pprob(n):
+    if is_prime(n):
+        return Fraction(2, 3)
+    return Fraction(1, 3)
+
+def croak(n):
+    return {'P' : Pprob(n),
+            'N' : 1-Pprob(n)}
+
+def jump(sqrs, let):
+    r = [0] * len(sqrs)
+    for i in range(1, len(sqrs)-1):
+        prob = croak(i+1)[let] * \
+               Fraction(1, 2) * sqrs[i]
+        r[i+1] += prob
+        r[i-1] += prob
+    r[1] += sqrs[0] * croak(1)[let]
+    r[-2] += sqrs[-1] * croak(len(sqrs))[let]
+    return r
+
+def start(n):
+    return [Fraction(1, n)] * n
+
+
+seq = "PPPPNNPPPNPPNPN"
+N = 500
+
+primes = set(primegen(N))
+sqrs = start(N)
+
+for letter in seq:
+    sqrs = jump(sqrs, letter)
+
+print(sum(sqrs))
diff --git a/346_strong_repunits.py b/346_strong_repunits.py
new file mode 100644
index 0000000..d32d9f7
--- /dev/null
+++ b/346_strong_repunits.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+
+bound = 10**12
+nums = set()
+nums.add(1)
+
+for base in range(2, int(bound**0.5)+1):
+    r = base**2 + base + 1
+    pw = 3
+    while r < bound:
+        if r not in nums:
+            nums.add(r)
+        r += base**pw
+        pw += 1
+
+print(sum(nums))
diff --git a/347_largest_integer.py b/347_largest_integer.py
new file mode 100644
index 0000000..9ca0a11
--- /dev/null
+++ b/347_largest_integer.py
@@ -0,0 +1,29 @@
+#!/usr/bin/env python3
+
+from lib import primegen
+from math import log
+
+def maxpw(n, logN):
+    return int(logN / log(n))
+
+def M(p, q, logN):
+    return max(q**pw * p**maxpw(p, logN - pw*log(q))
+               for pw in range(1, maxpw(q, logN - log(p))+1)
+              )
+
+N = 10**7
+logN = log(N)
+
+primes = list(primegen(N//2+1))
+
+s = 0
+
+for i, p in enumerate(primes):
+    if p > int(N**0.5):
+        break
+    for q in primes[i+1:]:
+        if p*q > N:
+            break
+        s += M(p, q, logN)
+
+print(s)
diff --git a/357_prime_generating_integers.py b/357_prime_generating_integers.py
new file mode 100644
index 0000000..d2fba4d
--- /dev/null
+++ b/357_prime_generating_integers.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+
+from itertools import combinations
+from operator import mul
+
+def candidates(sieve):
+    n = len(sieve)
+    for i in range(4, n, 4):
+        sieve[i] = True
+    for p in range(3, n, 2):
+        if sieve[p]:
+            continue
+        for i in range(p**2, n, 2*p):
+            sieve[i] = p
+
+        for i in range(2*p**2, n, 2*p**2):
+            sieve[i] = True
+        if not sieve[p-1]:
+            yield p-1
+
+def is_prime(n):
+    return not sieve[n] and n%2
+
+def factors(n):
+    r = []
+    while not n%2:
+        r.append(2)
+        n //= 2
+    while sieve[n] != False:
+        r.append(sieve[n])
+        n //= sieve[n]
+    return r + [n]
+
+def test_alt(n):
+    fctrs = [f for f in factors(n)]
+    for i in range(len(fctrs)//2):
+        for c in combinations(fctrs, i+1):
+            d = reduce(mul, c)
+            if not is_prime(d + n//d):
+                return False
+    return True
+
+def test(n, fctrs, p=1, i=0, depth=1):
+    if depth > len(fctrs) // 2:
+        return True
+    for j, f in enumerate(fctrs[i:]):
+        if not is_prime(p*f + n//(p*f)) or not test(n, fctrs, p*f, j+1, depth+1):
+            return False
+    return True
+
+N = 10**8 + 2
+sieve = [False] * N
+
+print(sum(c for c in candidates(sieve) if test(c, factors(c))) + 1)
diff --git a/381_prime-k_factorial.py b/381_prime-k_factorial.py
new file mode 100644
index 0000000..83d9616
--- /dev/null
+++ b/381_prime-k_factorial.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+
+from lib import primegen
+
+def S(p):
+    return (9 * pow(-24, p-2, p)) % p
+
+N = 10**8
+
+print(sum(S(p) for p in list(primegen(N))[2:]))
diff --git a/387_harshad_numbers.py b/387_harshad_numbers.py
new file mode 100644
index 0000000..cc60e50
--- /dev/null
+++ b/387_harshad_numbers.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+from random import randint
+
+DEPTH = 10
+
+def is_prime(n):
+    if n == 2:
+        return True
+    if not n % 2 or n < 2:
+        return False
+    for i in range(DEPTH):
+        a = randint(2, n-1)
+        if pow(a, n-1, n) != 1:
+            return False
+    return True
+
+def harshad_nums(N, n=0, s=0):
+    r = []
+    for i in range(n==0, 10):
+        m = 10*n + i
+        nws = s + i
+        if m > N:
+            break
+        if m % nws:
+            continue
+        r.append((m, nws))
+        r += harshad_nums(N, m, nws)
+    return r
+
+def main(N):
+    for n, s in harshad_nums(N // 10):
+        if not is_prime(n // s):
+            continue
+        for d in range(1, 10):
+            m = 10*n + d
+            if m >= N:
+                break
+            if is_prime(m):
+                yield m
+
+N = 10**14
+
+print(sum(main(N)))
diff --git a/429_sum_of_squares.py b/429_sum_of_squares.py
new file mode 100644
index 0000000..c7352c6
--- /dev/null
+++ b/429_sum_of_squares.py
@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+
+from lib import primegen
+
+def fpw(p, n):
+    pw = 0
+    pwofp = p
+    while pwofp <= n:
+        pw += n // pwofp
+        pwofp *= p
+    return pw
+
+N = 10**8
+mod = 10**9 + 9
+
+k = 1
+for p in primegen(N+1):
+    k *= pow(p, 2 * fpw(p, N), mod) + 1
+    k %= mod
+
+print(k)
diff --git a/493_under_the_rainbow.py b/493_under_the_rainbow.py
new file mode 100644
index 0000000..40fa0c4
--- /dev/null
+++ b/493_under_the_rainbow.py
@@ -0,0 +1,20 @@
+#!/usr/bin/env python3
+
+from fractions import Fraction
+
+def pick(state, depth, numofc, cache={}):
+    if sum(state) == depth:
+        return len(state) - state.count(0)
+    if state in cache:
+        return cache[state]
+    r = 0
+    for c, num in enumerate(state):
+        nws = sorted(state[:c] + (num + 1,) + state[c+1:])
+        nws = tuple(nws)
+        r += Fraction(numofc-num, numofc*len(state) - sum(state)) * \
+             pick(nws, depth, numofc, cache)
+    cache[state] = r
+    return r
+
+r = pick(7*(0,), 20, 10)
+print(float(round(r, 9)))
diff --git a/504_square_on_the_inside.py b/504_square_on_the_inside.py
new file mode 100644
index 0000000..7f102d8
--- /dev/null
+++ b/504_square_on_the_inside.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+
+from fractions import gcd
+from itertools import combinations_with_replacement as combs
+from math import sqrt
+
+
+def points(a, b, c, d):
+    return (gcds[a][b] + gcds[b][c] + gcds[c][d] + gcds[d][a])//2 + 1
+
+def is_square(n):
+    return sieve[n]
+
+def refs(a,b,c,d):
+    r = [ (a,b,c,d),
+          (d,a,b,c),
+          (c,d,a,b),
+          (b,c,d,a),
+          (a,d,c,b),
+          (b,a,d,c),
+          (c,b,a,d),
+          (d,c,b,a)]
+    return list(set(r))
+
+count = 0
+m = 100
+space = range(1, m+1)
+
+gcds = [None for n in range(m+1)]
+for a in range(1, m+1):
+    gcds[a] = [0] + [a*b - gcd(a,b) for b in range(1, m+1)]
+
+sieve = [False] * 2*m**2
+n = 1
+while n**2 < len(sieve):
+    sieve[n**2] = True
+    n += 1
+
+for a,b,c,d in combs(space, 4):
+    perms = [(a,b,c,d), (a,c,b,d), (a,c,d,b)]
+    wins = set()
+    for p in perms:
+        if not is_square(points(*p)):
+            continue
+        wins |= set(refs(*p))
+    count += len(wins)
+
+print(count)
diff --git a/518_prime_triples.py b/518_prime_triples.py
new file mode 100644
index 0000000..6b7182d
--- /dev/null
+++ b/518_prime_triples.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+
+def getprimes(n):
+    sieve = [True] * n
+    for p in range(3, n, 2):
+        if not sieve[p]:
+            continue
+        for i in range(p**2, n, 2*p):
+            sieve[i] = False
+    return sieve
+
+def is_prime(n):
+    if n == 2:
+        return True
+    return n % 2 and sieve[n]
+
+N = 10**8
+sieve = getprimes(N)
+primes = [2] + [n for n in range(3, N, 2) if sieve[n]]
+
+for n in range(2, int(N**0.5)+1):
+    b = n**2
+    if n % 2:
+        b *= 2
+    for i in range(b, N, b):
+        sieve[i] = n
+
+s = 0
+
+for p in primes:
+    d = sieve[p+1]
+    k = d + 1
+    while ((p+1) * k**2)/d**2 <= N:
+        if is_prime(((p+1)*k)//d - 1) and is_prime(((p+1)*k**2)//d**2 - 1):
+            s += p + ((p+1)*k)//d - 1 + ((p+1)*k**2)//d**2 - 1
+        k += 1
+
+print(s)
diff --git a/549_divisibility_of_factorials.py b/549_divisibility_of_factorials.py
new file mode 100644
index 0000000..4f91282
--- /dev/null
+++ b/549_divisibility_of_factorials.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+
+def modpw(n, p):
+    r = 0
+    while n >= p:
+        r += n // p
+        n //= p
+    return r
+
+n = 10**8 + 1
+sieve = [0] * n
+
+for p in range(2, n):
+    if sieve[p] > 0:
+        continue
+    minN = p
+    frompw = 1
+    topw   = 2
+    k = 1
+    while k*p < n:
+        for pw in range(frompw, topw):
+            k *= p
+            for m in range(k, n, k):
+                sieve[m] = max(sieve[m], minN)
+        minN += p
+        frompw = topw
+        topw   = modpw(minN, p) + 1
+
+print(sum(sieve[2:]))
+
diff --git a/data/digs.txt b/data/digs.txt
new file mode 100644
index 0000000..cab633c
--- /dev/null
+++ b/data/digs.txt
@@ -0,0 +1,30 @@
+ _
+| |
+|_|
+
+  |
+  |
+ _
+ _|
+|_
+ _
+ _|
+ _|
+
+|_|
+  |
+ _
+|_
+ _|
+ _
+|_
+|_|
+ _
+| |
+  |
+ _
+|_|
+|_|
+ _
+|_|
+ _|
diff --git a/lib.py b/lib.py
new file mode 100644
index 0000000..8bc7234
--- /dev/null
+++ b/lib.py
@@ -0,0 +1,68 @@
+#!/usr/bin/env python3
+
+from random import randint
+from collections import Counter
+
+FACTOR_SIEVE = []
+
+def primegen(n):
+    sieve = [False] * n
+    yield 2
+    for p in range(3, n, 2):
+        if sieve[p]:
+            continue
+        for i in range(p**2, n, 2*p):
+            sieve[i] = True
+        yield p
+
+def iter_factors(n):
+    sieve = init_factors(n)
+    for i in range(n):
+        yield pfactors(i, sieve)
+
+def init_factors(n):
+    sieve = [False] * n
+    for p in range(3, n, 2):
+        if sieve[p]:
+            continue
+        for i in range(p**2, n, 2*p):
+            sieve[i] = p
+    return sieve
+
+def pfactors(n, sieve=FACTOR_SIEVE):
+    if n < 1:
+        if not n:
+            return Counter()
+        n = abs(n)
+    fctrs = Counter()
+    while not n % 2:
+        fctrs[2] += 1
+        n //= 2
+    while sieve[n]:
+        fctrs[sieve[n]] += 1
+        n //= sieve[n]
+    if n != 1:
+        fctrs[n] += 1
+    return fctrs
+
+
+def miller_rabin(n, k=10):
+    if not (n % 2 and n % 3):
+        return False
+    d = n - 1
+    r = 0
+    while not r % 2:
+        d //= 2
+        r += 1
+    for i in range(k):
+        a = randint(2, n-2)
+        x = pow(a, d, n)
+        if x == 1 or x == n - 1:
+            continue
+        for j in range(r-1):
+            x = pow(x, 2, n)
+            if x == n - 1:
+                break
+        else:
+            return False
+    return True