Python script to find numerological coincidences in bible verses

numerology.py

# Find numerological coincidences in bible verses
# https://skeptics.stackexchange.com/questions/53659/are-the-mathematical-constants-e-and-%cf%80-encoded-in-the-bible

import itertools
import functools
import math
import random

# Source data, in this case Genesis 1:1
words = [[400, 10, 300, 1, 200, 2], [1, 200, 2], [40, 10, 5, 30, 1], [400, 1], [40, 10, 40, 300, 5], [400, 1, 6], [90, 200, 1, 5]]

# Helpers, just all letters flattened and a function to calculate product of list
flat = sum(words, [])
product = lambda x: functools.reduce(lambda a,b: a * b, x)

# List of constants, operators and terms that we use to find a match
constants = [math.pi, math.tau, math.e, math.sqrt(2), ( 1 + math.sqrt(5) ) / 2]
operators = ['+', '-', '*', '/']
terms = {
    'num_words': len(words),
    'num_letters': len(flat),
    'sum_letters': sum(flat),
    'prod_letters': product(flat),
    'sum_word_prods': sum([product(x) for x in words]),
    'prod_word_sums': product([sum(x) for x in words]),
    'sum_word_lens': sum([len(x) for x in words]),
    'prod_word_lens': product([len(x) for x in words]),
}

def evaluate_stack(expr):
    '''Evaluate an expression in RPN notation'''
    stack = []
    for term in expr:
        if term in '*/+-':
            b = stack.pop()
            a = stack.pop()
            
            if term == '*': stack.append(a * b)
            elif term == '/': stack.append(a / b)
            elif term == '+': stack.append(a + b)
            elif term == '-': stack.append(a - b)
        elif term in terms:
            stack.append(terms[term])
        else:
            stack.append(float(term))

    # Normalize to range 1 to 10
    result = abs(stack[-1])
    if result > 0:
        result *= 10**(-math.floor(math.log10(result)))
    return result

print("Terms:")
for k, v in terms.items():
    print("%s = %d" % (k, v))
print("\n")

# Now just randomly search through combinations
best = {}
termnames = list(terms.keys())
choices = operators + termnames
while True:
    # Form a random RPN expression
    expr = []
    expr.append(random.choice(termnames))
    expr.append(random.choice(termnames))
    for i in range(random.randint(1, 6)):
        expr.append(random.choice(choices))
    
    try:
        result = evaluate_stack(expr)
    except IndexError:
        continue
    except ZeroDivisionError:
        continue
    
    for c in constants:
        rel_error = round(abs(result - c) / c, 9)
        if rel_error < best.get(c, 0.0001):
            best[c] = rel_error
            print("%s = %0.5f ~ %0.5f, relative error %0.8f" % (str(expr), result, c, rel_error))