roll.py (ver5)

# There is a simple connection between natural numbers and DFAs on a singleton
# alphabet. The importance of this connection is easily illustrated in modular
# arithmetic.

# Def: A roll is a pair of natural numbers, the second one of which is nonzero.

# (Thus a roll is a fraction, but we won't make this connection.)

# To every roll [p, q] we associate a DFA which goes through p states and then
# goes into a loop of length q. Addition and multiplication of rolls, as well
# as of their states, is easily defined. The states of this DFA are
# conventionally taken to be the natural range [0, p+q).

# Good picture: https://pixabay.com/p-147631/?no_redirect

# (Actually, the theory of rolls should require no prior theory of natural
# numbers if a roll is defined more like a DFA directly, but that's a more
# beautiful story.)

# Here's some arithmetic aided by rolls.

import fractions
import collections
import functools
import operator

""" Compute what state the DFA of b will be in after natural input a. """
def nat_mod_roll(a, b):
    flat, loop = b
    return a if a < flat else (a - flat) % loop + flat

def prod(a, one=1):
    return functools.reduce(operator.mul, a, one)

"""
Compute many factorizations, primes, totient factorizations and totients.

The return value can be used as the 'cache' parameter when calling functions
that accept a cache. The size n of the cache must be strictly greater than the
cyclic part of the modulus or these functions may not function.
"""
def make_cache(n):
    # factorize all numbers below n
    facts = [collections.Counter() for _ in range(n)]
    for i in range(2, n):
        if not facts[i]:  # i is a prime; mark its multiples
            k = i
            while k < n:
                for j in range(k, n, k):
                    facts[j][i] += 1
                k *= i
    # find all primes below n
    primes = [i for i, j in enumerate(facts) if sum(j.values())==1]
    # compute euler_phi of all numbers below n, in factored form
    phi_facts = [sum((sum((facts[i] for _ in range(j-1)),
                          collections.Counter())+facts[i-1]
        for i, j in c.items()), collections.Counter()) for c in facts]
    # compute euler_phi of all numbers below n
    phis = [prod(j**k for j, k in i.items()) for i in phi_facts]
    phis[0] = 0
    return facts, primes, phi_facts, phis

"""
Compute the smallest roll whose DFA states can index the states that the DFA
of b ends up in after any input a^i.
"""
def roll_of_powers_small(a, b):
    powers = [1]
    while True:
        c = nat_mod_roll(a*powers[-1], b)
        if c in powers:
            return powers.index(c), len(powers) - powers.index(c)
        powers.append(c)

"""
Compute a small roll whose DFA states can index the states that the DFA of b
ends up in after any input a^i.
"""
def roll_of_powers_big(a, b, cache):
    if a == 0: return 1, 1
    if a == 1: return 0, 1
    m = next(iter(i for i in range(b[0]+b[1]) if a**i >= b[0]))
    n = 0
    d = b[1] // fractions.gcd(b[1], a**m % b[1])
    c = fractions.gcd(d, a % b[1])
    while c != 1:
        d = d // c
        c = fractions.gcd(d, c)
        n += 1
    _, _, _, phi = cache
    return m + n, phi[d]

def roll_of_powers(a, b, cache=None):
    if cache:
        return roll_of_powers_big(a, b, cache)
    else:
        return roll_of_powers_small(a, b)

""" Compute what state the DFA of b will be in after tower input a. """
def tower_mod_roll(a, b, cache=None):
    if len(a) == 0:
        return nat_mod_roll(1, b)
    return nat_mod_roll(a[0]**tower_mod_roll(a[1:], roll_of_powers(a[0], b, cache), cache), b)

def mod(a, b, cache=None):
    if isinstance(a, int): a = a,
    if isinstance(b, int): b = 0, b
    return tower_mod_roll(a, b, cache)

def test():
    arg_a = [2, 10], [100, 1000, 10000], [104, 2], [2, 3]
    arg_b = 100, 100000, 5, 112
    expect = 24, 0, 1, 8
    for a, b, good in zip(arg_a, arg_b, expect):
        c = mod(a, b)
        print("mod({}, {}) = {} -- ".format(a, b, c), end="")
        print("OK" if c == good else "FAIL (expected {})".format(good))

if __name__ == "__main__":
    test()