modpl_generator.py

#
# Manifestation of Permuted Diagonal Line generator
#
import numpy as np
from PIL import Image

def digital_reverse(n, length, base):
  r = 0
  for _ in range(length):
    r = base*r + n % base
    n /= base
  return r

assert digital_reverse(123, 3, 10) == 321
assert digital_reverse(54321, 5, 10) == 12345
assert digital_reverse(0x123, 3, 16) == 0x321
assert digital_reverse(0x54321, 5, 16) == 0x12345

def save_matrix(m, fn, modulo):
    image = Image.new("L", (m.shape[1], m.shape[0]), 0)
    pix = image.load()
    for y in xrange(m.shape[0]):
        for x in xrange(m.shape[1]):
             pix[x, y] = 255*m[y, x]/(modulo-1)

    image.save(fn, "png")

if __name__ == '__main__':
    # state size
    modulo = 16
    power = 2

    print "modulo&base =", modulo

    half = modulo**power
    size = 2*half
    print "state size", size

    # known state after applying rule 15
    reversible_state = np.zeros((size, size), np.int16)
    for i in xrange(half):
        r = digital_reverse(i, power, modulo)
        reversible_state[i, r] = 1
        reversible_state[i, r+half] = modulo-1
        reversible_state[i+half, r] = modulo-1
        reversible_state[i+half, r+half] = 1

    # partially known state
    reversed_state = np.empty((size, size), np.int16)

    # fill initial (known) values
    reversed_state[0,:] = 0 # first row
    reversed_state[:,0] = 0 # first column

    # find rest values of the state
    for i0 in xrange(size-1):
        i1 = i0 + 1
        for j0 in xrange(size-1):
            j1 = j0 + 1
            # Quadruple rule layout:
            #      j0 | j1
            #    +----+----+
            # i0 | r0 | r2 |
            # ---+----+----+
            # i1 | r8 | r4 |
            #    +----+----+
            # solve r4[t-1] from (r1[t-1] - r2[t-1] + r4[t-1] - r8[t-1]) == r1[t] (mod modulo)
            # where [t-1] is related to reversed_state, [t] is related to reversible_state
            reversed_state[i1, j1] = (reversible_state[i0, j0] - reversed_state[i0, j0] + reversed_state[i0, j1] + reversed_state[i1, j0]) % modulo

    # produce image
    save_matrix(reversed_state, "modpl-modulo%d-power%d.png"%(modulo, power), modulo)