Prime Christmas Tree Generator

tree.py

#!/usr/bin/python3
# Bart Massey

# Prime Christmas Tree printer, inspired by
# https://www.reddit.com/r/interestingasfuck/
#  comments/efhy5m/i_found_this_really_great_the_correct_alignment/

# This code is licensed under the "MIT license".  See
# https://opensource.org/licenses/MIT for license terms.

import random
import sys

# Miller-Rabin test below is recursive, and this number is
# going to be big. Should probably re-implement Miller-Rabin
# iteratively.
sys.setrecursionlimit(10000)

# Template Christmas tree from Reddit post.  Put . where
# digits are to be filled in.
tree = """
20191111111111111111111111111111111111
111111111111111111..111111111111111111
11111111111111111....11111111111111111
1111111111111111......1111111111111111
111111111111111........111111111111111
11111111111111...........1111111111111
111111111111..............111111111111
111111111111111........111111111111111
11111111111111..........11111111111111
1111111111111............1111111111111
11111111111................11111111111
1111111111..................1111111111
11111111......................11111111
111111..........................111111
1111111111111............1111111111111
11111111111...............111111111111
1111111111..................1111111111
11111111......................11111111
111111..........................111111
1111..............................1111
11..................................11
11111111111111111....11111111111111111
11111111111111111....11111111111111111
11111111111111111....11111111111111111
"""

# Miller-Rabin probabilistic primality test.
# Code based on a 1999 Nickle implementation by me.
def is_composite(n, d):
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
    
    def witness_exp(b, e, m):
        if e == 0:
            return (0, 1)
        if e == 1:
            return (b % m, 0)
        p, w = witness_exp(b, e // 2, m)
        if w != 0:
            return (p, w)
        t = (p ** 2) % m
        if t == 1 and p != 1 and p != m - 1:
            return (t, p)
        if e % 2 == 0:
            return (t, w)
        return ((t * b) % m, w)

    def witness(a, n):
        p, w = witness_exp(a, n - 1, n)
        if w != 0:
            return True
        if p != 1:
            return True
        return False

    for p in primes:
        if n % p == 0:
            return True
    for _ in range(d):
        a = 1 + random.randrange(n - 1)
        if witness(a, n):
            return True
    return False

# Repeatedly fill in the . characters in the template with
# random digits and see if the resulting number is prime.
tree_digits = ''.join(tree.split('\n'))
while True:
    td = list(tree_digits)
    for i, c in enumerate(td):
        if c == '.':
            td[i] = chr(ord('0') + random.randrange(10))
    tn = ''.join(td)
    if not is_composite(int(tn), 40):
        break

# Substitute the . characters in the original tree template
# to produce a prime-tree.
ti = 0
tl = list(tree)
for i, c in enumerate(tl):
    if c == '\n':
        continue
    if c.isnumeric():
        assert tn[ti] == c
        ti += 1
        continue
    assert c == '.'
    tl[i] = tn[ti]
    ti += 1

# Render the tree.
print(''.join(tl))